1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 11 PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian) 12 { 13 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 14 PetscInt j, k, n = A->rmap->n; 15 PetscScalar *v; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 21 if (!hermitian) { 22 for (k=0;k<n;k++) { 23 for (j=k;j<n;j++) { 24 v[j*mat->lda + k] = v[k*mat->lda + j]; 25 } 26 } 27 } else { 28 for (k=0;k<n;k++) { 29 for (j=k;j<n;j++) { 30 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 31 } 32 } 33 } 34 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 39 { 40 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 41 PetscErrorCode ierr; 42 PetscBLASInt info,n; 43 44 PetscFunctionBegin; 45 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 46 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 47 if (A->factortype == MAT_FACTOR_LU) { 48 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 49 if (!mat->fwork) { 50 mat->lfwork = n; 51 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 52 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 53 } 54 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 55 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 56 ierr = PetscFPTrapPop();CHKERRQ(ierr); 57 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 58 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 59 if (A->spd) { 60 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 61 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 62 ierr = PetscFPTrapPop();CHKERRQ(ierr); 63 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 64 #if defined(PETSC_USE_COMPLEX) 65 } else if (A->hermitian) { 66 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 67 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 70 ierr = PetscFPTrapPop();CHKERRQ(ierr); 71 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 72 #endif 73 } else { /* symmetric case */ 74 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 75 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 76 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 77 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 78 ierr = PetscFPTrapPop();CHKERRQ(ierr); 79 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr); 80 } 81 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1); 82 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 84 85 A->ops->solve = NULL; 86 A->ops->matsolve = NULL; 87 A->ops->solvetranspose = NULL; 88 A->ops->matsolvetranspose = NULL; 89 A->ops->solveadd = NULL; 90 A->ops->solvetransposeadd = NULL; 91 A->factortype = MAT_FACTOR_NONE; 92 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 97 { 98 PetscErrorCode ierr; 99 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 100 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 101 PetscScalar *slot,*bb,*v; 102 const PetscScalar *xx; 103 104 PetscFunctionBegin; 105 if (PetscDefined(USE_DEBUG)) { 106 for (i=0; i<N; i++) { 107 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 108 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 109 if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n); 110 } 111 } 112 if (!N) PetscFunctionReturn(0); 113 114 /* fix right hand side if needed */ 115 if (x && b) { 116 Vec xt; 117 118 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 119 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr); 120 ierr = VecCopy(x,xt);CHKERRQ(ierr); 121 ierr = VecScale(xt,-1.0);CHKERRQ(ierr); 122 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr); 123 ierr = VecDestroy(&xt);CHKERRQ(ierr); 124 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 125 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 126 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 127 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 128 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 129 } 130 131 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 132 for (i=0; i<N; i++) { 133 slot = v + rows[i]*m; 134 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 135 } 136 for (i=0; i<N; i++) { 137 slot = v + rows[i]; 138 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 139 } 140 if (diag != 0.0) { 141 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 142 for (i=0; i<N; i++) { 143 slot = v + (m+1)*rows[i]; 144 *slot = diag; 145 } 146 } 147 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 152 { 153 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 if (c->ptapwork) { 158 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr); 159 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr); 160 } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first"); 161 PetscFunctionReturn(0); 162 } 163 164 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C) 165 { 166 Mat_SeqDense *c; 167 PetscBool cisdense; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 ierr = MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);CHKERRQ(ierr); 172 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 173 if (!cisdense) { 174 PetscBool flg; 175 176 ierr = PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 177 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 178 } 179 ierr = MatSetUp(C);CHKERRQ(ierr); 180 c = (Mat_SeqDense*)C->data; 181 ierr = MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);CHKERRQ(ierr); 182 ierr = MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);CHKERRQ(ierr); 183 ierr = MatSetType(c->ptapwork,((PetscObject)C)->type_name);CHKERRQ(ierr); 184 ierr = MatSetUp(c->ptapwork);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 188 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 189 { 190 Mat B = NULL; 191 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 192 Mat_SeqDense *b; 193 PetscErrorCode ierr; 194 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 195 MatScalar *av=a->a; 196 PetscBool isseqdense; 197 198 PetscFunctionBegin; 199 if (reuse == MAT_REUSE_MATRIX) { 200 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 201 if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 202 } 203 if (reuse != MAT_REUSE_MATRIX) { 204 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 205 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 206 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 207 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr); 208 b = (Mat_SeqDense*)(B->data); 209 } else { 210 b = (Mat_SeqDense*)((*newmat)->data); 211 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr); 212 } 213 for (i=0; i<m; i++) { 214 PetscInt j; 215 for (j=0;j<ai[1]-ai[0];j++) { 216 b->v[*aj*m+i] = *av; 217 aj++; 218 av++; 219 } 220 ai++; 221 } 222 223 if (reuse == MAT_INPLACE_MATRIX) { 224 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 225 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 226 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 227 } else { 228 if (B) *newmat = B; 229 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 230 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 231 } 232 PetscFunctionReturn(0); 233 } 234 235 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 236 { 237 Mat B; 238 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 239 PetscErrorCode ierr; 240 PetscInt i, j; 241 PetscInt *rows, *nnz; 242 MatScalar *aa = a->v, *vals; 243 244 PetscFunctionBegin; 245 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 246 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 247 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 248 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 249 for (j=0; j<A->cmap->n; j++) { 250 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i]; 251 aa += a->lda; 252 } 253 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 254 aa = a->v; 255 for (j=0; j<A->cmap->n; j++) { 256 PetscInt numRows = 0; 257 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];} 258 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 259 aa += a->lda; 260 } 261 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 262 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 263 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 264 265 if (reuse == MAT_INPLACE_MATRIX) { 266 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 267 } else { 268 *newmat = B; 269 } 270 PetscFunctionReturn(0); 271 } 272 273 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 274 { 275 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 276 const PetscScalar *xv; 277 PetscScalar *yv; 278 PetscBLASInt N,m,ldax,lday,one = 1; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 ierr = MatDenseGetArrayRead(X,&xv);CHKERRQ(ierr); 283 ierr = MatDenseGetArray(Y,&yv);CHKERRQ(ierr); 284 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 285 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 286 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 287 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 288 if (ldax>m || lday>m) { 289 PetscInt j; 290 291 for (j=0; j<X->cmap->n; j++) { 292 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one)); 293 } 294 } else { 295 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one)); 296 } 297 ierr = MatDenseRestoreArrayRead(X,&xv);CHKERRQ(ierr); 298 ierr = MatDenseRestoreArray(Y,&yv);CHKERRQ(ierr); 299 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 304 { 305 PetscLogDouble N = A->rmap->n*A->cmap->n; 306 307 PetscFunctionBegin; 308 info->block_size = 1.0; 309 info->nz_allocated = N; 310 info->nz_used = N; 311 info->nz_unneeded = 0; 312 info->assemblies = A->num_ass; 313 info->mallocs = 0; 314 info->memory = ((PetscObject)A)->mem; 315 info->fill_ratio_given = 0; 316 info->fill_ratio_needed = 0; 317 info->factor_mallocs = 0; 318 PetscFunctionReturn(0); 319 } 320 321 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 322 { 323 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 324 PetscScalar *v; 325 PetscErrorCode ierr; 326 PetscBLASInt one = 1,j,nz,lda; 327 328 PetscFunctionBegin; 329 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 330 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 331 if (lda>A->rmap->n) { 332 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 333 for (j=0; j<A->cmap->n; j++) { 334 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one)); 335 } 336 } else { 337 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 338 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one)); 339 } 340 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 341 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 342 PetscFunctionReturn(0); 343 } 344 345 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 346 { 347 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 348 PetscInt i,j,m = A->rmap->n,N = a->lda; 349 const PetscScalar *v; 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 *fl = PETSC_FALSE; 354 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 355 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 356 for (i=0; i<m; i++) { 357 for (j=i; j<m; j++) { 358 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) { 359 goto restore; 360 } 361 } 362 } 363 *fl = PETSC_TRUE; 364 restore: 365 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 366 PetscFunctionReturn(0); 367 } 368 369 static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 370 { 371 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 372 PetscInt i,j,m = A->rmap->n,N = a->lda; 373 const PetscScalar *v; 374 PetscErrorCode ierr; 375 376 PetscFunctionBegin; 377 *fl = PETSC_FALSE; 378 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 379 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 380 for (i=0; i<m; i++) { 381 for (j=i; j<m; j++) { 382 if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) { 383 goto restore; 384 } 385 } 386 } 387 *fl = PETSC_TRUE; 388 restore: 389 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 394 { 395 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 396 PetscErrorCode ierr; 397 PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda; 398 399 PetscFunctionBegin; 400 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 401 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 402 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */ 403 ierr = MatDenseSetLDA(newi,lda);CHKERRQ(ierr); 404 } 405 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 406 if (cpvalues == MAT_COPY_VALUES) { 407 const PetscScalar *av; 408 PetscScalar *v; 409 410 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 411 ierr = MatDenseGetArray(newi,&v);CHKERRQ(ierr); 412 ierr = MatDenseGetLDA(newi,&nlda);CHKERRQ(ierr); 413 m = A->rmap->n; 414 if (lda>m || nlda>m) { 415 for (j=0; j<A->cmap->n; j++) { 416 ierr = PetscArraycpy(v+j*nlda,av+j*lda,m);CHKERRQ(ierr); 417 } 418 } else { 419 ierr = PetscArraycpy(v,av,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 420 } 421 ierr = MatDenseRestoreArray(newi,&v);CHKERRQ(ierr); 422 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 423 } 424 PetscFunctionReturn(0); 425 } 426 427 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 428 { 429 PetscErrorCode ierr; 430 431 PetscFunctionBegin; 432 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 433 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 434 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 435 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 436 PetscFunctionReturn(0); 437 } 438 439 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 440 { 441 MatFactorInfo info; 442 PetscErrorCode ierr; 443 444 PetscFunctionBegin; 445 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 446 ierr = (*fact->ops->lufactor)(fact,0,0,&info);CHKERRQ(ierr); 447 PetscFunctionReturn(0); 448 } 449 450 static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 451 { 452 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 453 PetscErrorCode ierr; 454 const PetscScalar *x; 455 PetscScalar *y; 456 PetscBLASInt one = 1,info,m; 457 458 PetscFunctionBegin; 459 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 460 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 461 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 462 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 463 if (A->factortype == MAT_FACTOR_LU) { 464 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 465 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 466 ierr = PetscFPTrapPop();CHKERRQ(ierr); 467 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 468 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 469 if (A->spd) { 470 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 471 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 472 ierr = PetscFPTrapPop();CHKERRQ(ierr); 473 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 474 #if defined(PETSC_USE_COMPLEX) 475 } else if (A->hermitian) { 476 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 477 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 478 ierr = PetscFPTrapPop();CHKERRQ(ierr); 479 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 480 #endif 481 } else { /* symmetric case */ 482 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 483 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 484 ierr = PetscFPTrapPop();CHKERRQ(ierr); 485 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 486 } 487 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 488 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 489 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 490 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 491 PetscFunctionReturn(0); 492 } 493 494 static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 495 { 496 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 497 PetscErrorCode ierr; 498 const PetscScalar *b; 499 PetscScalar *x; 500 PetscInt n; 501 PetscBLASInt nrhs,info,m; 502 503 PetscFunctionBegin; 504 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 505 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 506 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 507 ierr = MatDenseGetArrayRead(B,&b);CHKERRQ(ierr); 508 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 509 510 ierr = PetscArraycpy(x,b,m*nrhs);CHKERRQ(ierr); 511 512 if (A->factortype == MAT_FACTOR_LU) { 513 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 514 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 515 ierr = PetscFPTrapPop();CHKERRQ(ierr); 516 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 517 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 518 if (A->spd) { 519 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 520 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 521 ierr = PetscFPTrapPop();CHKERRQ(ierr); 522 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 523 #if defined(PETSC_USE_COMPLEX) 524 } else if (A->hermitian) { 525 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 526 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 527 ierr = PetscFPTrapPop();CHKERRQ(ierr); 528 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 529 #endif 530 } else { /* symmetric case */ 531 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 532 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 533 ierr = PetscFPTrapPop();CHKERRQ(ierr); 534 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 535 } 536 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 537 538 ierr = MatDenseRestoreArrayRead(B,&b);CHKERRQ(ierr); 539 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 540 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 541 PetscFunctionReturn(0); 542 } 543 544 static PetscErrorCode MatConjugate_SeqDense(Mat); 545 546 static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 547 { 548 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 549 PetscErrorCode ierr; 550 const PetscScalar *x; 551 PetscScalar *y; 552 PetscBLASInt one = 1,info,m; 553 554 PetscFunctionBegin; 555 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 556 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 557 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 558 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 559 if (A->factortype == MAT_FACTOR_LU) { 560 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 561 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 562 ierr = PetscFPTrapPop();CHKERRQ(ierr); 563 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 564 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 565 if (A->spd) { 566 #if defined(PETSC_USE_COMPLEX) 567 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 568 #endif 569 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 570 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 571 ierr = PetscFPTrapPop();CHKERRQ(ierr); 572 #if defined(PETSC_USE_COMPLEX) 573 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 574 #endif 575 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 576 #if defined(PETSC_USE_COMPLEX) 577 } else if (A->hermitian) { 578 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 579 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 580 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 581 ierr = PetscFPTrapPop();CHKERRQ(ierr); 582 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 583 #endif 584 } else { /* symmetric case */ 585 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 586 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 587 ierr = PetscFPTrapPop();CHKERRQ(ierr); 588 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 589 } 590 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 591 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 592 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 593 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 594 PetscFunctionReturn(0); 595 } 596 597 /* ---------------------------------------------------------------*/ 598 /* COMMENT: I have chosen to hide row permutation in the pivots, 599 rather than put it in the Mat->row slot.*/ 600 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 601 { 602 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 603 PetscErrorCode ierr; 604 PetscBLASInt n,m,info; 605 606 PetscFunctionBegin; 607 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 608 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 609 if (!mat->pivots) { 610 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 611 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 612 } 613 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 614 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 615 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 616 ierr = PetscFPTrapPop();CHKERRQ(ierr); 617 618 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 619 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 620 621 A->ops->solve = MatSolve_SeqDense; 622 A->ops->matsolve = MatMatSolve_SeqDense; 623 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 624 A->factortype = MAT_FACTOR_LU; 625 626 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 627 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 628 629 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 630 PetscFunctionReturn(0); 631 } 632 633 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 634 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 635 { 636 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 637 PetscErrorCode ierr; 638 PetscBLASInt info,n; 639 640 PetscFunctionBegin; 641 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 642 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 643 if (A->spd) { 644 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 645 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 646 ierr = PetscFPTrapPop();CHKERRQ(ierr); 647 #if defined(PETSC_USE_COMPLEX) 648 } else if (A->hermitian) { 649 if (!mat->pivots) { 650 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 651 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 652 } 653 if (!mat->fwork) { 654 PetscScalar dummy; 655 656 mat->lfwork = -1; 657 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 658 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 659 ierr = PetscFPTrapPop();CHKERRQ(ierr); 660 mat->lfwork = (PetscInt)PetscRealPart(dummy); 661 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 662 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 663 } 664 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 665 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 666 ierr = PetscFPTrapPop();CHKERRQ(ierr); 667 #endif 668 } else { /* symmetric case */ 669 if (!mat->pivots) { 670 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 671 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 672 } 673 if (!mat->fwork) { 674 PetscScalar dummy; 675 676 mat->lfwork = -1; 677 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 678 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 679 ierr = PetscFPTrapPop();CHKERRQ(ierr); 680 mat->lfwork = (PetscInt)PetscRealPart(dummy); 681 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 682 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 683 } 684 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 685 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 686 ierr = PetscFPTrapPop();CHKERRQ(ierr); 687 } 688 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 689 690 A->ops->solve = MatSolve_SeqDense; 691 A->ops->matsolve = MatMatSolve_SeqDense; 692 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 693 A->factortype = MAT_FACTOR_CHOLESKY; 694 695 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 696 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 697 698 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 699 PetscFunctionReturn(0); 700 } 701 702 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 703 { 704 PetscErrorCode ierr; 705 MatFactorInfo info; 706 707 PetscFunctionBegin; 708 info.fill = 1.0; 709 710 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 711 ierr = (*fact->ops->choleskyfactor)(fact,0,&info);CHKERRQ(ierr); 712 PetscFunctionReturn(0); 713 } 714 715 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 716 { 717 PetscFunctionBegin; 718 fact->assembled = PETSC_TRUE; 719 fact->preallocated = PETSC_TRUE; 720 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 721 fact->ops->solve = MatSolve_SeqDense; 722 fact->ops->matsolve = MatMatSolve_SeqDense; 723 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 724 PetscFunctionReturn(0); 725 } 726 727 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 728 { 729 PetscFunctionBegin; 730 fact->preallocated = PETSC_TRUE; 731 fact->assembled = PETSC_TRUE; 732 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 733 fact->ops->solve = MatSolve_SeqDense; 734 fact->ops->matsolve = MatMatSolve_SeqDense; 735 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 736 PetscFunctionReturn(0); 737 } 738 739 /* uses LAPACK */ 740 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 741 { 742 PetscErrorCode ierr; 743 744 PetscFunctionBegin; 745 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 746 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 747 ierr = MatSetType(*fact,MATDENSE);CHKERRQ(ierr); 748 if (ftype == MAT_FACTOR_LU) { 749 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 750 } else { 751 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 752 } 753 (*fact)->factortype = ftype; 754 755 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 756 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 /* ------------------------------------------------------------------*/ 761 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 762 { 763 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 764 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 765 const PetscScalar *b; 766 PetscErrorCode ierr; 767 PetscInt m = A->rmap->n,i; 768 PetscBLASInt o = 1,bm; 769 770 PetscFunctionBegin; 771 #if defined(PETSC_HAVE_CUDA) 772 if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented"); 773 #endif 774 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 775 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 776 if (flag & SOR_ZERO_INITIAL_GUESS) { 777 /* this is a hack fix, should have another version without the second BLASdotu */ 778 ierr = VecSet(xx,zero);CHKERRQ(ierr); 779 } 780 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 781 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 782 its = its*lits; 783 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 784 while (its--) { 785 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 786 for (i=0; i<m; i++) { 787 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 788 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 789 } 790 } 791 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 792 for (i=m-1; i>=0; i--) { 793 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 794 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 795 } 796 } 797 } 798 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 799 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 800 PetscFunctionReturn(0); 801 } 802 803 /* -----------------------------------------------------------------*/ 804 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 805 { 806 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 807 const PetscScalar *v = mat->v,*x; 808 PetscScalar *y; 809 PetscErrorCode ierr; 810 PetscBLASInt m, n,_One=1; 811 PetscScalar _DOne=1.0,_DZero=0.0; 812 813 PetscFunctionBegin; 814 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 815 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 816 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 817 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 818 if (!A->rmap->n || !A->cmap->n) { 819 PetscBLASInt i; 820 for (i=0; i<n; i++) y[i] = 0.0; 821 } else { 822 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 823 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 824 } 825 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 826 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 827 PetscFunctionReturn(0); 828 } 829 830 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 831 { 832 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 833 PetscScalar *y,_DOne=1.0,_DZero=0.0; 834 PetscErrorCode ierr; 835 PetscBLASInt m, n, _One=1; 836 const PetscScalar *v = mat->v,*x; 837 838 PetscFunctionBegin; 839 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 840 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 841 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 842 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 843 if (!A->rmap->n || !A->cmap->n) { 844 PetscBLASInt i; 845 for (i=0; i<m; i++) y[i] = 0.0; 846 } else { 847 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 848 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 849 } 850 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 851 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854 855 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 856 { 857 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 858 const PetscScalar *v = mat->v,*x; 859 PetscScalar *y,_DOne=1.0; 860 PetscErrorCode ierr; 861 PetscBLASInt m, n, _One=1; 862 863 PetscFunctionBegin; 864 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 865 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 866 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 867 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 868 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 869 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 870 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 871 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 872 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 873 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 874 PetscFunctionReturn(0); 875 } 876 877 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 878 { 879 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 880 const PetscScalar *v = mat->v,*x; 881 PetscScalar *y; 882 PetscErrorCode ierr; 883 PetscBLASInt m, n, _One=1; 884 PetscScalar _DOne=1.0; 885 886 PetscFunctionBegin; 887 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 888 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 889 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 890 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 891 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 892 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 893 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 894 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 895 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 896 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 897 PetscFunctionReturn(0); 898 } 899 900 /* -----------------------------------------------------------------*/ 901 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 902 { 903 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 904 PetscErrorCode ierr; 905 PetscInt i; 906 907 PetscFunctionBegin; 908 *ncols = A->cmap->n; 909 if (cols) { 910 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 911 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 912 } 913 if (vals) { 914 const PetscScalar *v; 915 916 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 917 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 918 v += row; 919 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 920 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 921 } 922 PetscFunctionReturn(0); 923 } 924 925 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 926 { 927 PetscErrorCode ierr; 928 929 PetscFunctionBegin; 930 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 931 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 932 PetscFunctionReturn(0); 933 } 934 /* ----------------------------------------------------------------*/ 935 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 936 { 937 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 938 PetscScalar *av; 939 PetscInt i,j,idx=0; 940 #if defined(PETSC_HAVE_CUDA) 941 PetscOffloadMask oldf; 942 #endif 943 PetscErrorCode ierr; 944 945 PetscFunctionBegin; 946 ierr = MatDenseGetArray(A,&av);CHKERRQ(ierr); 947 if (!mat->roworiented) { 948 if (addv == INSERT_VALUES) { 949 for (j=0; j<n; j++) { 950 if (indexn[j] < 0) {idx += m; continue;} 951 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 952 for (i=0; i<m; i++) { 953 if (indexm[i] < 0) {idx++; continue;} 954 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 955 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 956 } 957 } 958 } else { 959 for (j=0; j<n; j++) { 960 if (indexn[j] < 0) {idx += m; continue;} 961 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 962 for (i=0; i<m; i++) { 963 if (indexm[i] < 0) {idx++; continue;} 964 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 965 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 966 } 967 } 968 } 969 } else { 970 if (addv == INSERT_VALUES) { 971 for (i=0; i<m; i++) { 972 if (indexm[i] < 0) { idx += n; continue;} 973 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 974 for (j=0; j<n; j++) { 975 if (indexn[j] < 0) { idx++; continue;} 976 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 977 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 978 } 979 } 980 } else { 981 for (i=0; i<m; i++) { 982 if (indexm[i] < 0) { idx += n; continue;} 983 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 984 for (j=0; j<n; j++) { 985 if (indexn[j] < 0) { idx++; continue;} 986 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 987 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 988 } 989 } 990 } 991 } 992 /* hack to prevent unneeded copy to the GPU while returning the array */ 993 #if defined(PETSC_HAVE_CUDA) 994 oldf = A->offloadmask; 995 A->offloadmask = PETSC_OFFLOAD_GPU; 996 #endif 997 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr); 998 #if defined(PETSC_HAVE_CUDA) 999 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 1000 #endif 1001 PetscFunctionReturn(0); 1002 } 1003 1004 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 1005 { 1006 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1007 const PetscScalar *vv; 1008 PetscInt i,j; 1009 PetscErrorCode ierr; 1010 1011 PetscFunctionBegin; 1012 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1013 /* row-oriented output */ 1014 for (i=0; i<m; i++) { 1015 if (indexm[i] < 0) {v += n;continue;} 1016 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 1017 for (j=0; j<n; j++) { 1018 if (indexn[j] < 0) {v++; continue;} 1019 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 1020 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1021 } 1022 } 1023 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1024 PetscFunctionReturn(0); 1025 } 1026 1027 /* -----------------------------------------------------------------*/ 1028 1029 PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer) 1030 { 1031 PetscErrorCode ierr; 1032 PetscBool skipHeader; 1033 PetscViewerFormat format; 1034 PetscInt header[4],M,N,m,lda,i,j,k; 1035 const PetscScalar *v; 1036 PetscScalar *vwork; 1037 1038 PetscFunctionBegin; 1039 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1040 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1041 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1042 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1043 1044 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1045 1046 /* write matrix header */ 1047 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1048 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1049 if (!skipHeader) {ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr);} 1050 1051 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1052 if (format != PETSC_VIEWER_NATIVE) { 1053 PetscInt nnz = m*N, *iwork; 1054 /* store row lengths for each row */ 1055 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr); 1056 for (i=0; i<m; i++) iwork[i] = N; 1057 ierr = PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1058 /* store column indices (zero start index) */ 1059 for (k=0, i=0; i<m; i++) 1060 for (j=0; j<N; j++, k++) 1061 iwork[k] = j; 1062 ierr = PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1063 ierr = PetscFree(iwork);CHKERRQ(ierr); 1064 } 1065 /* store matrix values as a dense matrix in row major order */ 1066 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr); 1067 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr); 1068 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1069 for (k=0, i=0; i<m; i++) 1070 for (j=0; j<N; j++, k++) 1071 vwork[k] = v[i+lda*j]; 1072 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr); 1073 ierr = PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1074 ierr = PetscFree(vwork);CHKERRQ(ierr); 1075 PetscFunctionReturn(0); 1076 } 1077 1078 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer) 1079 { 1080 PetscErrorCode ierr; 1081 PetscBool skipHeader; 1082 PetscInt header[4],M,N,m,nz,lda,i,j,k; 1083 PetscInt rows,cols; 1084 PetscScalar *v,*vwork; 1085 1086 PetscFunctionBegin; 1087 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1088 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1089 1090 if (!skipHeader) { 1091 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1092 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1093 M = header[1]; N = header[2]; 1094 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 1095 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 1096 nz = header[3]; 1097 if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz); 1098 } else { 1099 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1100 if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix"); 1101 nz = MATRIX_BINARY_FORMAT_DENSE; 1102 } 1103 1104 /* setup global sizes if not set */ 1105 if (mat->rmap->N < 0) mat->rmap->N = M; 1106 if (mat->cmap->N < 0) mat->cmap->N = N; 1107 ierr = MatSetUp(mat);CHKERRQ(ierr); 1108 /* check if global sizes are correct */ 1109 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1110 if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 1111 1112 ierr = MatGetSize(mat,NULL,&N);CHKERRQ(ierr); 1113 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1114 ierr = MatDenseGetArray(mat,&v);CHKERRQ(ierr); 1115 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1116 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1117 PetscInt nnz = m*N; 1118 /* read in matrix values */ 1119 ierr = PetscMalloc1(nnz,&vwork);CHKERRQ(ierr); 1120 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1121 /* store values in column major order */ 1122 for (j=0; j<N; j++) 1123 for (i=0; i<m; i++) 1124 v[i+lda*j] = vwork[i*N+j]; 1125 ierr = PetscFree(vwork);CHKERRQ(ierr); 1126 } else { /* matrix in file is sparse format */ 1127 PetscInt nnz = 0, *rlens, *icols; 1128 /* read in row lengths */ 1129 ierr = PetscMalloc1(m,&rlens);CHKERRQ(ierr); 1130 ierr = PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1131 for (i=0; i<m; i++) nnz += rlens[i]; 1132 /* read in column indices and values */ 1133 ierr = PetscMalloc2(nnz,&icols,nnz,&vwork);CHKERRQ(ierr); 1134 ierr = PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1135 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1136 /* store values in column major order */ 1137 for (k=0, i=0; i<m; i++) 1138 for (j=0; j<rlens[i]; j++, k++) 1139 v[i+lda*icols[k]] = vwork[k]; 1140 ierr = PetscFree(rlens);CHKERRQ(ierr); 1141 ierr = PetscFree2(icols,vwork);CHKERRQ(ierr); 1142 } 1143 ierr = MatDenseRestoreArray(mat,&v);CHKERRQ(ierr); 1144 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1145 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1146 PetscFunctionReturn(0); 1147 } 1148 1149 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1150 { 1151 PetscBool isbinary, ishdf5; 1152 PetscErrorCode ierr; 1153 1154 PetscFunctionBegin; 1155 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1156 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1157 /* force binary viewer to load .info file if it has not yet done so */ 1158 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1159 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1161 if (isbinary) { 1162 ierr = MatLoad_Dense_Binary(newMat,viewer);CHKERRQ(ierr); 1163 } else if (ishdf5) { 1164 #if defined(PETSC_HAVE_HDF5) 1165 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1166 #else 1167 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1168 #endif 1169 } else { 1170 SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 1171 } 1172 PetscFunctionReturn(0); 1173 } 1174 1175 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1176 { 1177 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1178 PetscErrorCode ierr; 1179 PetscInt i,j; 1180 const char *name; 1181 PetscScalar *v,*av; 1182 PetscViewerFormat format; 1183 #if defined(PETSC_USE_COMPLEX) 1184 PetscBool allreal = PETSC_TRUE; 1185 #endif 1186 1187 PetscFunctionBegin; 1188 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1189 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1190 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1191 PetscFunctionReturn(0); /* do nothing for now */ 1192 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1193 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1194 for (i=0; i<A->rmap->n; i++) { 1195 v = av + i; 1196 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1197 for (j=0; j<A->cmap->n; j++) { 1198 #if defined(PETSC_USE_COMPLEX) 1199 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1200 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1201 } else if (PetscRealPart(*v)) { 1202 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1203 } 1204 #else 1205 if (*v) { 1206 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1207 } 1208 #endif 1209 v += a->lda; 1210 } 1211 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1212 } 1213 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1214 } else { 1215 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1216 #if defined(PETSC_USE_COMPLEX) 1217 /* determine if matrix has all real values */ 1218 v = av; 1219 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1220 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1221 } 1222 #endif 1223 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1224 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1225 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1226 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1228 } 1229 1230 for (i=0; i<A->rmap->n; i++) { 1231 v = av + i; 1232 for (j=0; j<A->cmap->n; j++) { 1233 #if defined(PETSC_USE_COMPLEX) 1234 if (allreal) { 1235 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1236 } else { 1237 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1238 } 1239 #else 1240 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1241 #endif 1242 v += a->lda; 1243 } 1244 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1245 } 1246 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1247 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1248 } 1249 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1250 } 1251 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1252 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1253 PetscFunctionReturn(0); 1254 } 1255 1256 #include <petscdraw.h> 1257 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1258 { 1259 Mat A = (Mat) Aa; 1260 PetscErrorCode ierr; 1261 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1262 int color = PETSC_DRAW_WHITE; 1263 const PetscScalar *v; 1264 PetscViewer viewer; 1265 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1266 PetscViewerFormat format; 1267 1268 PetscFunctionBegin; 1269 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1270 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1271 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1272 1273 /* Loop over matrix elements drawing boxes */ 1274 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 1275 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1276 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1277 /* Blue for negative and Red for positive */ 1278 for (j = 0; j < n; j++) { 1279 x_l = j; x_r = x_l + 1.0; 1280 for (i = 0; i < m; i++) { 1281 y_l = m - i - 1.0; 1282 y_r = y_l + 1.0; 1283 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1284 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1285 else continue; 1286 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1287 } 1288 } 1289 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1290 } else { 1291 /* use contour shading to indicate magnitude of values */ 1292 /* first determine max of all nonzero values */ 1293 PetscReal minv = 0.0, maxv = 0.0; 1294 PetscDraw popup; 1295 1296 for (i=0; i < m*n; i++) { 1297 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1298 } 1299 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1300 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1301 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1302 1303 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1304 for (j=0; j<n; j++) { 1305 x_l = j; 1306 x_r = x_l + 1.0; 1307 for (i=0; i<m; i++) { 1308 y_l = m - i - 1.0; 1309 y_r = y_l + 1.0; 1310 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1311 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1312 } 1313 } 1314 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1315 } 1316 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 1317 PetscFunctionReturn(0); 1318 } 1319 1320 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1321 { 1322 PetscDraw draw; 1323 PetscBool isnull; 1324 PetscReal xr,yr,xl,yl,h,w; 1325 PetscErrorCode ierr; 1326 1327 PetscFunctionBegin; 1328 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1329 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1330 if (isnull) PetscFunctionReturn(0); 1331 1332 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1333 xr += w; yr += h; xl = -w; yl = -h; 1334 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1335 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1336 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1337 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1338 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1339 PetscFunctionReturn(0); 1340 } 1341 1342 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1343 { 1344 PetscErrorCode ierr; 1345 PetscBool iascii,isbinary,isdraw; 1346 1347 PetscFunctionBegin; 1348 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1349 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1350 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1351 1352 if (iascii) { 1353 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1354 } else if (isbinary) { 1355 ierr = MatView_Dense_Binary(A,viewer);CHKERRQ(ierr); 1356 } else if (isdraw) { 1357 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1358 } 1359 PetscFunctionReturn(0); 1360 } 1361 1362 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array) 1363 { 1364 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1365 1366 PetscFunctionBegin; 1367 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1368 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1369 if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first"); 1370 a->unplacedarray = a->v; 1371 a->unplaced_user_alloc = a->user_alloc; 1372 a->v = (PetscScalar*) array; 1373 a->user_alloc = PETSC_TRUE; 1374 #if defined(PETSC_HAVE_CUDA) 1375 A->offloadmask = PETSC_OFFLOAD_CPU; 1376 #endif 1377 PetscFunctionReturn(0); 1378 } 1379 1380 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1381 { 1382 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1383 1384 PetscFunctionBegin; 1385 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1386 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1387 a->v = a->unplacedarray; 1388 a->user_alloc = a->unplaced_user_alloc; 1389 a->unplacedarray = NULL; 1390 #if defined(PETSC_HAVE_CUDA) 1391 A->offloadmask = PETSC_OFFLOAD_CPU; 1392 #endif 1393 PetscFunctionReturn(0); 1394 } 1395 1396 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array) 1397 { 1398 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1399 PetscErrorCode ierr; 1400 1401 PetscFunctionBegin; 1402 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1403 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1404 if (!a->user_alloc) { ierr = PetscFree(a->v);CHKERRQ(ierr); } 1405 a->v = (PetscScalar*) array; 1406 a->user_alloc = PETSC_FALSE; 1407 #if defined(PETSC_HAVE_CUDA) 1408 A->offloadmask = PETSC_OFFLOAD_CPU; 1409 #endif 1410 PetscFunctionReturn(0); 1411 } 1412 1413 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1414 { 1415 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1416 PetscErrorCode ierr; 1417 1418 PetscFunctionBegin; 1419 #if defined(PETSC_USE_LOG) 1420 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1421 #endif 1422 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1423 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1424 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1425 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1426 if (!l->unplaced_user_alloc) {ierr = PetscFree(l->unplacedarray);CHKERRQ(ierr);} 1427 if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1428 if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1429 ierr = VecDestroy(&l->cvec);CHKERRQ(ierr); 1430 ierr = MatDestroy(&l->cmat);CHKERRQ(ierr); 1431 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1432 1433 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1434 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1435 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);CHKERRQ(ierr); 1436 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1437 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1440 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);CHKERRQ(ierr); 1441 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr); 1442 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr); 1443 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);CHKERRQ(ierr); 1444 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);CHKERRQ(ierr); 1445 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1446 #if defined(PETSC_HAVE_ELEMENTAL) 1447 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1448 #endif 1449 #if defined(PETSC_HAVE_CUDA) 1450 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1451 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr); 1452 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr); 1453 #endif 1454 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1455 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1456 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr); 1457 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1458 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1459 1460 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1461 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);CHKERRQ(ierr); 1464 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);CHKERRQ(ierr); 1469 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);CHKERRQ(ierr); 1470 PetscFunctionReturn(0); 1471 } 1472 1473 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1474 { 1475 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1476 PetscErrorCode ierr; 1477 PetscInt k,j,m,n,M; 1478 PetscScalar *v,tmp; 1479 1480 PetscFunctionBegin; 1481 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1482 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1483 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1484 for (j=0; j<m; j++) { 1485 for (k=0; k<j; k++) { 1486 tmp = v[j + k*M]; 1487 v[j + k*M] = v[k + j*M]; 1488 v[k + j*M] = tmp; 1489 } 1490 } 1491 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1492 } else { /* out-of-place transpose */ 1493 Mat tmat; 1494 Mat_SeqDense *tmatd; 1495 PetscScalar *v2; 1496 PetscInt M2; 1497 1498 if (reuse != MAT_REUSE_MATRIX) { 1499 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1500 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1501 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1502 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1503 } else tmat = *matout; 1504 1505 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1506 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1507 tmatd = (Mat_SeqDense*)tmat->data; 1508 M2 = tmatd->lda; 1509 for (j=0; j<n; j++) { 1510 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1511 } 1512 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1513 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1514 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1515 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1516 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1517 else { 1518 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1519 } 1520 } 1521 PetscFunctionReturn(0); 1522 } 1523 1524 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1525 { 1526 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1527 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1528 PetscInt i; 1529 const PetscScalar *v1,*v2; 1530 PetscErrorCode ierr; 1531 1532 PetscFunctionBegin; 1533 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1534 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1535 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1536 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1537 for (i=0; i<A1->cmap->n; i++) { 1538 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1539 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1540 v1 += mat1->lda; 1541 v2 += mat2->lda; 1542 } 1543 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1544 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1545 *flg = PETSC_TRUE; 1546 PetscFunctionReturn(0); 1547 } 1548 1549 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1550 { 1551 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1552 PetscInt i,n,len; 1553 PetscScalar *x; 1554 const PetscScalar *vv; 1555 PetscErrorCode ierr; 1556 1557 PetscFunctionBegin; 1558 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1559 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1560 len = PetscMin(A->rmap->n,A->cmap->n); 1561 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1562 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1563 for (i=0; i<len; i++) { 1564 x[i] = vv[i*mat->lda + i]; 1565 } 1566 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1567 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1568 PetscFunctionReturn(0); 1569 } 1570 1571 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1572 { 1573 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1574 const PetscScalar *l,*r; 1575 PetscScalar x,*v,*vv; 1576 PetscErrorCode ierr; 1577 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1578 1579 PetscFunctionBegin; 1580 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1581 if (ll) { 1582 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1583 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1584 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1585 for (i=0; i<m; i++) { 1586 x = l[i]; 1587 v = vv + i; 1588 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1589 } 1590 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1591 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1592 } 1593 if (rr) { 1594 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1595 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1596 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1597 for (i=0; i<n; i++) { 1598 x = r[i]; 1599 v = vv + i*mat->lda; 1600 for (j=0; j<m; j++) (*v++) *= x; 1601 } 1602 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1603 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1604 } 1605 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1606 PetscFunctionReturn(0); 1607 } 1608 1609 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1610 { 1611 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1612 PetscScalar *v,*vv; 1613 PetscReal sum = 0.0; 1614 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1615 PetscErrorCode ierr; 1616 1617 PetscFunctionBegin; 1618 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1619 v = vv; 1620 if (type == NORM_FROBENIUS) { 1621 if (lda>m) { 1622 for (j=0; j<A->cmap->n; j++) { 1623 v = vv+j*lda; 1624 for (i=0; i<m; i++) { 1625 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1626 } 1627 } 1628 } else { 1629 #if defined(PETSC_USE_REAL___FP16) 1630 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1631 *nrm = BLASnrm2_(&cnt,v,&one); 1632 } 1633 #else 1634 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1635 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1636 } 1637 } 1638 *nrm = PetscSqrtReal(sum); 1639 #endif 1640 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1641 } else if (type == NORM_1) { 1642 *nrm = 0.0; 1643 for (j=0; j<A->cmap->n; j++) { 1644 v = vv + j*mat->lda; 1645 sum = 0.0; 1646 for (i=0; i<A->rmap->n; i++) { 1647 sum += PetscAbsScalar(*v); v++; 1648 } 1649 if (sum > *nrm) *nrm = sum; 1650 } 1651 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1652 } else if (type == NORM_INFINITY) { 1653 *nrm = 0.0; 1654 for (j=0; j<A->rmap->n; j++) { 1655 v = vv + j; 1656 sum = 0.0; 1657 for (i=0; i<A->cmap->n; i++) { 1658 sum += PetscAbsScalar(*v); v += mat->lda; 1659 } 1660 if (sum > *nrm) *nrm = sum; 1661 } 1662 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1663 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1664 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1665 PetscFunctionReturn(0); 1666 } 1667 1668 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1669 { 1670 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1671 PetscErrorCode ierr; 1672 1673 PetscFunctionBegin; 1674 switch (op) { 1675 case MAT_ROW_ORIENTED: 1676 aij->roworiented = flg; 1677 break; 1678 case MAT_NEW_NONZERO_LOCATIONS: 1679 case MAT_NEW_NONZERO_LOCATION_ERR: 1680 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1681 case MAT_NEW_DIAGONALS: 1682 case MAT_KEEP_NONZERO_PATTERN: 1683 case MAT_IGNORE_OFF_PROC_ENTRIES: 1684 case MAT_USE_HASH_TABLE: 1685 case MAT_IGNORE_ZERO_ENTRIES: 1686 case MAT_IGNORE_LOWER_TRIANGULAR: 1687 case MAT_SORTED_FULL: 1688 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1689 break; 1690 case MAT_SPD: 1691 case MAT_SYMMETRIC: 1692 case MAT_STRUCTURALLY_SYMMETRIC: 1693 case MAT_HERMITIAN: 1694 case MAT_SYMMETRY_ETERNAL: 1695 /* These options are handled directly by MatSetOption() */ 1696 break; 1697 default: 1698 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1699 } 1700 PetscFunctionReturn(0); 1701 } 1702 1703 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1704 { 1705 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1706 PetscErrorCode ierr; 1707 PetscInt lda=l->lda,m=A->rmap->n,j; 1708 PetscScalar *v; 1709 1710 PetscFunctionBegin; 1711 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1712 if (lda>m) { 1713 for (j=0; j<A->cmap->n; j++) { 1714 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1715 } 1716 } else { 1717 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1718 } 1719 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1720 PetscFunctionReturn(0); 1721 } 1722 1723 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1724 { 1725 PetscErrorCode ierr; 1726 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1727 PetscInt m = l->lda, n = A->cmap->n, i,j; 1728 PetscScalar *slot,*bb,*v; 1729 const PetscScalar *xx; 1730 1731 PetscFunctionBegin; 1732 if (PetscDefined(USE_DEBUG)) { 1733 for (i=0; i<N; i++) { 1734 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1735 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 1736 } 1737 } 1738 if (!N) PetscFunctionReturn(0); 1739 1740 /* fix right hand side if needed */ 1741 if (x && b) { 1742 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1743 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1744 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1745 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1746 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1747 } 1748 1749 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1750 for (i=0; i<N; i++) { 1751 slot = v + rows[i]; 1752 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1753 } 1754 if (diag != 0.0) { 1755 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1756 for (i=0; i<N; i++) { 1757 slot = v + (m+1)*rows[i]; 1758 *slot = diag; 1759 } 1760 } 1761 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1762 PetscFunctionReturn(0); 1763 } 1764 1765 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1766 { 1767 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1768 1769 PetscFunctionBegin; 1770 *lda = mat->lda; 1771 PetscFunctionReturn(0); 1772 } 1773 1774 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 1775 { 1776 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1777 1778 PetscFunctionBegin; 1779 if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1780 *array = mat->v; 1781 PetscFunctionReturn(0); 1782 } 1783 1784 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 1785 { 1786 PetscFunctionBegin; 1787 *array = NULL; 1788 PetscFunctionReturn(0); 1789 } 1790 1791 /*@C 1792 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1793 1794 Not collective 1795 1796 Input Parameter: 1797 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1798 1799 Output Parameter: 1800 . lda - the leading dimension 1801 1802 Level: intermediate 1803 1804 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA() 1805 @*/ 1806 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1807 { 1808 PetscErrorCode ierr; 1809 1810 PetscFunctionBegin; 1811 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1812 PetscValidPointer(lda,2); 1813 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1814 PetscFunctionReturn(0); 1815 } 1816 1817 /*@C 1818 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 1819 1820 Not collective 1821 1822 Input Parameter: 1823 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1824 - lda - the leading dimension 1825 1826 Level: intermediate 1827 1828 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA() 1829 @*/ 1830 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 1831 { 1832 PetscErrorCode ierr; 1833 1834 PetscFunctionBegin; 1835 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1836 ierr = PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));CHKERRQ(ierr); 1837 PetscFunctionReturn(0); 1838 } 1839 1840 /*@C 1841 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 1842 1843 Logically Collective on Mat 1844 1845 Input Parameter: 1846 . mat - a dense matrix 1847 1848 Output Parameter: 1849 . array - pointer to the data 1850 1851 Level: intermediate 1852 1853 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1854 @*/ 1855 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1856 { 1857 PetscErrorCode ierr; 1858 1859 PetscFunctionBegin; 1860 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1861 PetscValidPointer(array,2); 1862 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1863 PetscFunctionReturn(0); 1864 } 1865 1866 /*@C 1867 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1868 1869 Logically Collective on Mat 1870 1871 Input Parameters: 1872 + mat - a dense matrix 1873 - array - pointer to the data 1874 1875 Level: intermediate 1876 1877 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1878 @*/ 1879 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1880 { 1881 PetscErrorCode ierr; 1882 1883 PetscFunctionBegin; 1884 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1885 PetscValidPointer(array,2); 1886 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1887 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1888 #if defined(PETSC_HAVE_CUDA) 1889 A->offloadmask = PETSC_OFFLOAD_CPU; 1890 #endif 1891 PetscFunctionReturn(0); 1892 } 1893 1894 /*@C 1895 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 1896 1897 Not Collective 1898 1899 Input Parameter: 1900 . mat - a dense matrix 1901 1902 Output Parameter: 1903 . array - pointer to the data 1904 1905 Level: intermediate 1906 1907 .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1908 @*/ 1909 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1910 { 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1915 PetscValidPointer(array,2); 1916 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1917 PetscFunctionReturn(0); 1918 } 1919 1920 /*@C 1921 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 1922 1923 Not Collective 1924 1925 Input Parameters: 1926 + mat - a dense matrix 1927 - array - pointer to the data 1928 1929 Level: intermediate 1930 1931 .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1932 @*/ 1933 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1934 { 1935 PetscErrorCode ierr; 1936 1937 PetscFunctionBegin; 1938 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1939 PetscValidPointer(array,2); 1940 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1941 PetscFunctionReturn(0); 1942 } 1943 1944 /*@C 1945 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 1946 1947 Not Collective 1948 1949 Input Parameter: 1950 . mat - a dense matrix 1951 1952 Output Parameter: 1953 . array - pointer to the data 1954 1955 Level: intermediate 1956 1957 .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1958 @*/ 1959 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 1960 { 1961 PetscErrorCode ierr; 1962 1963 PetscFunctionBegin; 1964 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1965 PetscValidPointer(array,2); 1966 ierr = PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1967 PetscFunctionReturn(0); 1968 } 1969 1970 /*@C 1971 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 1972 1973 Not Collective 1974 1975 Input Parameters: 1976 + mat - a dense matrix 1977 - array - pointer to the data 1978 1979 Level: intermediate 1980 1981 .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1982 @*/ 1983 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 1984 { 1985 PetscErrorCode ierr; 1986 1987 PetscFunctionBegin; 1988 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1989 PetscValidPointer(array,2); 1990 ierr = PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1991 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1992 #if defined(PETSC_HAVE_CUDA) 1993 A->offloadmask = PETSC_OFFLOAD_CPU; 1994 #endif 1995 PetscFunctionReturn(0); 1996 } 1997 1998 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1999 { 2000 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2001 PetscErrorCode ierr; 2002 PetscInt i,j,nrows,ncols,blda; 2003 const PetscInt *irow,*icol; 2004 PetscScalar *av,*bv,*v = mat->v; 2005 Mat newmat; 2006 2007 PetscFunctionBegin; 2008 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2009 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2010 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2011 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2012 2013 /* Check submatrixcall */ 2014 if (scall == MAT_REUSE_MATRIX) { 2015 PetscInt n_cols,n_rows; 2016 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2017 if (n_rows != nrows || n_cols != ncols) { 2018 /* resize the result matrix to match number of requested rows/columns */ 2019 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2020 } 2021 newmat = *B; 2022 } else { 2023 /* Create and fill new matrix */ 2024 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 2025 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2026 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 2027 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 2028 } 2029 2030 /* Now extract the data pointers and do the copy,column at a time */ 2031 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 2032 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 2033 for (i=0; i<ncols; i++) { 2034 av = v + mat->lda*icol[i]; 2035 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2036 bv += blda; 2037 } 2038 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 2039 2040 /* Assemble the matrices so that the correct flags are set */ 2041 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2042 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2043 2044 /* Free work space */ 2045 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2046 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2047 *B = newmat; 2048 PetscFunctionReturn(0); 2049 } 2050 2051 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2052 { 2053 PetscErrorCode ierr; 2054 PetscInt i; 2055 2056 PetscFunctionBegin; 2057 if (scall == MAT_INITIAL_MATRIX) { 2058 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2059 } 2060 2061 for (i=0; i<n; i++) { 2062 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2063 } 2064 PetscFunctionReturn(0); 2065 } 2066 2067 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2068 { 2069 PetscFunctionBegin; 2070 PetscFunctionReturn(0); 2071 } 2072 2073 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2074 { 2075 PetscFunctionBegin; 2076 PetscFunctionReturn(0); 2077 } 2078 2079 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2080 { 2081 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2082 PetscErrorCode ierr; 2083 const PetscScalar *va; 2084 PetscScalar *vb; 2085 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2086 2087 PetscFunctionBegin; 2088 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2089 if (A->ops->copy != B->ops->copy) { 2090 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2091 PetscFunctionReturn(0); 2092 } 2093 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2094 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2095 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2096 if (lda1>m || lda2>m) { 2097 for (j=0; j<n; j++) { 2098 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2099 } 2100 } else { 2101 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2102 } 2103 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2104 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2105 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2106 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2111 { 2112 PetscErrorCode ierr; 2113 2114 PetscFunctionBegin; 2115 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2116 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2117 if (!A->preallocated) { 2118 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2119 } 2120 PetscFunctionReturn(0); 2121 } 2122 2123 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2124 { 2125 PetscInt i,nz = A->rmap->n*A->cmap->n; 2126 PetscScalar *aa; 2127 PetscErrorCode ierr; 2128 2129 PetscFunctionBegin; 2130 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2131 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2132 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2133 PetscFunctionReturn(0); 2134 } 2135 2136 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2137 { 2138 PetscInt i,nz = A->rmap->n*A->cmap->n; 2139 PetscScalar *aa; 2140 PetscErrorCode ierr; 2141 2142 PetscFunctionBegin; 2143 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2144 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2145 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2146 PetscFunctionReturn(0); 2147 } 2148 2149 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2150 { 2151 PetscInt i,nz = A->rmap->n*A->cmap->n; 2152 PetscScalar *aa; 2153 PetscErrorCode ierr; 2154 2155 PetscFunctionBegin; 2156 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2157 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2158 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2159 PetscFunctionReturn(0); 2160 } 2161 2162 /* ----------------------------------------------------------------*/ 2163 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2164 { 2165 PetscErrorCode ierr; 2166 PetscInt m=A->rmap->n,n=B->cmap->n; 2167 PetscBool cisdense; 2168 2169 PetscFunctionBegin; 2170 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2171 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2172 if (!cisdense) { 2173 PetscBool flg; 2174 2175 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2176 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2177 } 2178 ierr = MatSetUp(C);CHKERRQ(ierr); 2179 PetscFunctionReturn(0); 2180 } 2181 2182 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2183 { 2184 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2185 PetscBLASInt m,n,k; 2186 const PetscScalar *av,*bv; 2187 PetscScalar *cv; 2188 PetscScalar _DOne=1.0,_DZero=0.0; 2189 PetscErrorCode ierr; 2190 2191 PetscFunctionBegin; 2192 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2193 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2194 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2195 if (!m || !n || !k) PetscFunctionReturn(0); 2196 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2197 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2198 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2199 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2200 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2201 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2202 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2203 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2208 { 2209 PetscErrorCode ierr; 2210 PetscInt m=A->rmap->n,n=B->rmap->n; 2211 PetscBool cisdense; 2212 2213 PetscFunctionBegin; 2214 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2215 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2216 if (!cisdense) { 2217 PetscBool flg; 2218 2219 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2220 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2221 } 2222 ierr = MatSetUp(C);CHKERRQ(ierr); 2223 PetscFunctionReturn(0); 2224 } 2225 2226 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2227 { 2228 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2229 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2230 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2231 const PetscScalar *av,*bv; 2232 PetscScalar *cv; 2233 PetscBLASInt m,n,k; 2234 PetscScalar _DOne=1.0,_DZero=0.0; 2235 PetscErrorCode ierr; 2236 2237 PetscFunctionBegin; 2238 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2239 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2240 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2241 if (!m || !n || !k) PetscFunctionReturn(0); 2242 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2243 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2244 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2245 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2246 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2247 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2248 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2249 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2250 PetscFunctionReturn(0); 2251 } 2252 2253 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2254 { 2255 PetscErrorCode ierr; 2256 PetscInt m=A->cmap->n,n=B->cmap->n; 2257 PetscBool cisdense; 2258 2259 PetscFunctionBegin; 2260 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2261 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2262 if (!cisdense) { 2263 PetscBool flg; 2264 2265 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2266 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2267 } 2268 ierr = MatSetUp(C);CHKERRQ(ierr); 2269 PetscFunctionReturn(0); 2270 } 2271 2272 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2273 { 2274 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2275 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2276 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2277 const PetscScalar *av,*bv; 2278 PetscScalar *cv; 2279 PetscBLASInt m,n,k; 2280 PetscScalar _DOne=1.0,_DZero=0.0; 2281 PetscErrorCode ierr; 2282 2283 PetscFunctionBegin; 2284 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2285 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2286 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2287 if (!m || !n || !k) PetscFunctionReturn(0); 2288 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2289 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2290 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2291 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2292 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2293 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2294 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2295 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2296 PetscFunctionReturn(0); 2297 } 2298 2299 /* ----------------------------------------------- */ 2300 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2301 { 2302 PetscFunctionBegin; 2303 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2304 C->ops->productsymbolic = MatProductSymbolic_AB; 2305 PetscFunctionReturn(0); 2306 } 2307 2308 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2309 { 2310 PetscFunctionBegin; 2311 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2312 C->ops->productsymbolic = MatProductSymbolic_AtB; 2313 PetscFunctionReturn(0); 2314 } 2315 2316 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2317 { 2318 PetscFunctionBegin; 2319 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2320 C->ops->productsymbolic = MatProductSymbolic_ABt; 2321 PetscFunctionReturn(0); 2322 } 2323 2324 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2325 { 2326 PetscErrorCode ierr; 2327 Mat_Product *product = C->product; 2328 2329 PetscFunctionBegin; 2330 switch (product->type) { 2331 case MATPRODUCT_AB: 2332 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2333 break; 2334 case MATPRODUCT_AtB: 2335 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2336 break; 2337 case MATPRODUCT_ABt: 2338 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2339 break; 2340 default: 2341 break; 2342 } 2343 PetscFunctionReturn(0); 2344 } 2345 /* ----------------------------------------------- */ 2346 2347 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2348 { 2349 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2350 PetscErrorCode ierr; 2351 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2352 PetscScalar *x; 2353 const PetscScalar *aa; 2354 2355 PetscFunctionBegin; 2356 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2357 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2358 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2359 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2360 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2361 for (i=0; i<m; i++) { 2362 x[i] = aa[i]; if (idx) idx[i] = 0; 2363 for (j=1; j<n; j++) { 2364 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2365 } 2366 } 2367 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2368 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2369 PetscFunctionReturn(0); 2370 } 2371 2372 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2373 { 2374 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2375 PetscErrorCode ierr; 2376 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2377 PetscScalar *x; 2378 PetscReal atmp; 2379 const PetscScalar *aa; 2380 2381 PetscFunctionBegin; 2382 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2383 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2384 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2385 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2386 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2387 for (i=0; i<m; i++) { 2388 x[i] = PetscAbsScalar(aa[i]); 2389 for (j=1; j<n; j++) { 2390 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2391 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2392 } 2393 } 2394 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2395 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2396 PetscFunctionReturn(0); 2397 } 2398 2399 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2400 { 2401 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2402 PetscErrorCode ierr; 2403 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2404 PetscScalar *x; 2405 const PetscScalar *aa; 2406 2407 PetscFunctionBegin; 2408 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2409 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2410 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2411 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2412 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2413 for (i=0; i<m; i++) { 2414 x[i] = aa[i]; if (idx) idx[i] = 0; 2415 for (j=1; j<n; j++) { 2416 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2417 } 2418 } 2419 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2420 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2421 PetscFunctionReturn(0); 2422 } 2423 2424 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2425 { 2426 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2427 PetscErrorCode ierr; 2428 PetscScalar *x; 2429 const PetscScalar *aa; 2430 2431 PetscFunctionBegin; 2432 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2433 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2434 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2435 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2436 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2437 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2438 PetscFunctionReturn(0); 2439 } 2440 2441 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2442 { 2443 PetscErrorCode ierr; 2444 PetscInt i,j,m,n; 2445 const PetscScalar *a; 2446 2447 PetscFunctionBegin; 2448 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2449 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2450 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2451 if (type == NORM_2) { 2452 for (i=0; i<n; i++) { 2453 for (j=0; j<m; j++) { 2454 norms[i] += PetscAbsScalar(a[j]*a[j]); 2455 } 2456 a += m; 2457 } 2458 } else if (type == NORM_1) { 2459 for (i=0; i<n; i++) { 2460 for (j=0; j<m; j++) { 2461 norms[i] += PetscAbsScalar(a[j]); 2462 } 2463 a += m; 2464 } 2465 } else if (type == NORM_INFINITY) { 2466 for (i=0; i<n; i++) { 2467 for (j=0; j<m; j++) { 2468 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2469 } 2470 a += m; 2471 } 2472 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2473 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2474 if (type == NORM_2) { 2475 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2476 } 2477 PetscFunctionReturn(0); 2478 } 2479 2480 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2481 { 2482 PetscErrorCode ierr; 2483 PetscScalar *a; 2484 PetscInt lda,m,n,i,j; 2485 2486 PetscFunctionBegin; 2487 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2488 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2489 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2490 for (j=0; j<n; j++) { 2491 for (i=0; i<m; i++) { 2492 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2493 } 2494 } 2495 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2496 PetscFunctionReturn(0); 2497 } 2498 2499 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2500 { 2501 PetscFunctionBegin; 2502 *missing = PETSC_FALSE; 2503 PetscFunctionReturn(0); 2504 } 2505 2506 /* vals is not const */ 2507 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2508 { 2509 PetscErrorCode ierr; 2510 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2511 PetscScalar *v; 2512 2513 PetscFunctionBegin; 2514 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2515 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2516 *vals = v+col*a->lda; 2517 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2518 PetscFunctionReturn(0); 2519 } 2520 2521 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2522 { 2523 PetscFunctionBegin; 2524 *vals = 0; /* user cannot accidently use the array later */ 2525 PetscFunctionReturn(0); 2526 } 2527 2528 /* -------------------------------------------------------------------*/ 2529 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2530 MatGetRow_SeqDense, 2531 MatRestoreRow_SeqDense, 2532 MatMult_SeqDense, 2533 /* 4*/ MatMultAdd_SeqDense, 2534 MatMultTranspose_SeqDense, 2535 MatMultTransposeAdd_SeqDense, 2536 0, 2537 0, 2538 0, 2539 /* 10*/ 0, 2540 MatLUFactor_SeqDense, 2541 MatCholeskyFactor_SeqDense, 2542 MatSOR_SeqDense, 2543 MatTranspose_SeqDense, 2544 /* 15*/ MatGetInfo_SeqDense, 2545 MatEqual_SeqDense, 2546 MatGetDiagonal_SeqDense, 2547 MatDiagonalScale_SeqDense, 2548 MatNorm_SeqDense, 2549 /* 20*/ MatAssemblyBegin_SeqDense, 2550 MatAssemblyEnd_SeqDense, 2551 MatSetOption_SeqDense, 2552 MatZeroEntries_SeqDense, 2553 /* 24*/ MatZeroRows_SeqDense, 2554 0, 2555 0, 2556 0, 2557 0, 2558 /* 29*/ MatSetUp_SeqDense, 2559 0, 2560 0, 2561 0, 2562 0, 2563 /* 34*/ MatDuplicate_SeqDense, 2564 0, 2565 0, 2566 0, 2567 0, 2568 /* 39*/ MatAXPY_SeqDense, 2569 MatCreateSubMatrices_SeqDense, 2570 0, 2571 MatGetValues_SeqDense, 2572 MatCopy_SeqDense, 2573 /* 44*/ MatGetRowMax_SeqDense, 2574 MatScale_SeqDense, 2575 MatShift_Basic, 2576 0, 2577 MatZeroRowsColumns_SeqDense, 2578 /* 49*/ MatSetRandom_SeqDense, 2579 0, 2580 0, 2581 0, 2582 0, 2583 /* 54*/ 0, 2584 0, 2585 0, 2586 0, 2587 0, 2588 /* 59*/ 0, 2589 MatDestroy_SeqDense, 2590 MatView_SeqDense, 2591 0, 2592 0, 2593 /* 64*/ 0, 2594 0, 2595 0, 2596 0, 2597 0, 2598 /* 69*/ MatGetRowMaxAbs_SeqDense, 2599 0, 2600 0, 2601 0, 2602 0, 2603 /* 74*/ 0, 2604 0, 2605 0, 2606 0, 2607 0, 2608 /* 79*/ 0, 2609 0, 2610 0, 2611 0, 2612 /* 83*/ MatLoad_SeqDense, 2613 MatIsSymmetric_SeqDense, 2614 MatIsHermitian_SeqDense, 2615 0, 2616 0, 2617 0, 2618 /* 89*/ 0, 2619 0, 2620 MatMatMultNumeric_SeqDense_SeqDense, 2621 0, 2622 0, 2623 /* 94*/ 0, 2624 0, 2625 0, 2626 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2627 0, 2628 /* 99*/ MatProductSetFromOptions_SeqDense, 2629 0, 2630 0, 2631 MatConjugate_SeqDense, 2632 0, 2633 /*104*/ 0, 2634 MatRealPart_SeqDense, 2635 MatImaginaryPart_SeqDense, 2636 0, 2637 0, 2638 /*109*/ 0, 2639 0, 2640 MatGetRowMin_SeqDense, 2641 MatGetColumnVector_SeqDense, 2642 MatMissingDiagonal_SeqDense, 2643 /*114*/ 0, 2644 0, 2645 0, 2646 0, 2647 0, 2648 /*119*/ 0, 2649 0, 2650 0, 2651 0, 2652 0, 2653 /*124*/ 0, 2654 MatGetColumnNorms_SeqDense, 2655 0, 2656 0, 2657 0, 2658 /*129*/ 0, 2659 0, 2660 0, 2661 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2662 0, 2663 /*134*/ 0, 2664 0, 2665 0, 2666 0, 2667 0, 2668 /*139*/ 0, 2669 0, 2670 0, 2671 0, 2672 0, 2673 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2674 /*145*/ 0, 2675 0, 2676 0 2677 }; 2678 2679 /*@C 2680 MatCreateSeqDense - Creates a sequential dense matrix that 2681 is stored in column major order (the usual Fortran 77 manner). Many 2682 of the matrix operations use the BLAS and LAPACK routines. 2683 2684 Collective 2685 2686 Input Parameters: 2687 + comm - MPI communicator, set to PETSC_COMM_SELF 2688 . m - number of rows 2689 . n - number of columns 2690 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2691 to control all matrix memory allocation. 2692 2693 Output Parameter: 2694 . A - the matrix 2695 2696 Notes: 2697 The data input variable is intended primarily for Fortran programmers 2698 who wish to allocate their own matrix memory space. Most users should 2699 set data=NULL. 2700 2701 Level: intermediate 2702 2703 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2704 @*/ 2705 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2706 { 2707 PetscErrorCode ierr; 2708 2709 PetscFunctionBegin; 2710 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2711 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2712 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2713 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2714 PetscFunctionReturn(0); 2715 } 2716 2717 /*@C 2718 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2719 2720 Collective 2721 2722 Input Parameters: 2723 + B - the matrix 2724 - data - the array (or NULL) 2725 2726 Notes: 2727 The data input variable is intended primarily for Fortran programmers 2728 who wish to allocate their own matrix memory space. Most users should 2729 need not call this routine. 2730 2731 Level: intermediate 2732 2733 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA() 2734 2735 @*/ 2736 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2737 { 2738 PetscErrorCode ierr; 2739 2740 PetscFunctionBegin; 2741 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2742 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2743 PetscFunctionReturn(0); 2744 } 2745 2746 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2747 { 2748 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2749 PetscErrorCode ierr; 2750 2751 PetscFunctionBegin; 2752 if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2753 B->preallocated = PETSC_TRUE; 2754 2755 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2756 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2757 2758 if (b->lda <= 0) b->lda = B->rmap->n; 2759 2760 ierr = PetscIntMultError(b->lda,B->cmap->n,NULL);CHKERRQ(ierr); 2761 if (!data) { /* petsc-allocated storage */ 2762 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2763 ierr = PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);CHKERRQ(ierr); 2764 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2765 2766 b->user_alloc = PETSC_FALSE; 2767 } else { /* user-allocated storage */ 2768 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2769 b->v = data; 2770 b->user_alloc = PETSC_TRUE; 2771 } 2772 B->assembled = PETSC_TRUE; 2773 PetscFunctionReturn(0); 2774 } 2775 2776 #if defined(PETSC_HAVE_ELEMENTAL) 2777 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2778 { 2779 Mat mat_elemental; 2780 PetscErrorCode ierr; 2781 const PetscScalar *array; 2782 PetscScalar *v_colwise; 2783 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2784 2785 PetscFunctionBegin; 2786 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2787 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2788 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2789 k = 0; 2790 for (j=0; j<N; j++) { 2791 cols[j] = j; 2792 for (i=0; i<M; i++) { 2793 v_colwise[j*M+i] = array[k++]; 2794 } 2795 } 2796 for (i=0; i<M; i++) { 2797 rows[i] = i; 2798 } 2799 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2800 2801 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2802 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2803 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2804 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2805 2806 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2807 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2808 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2809 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2810 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2811 2812 if (reuse == MAT_INPLACE_MATRIX) { 2813 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2814 } else { 2815 *newmat = mat_elemental; 2816 } 2817 PetscFunctionReturn(0); 2818 } 2819 #endif 2820 2821 static PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 2822 { 2823 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2824 2825 PetscFunctionBegin; 2826 if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 2827 b->lda = lda; 2828 PetscFunctionReturn(0); 2829 } 2830 2831 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2832 { 2833 PetscErrorCode ierr; 2834 PetscMPIInt size; 2835 2836 PetscFunctionBegin; 2837 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2838 if (size == 1) { 2839 if (scall == MAT_INITIAL_MATRIX) { 2840 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2841 } else { 2842 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2843 } 2844 } else { 2845 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2846 } 2847 PetscFunctionReturn(0); 2848 } 2849 2850 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2851 { 2852 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2853 PetscErrorCode ierr; 2854 2855 PetscFunctionBegin; 2856 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2857 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2858 if (!a->cvec) { 2859 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2860 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2861 } 2862 a->vecinuse = col + 1; 2863 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2864 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2865 *v = a->cvec; 2866 PetscFunctionReturn(0); 2867 } 2868 2869 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2870 { 2871 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2872 PetscErrorCode ierr; 2873 2874 PetscFunctionBegin; 2875 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2876 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2877 a->vecinuse = 0; 2878 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2879 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2880 *v = NULL; 2881 PetscFunctionReturn(0); 2882 } 2883 2884 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2885 { 2886 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2887 PetscErrorCode ierr; 2888 2889 PetscFunctionBegin; 2890 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2891 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2892 if (!a->cvec) { 2893 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2894 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2895 } 2896 a->vecinuse = col + 1; 2897 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2898 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2899 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr); 2900 *v = a->cvec; 2901 PetscFunctionReturn(0); 2902 } 2903 2904 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2905 { 2906 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2907 PetscErrorCode ierr; 2908 2909 PetscFunctionBegin; 2910 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2911 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2912 a->vecinuse = 0; 2913 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2914 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr); 2915 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2916 *v = NULL; 2917 PetscFunctionReturn(0); 2918 } 2919 2920 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2921 { 2922 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2923 PetscErrorCode ierr; 2924 2925 PetscFunctionBegin; 2926 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2927 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2928 if (!a->cvec) { 2929 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2930 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2931 } 2932 a->vecinuse = col + 1; 2933 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2934 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2935 *v = a->cvec; 2936 PetscFunctionReturn(0); 2937 } 2938 2939 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2940 { 2941 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2942 PetscErrorCode ierr; 2943 2944 PetscFunctionBegin; 2945 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2946 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2947 a->vecinuse = 0; 2948 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2949 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2950 *v = NULL; 2951 PetscFunctionReturn(0); 2952 } 2953 2954 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 2955 { 2956 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2957 PetscErrorCode ierr; 2958 2959 PetscFunctionBegin; 2960 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2961 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2962 if (a->cmat && cend-cbegin != a->cmat->cmap->N) { 2963 ierr = MatDestroy(&a->cmat);CHKERRQ(ierr); 2964 } 2965 if (!a->cmat) { 2966 ierr = MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);CHKERRQ(ierr); 2967 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);CHKERRQ(ierr); 2968 } else { 2969 ierr = MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);CHKERRQ(ierr); 2970 } 2971 ierr = MatDenseSetLDA(a->cmat,a->lda);CHKERRQ(ierr); 2972 a->matinuse = cbegin + 1; 2973 *v = a->cmat; 2974 PetscFunctionReturn(0); 2975 } 2976 2977 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 2978 { 2979 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2980 PetscErrorCode ierr; 2981 2982 PetscFunctionBegin; 2983 if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 2984 if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 2985 if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 2986 a->matinuse = 0; 2987 ierr = MatDenseResetArray(a->cmat);CHKERRQ(ierr); 2988 *v = NULL; 2989 PetscFunctionReturn(0); 2990 } 2991 2992 /*MC 2993 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2994 2995 Options Database Keys: 2996 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2997 2998 Level: beginner 2999 3000 .seealso: MatCreateSeqDense() 3001 3002 M*/ 3003 PetscErrorCode MatCreate_SeqDense(Mat B) 3004 { 3005 Mat_SeqDense *b; 3006 PetscErrorCode ierr; 3007 PetscMPIInt size; 3008 3009 PetscFunctionBegin; 3010 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3011 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3012 3013 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3014 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3015 B->data = (void*)b; 3016 3017 b->roworiented = PETSC_TRUE; 3018 3019 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 3020 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);CHKERRQ(ierr); 3021 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3022 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3023 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 3024 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 3025 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);CHKERRQ(ierr); 3026 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3027 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3028 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3029 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3030 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 3031 #if defined(PETSC_HAVE_ELEMENTAL) 3032 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 3033 #endif 3034 #if defined(PETSC_HAVE_CUDA) 3035 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 3036 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3037 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3038 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3039 #endif 3040 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 3041 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 3042 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3043 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3044 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3045 3046 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 3047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 3048 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr); 3049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr); 3050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr); 3051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr); 3052 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr); 3053 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr); 3054 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);CHKERRQ(ierr); 3055 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);CHKERRQ(ierr); 3056 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 3057 PetscFunctionReturn(0); 3058 } 3059 3060 /*@C 3061 MatDenseGetColumn - gives access to a column of a dense matrix. This is only the local part of the column. You MUST call MatDenseRestoreColumn() to avoid memory bleeding. 3062 3063 Not Collective 3064 3065 Input Parameters: 3066 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3067 - col - column index 3068 3069 Output Parameter: 3070 . vals - pointer to the data 3071 3072 Level: intermediate 3073 3074 .seealso: MatDenseRestoreColumn() 3075 @*/ 3076 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3077 { 3078 PetscErrorCode ierr; 3079 3080 PetscFunctionBegin; 3081 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3082 PetscValidLogicalCollectiveInt(A,col,2); 3083 PetscValidPointer(vals,3); 3084 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 3085 PetscFunctionReturn(0); 3086 } 3087 3088 /*@C 3089 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3090 3091 Not Collective 3092 3093 Input Parameter: 3094 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3095 3096 Output Parameter: 3097 . vals - pointer to the data 3098 3099 Level: intermediate 3100 3101 .seealso: MatDenseGetColumn() 3102 @*/ 3103 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3104 { 3105 PetscErrorCode ierr; 3106 3107 PetscFunctionBegin; 3108 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3109 PetscValidPointer(vals,2); 3110 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 3111 PetscFunctionReturn(0); 3112 } 3113 3114 /*@C 3115 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3116 3117 Collective 3118 3119 Input Parameters: 3120 + mat - the Mat object 3121 - col - the column index 3122 3123 Output Parameter: 3124 . v - the vector 3125 3126 Notes: 3127 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3128 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3129 3130 Level: intermediate 3131 3132 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3133 @*/ 3134 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3135 { 3136 PetscErrorCode ierr; 3137 3138 PetscFunctionBegin; 3139 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3140 PetscValidType(A,1); 3141 PetscValidLogicalCollectiveInt(A,col,2); 3142 PetscValidPointer(v,3); 3143 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3144 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3145 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3146 PetscFunctionReturn(0); 3147 } 3148 3149 /*@C 3150 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3151 3152 Collective 3153 3154 Input Parameters: 3155 + mat - the Mat object 3156 . col - the column index 3157 - v - the Vec object 3158 3159 Level: intermediate 3160 3161 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3162 @*/ 3163 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3164 { 3165 PetscErrorCode ierr; 3166 3167 PetscFunctionBegin; 3168 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3169 PetscValidType(A,1); 3170 PetscValidLogicalCollectiveInt(A,col,2); 3171 PetscValidPointer(v,3); 3172 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3173 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3174 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3175 PetscFunctionReturn(0); 3176 } 3177 3178 /*@C 3179 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3180 3181 Collective 3182 3183 Input Parameters: 3184 + mat - the Mat object 3185 - col - the column index 3186 3187 Output Parameter: 3188 . v - the vector 3189 3190 Notes: 3191 The vector is owned by PETSc and users cannot modify it. 3192 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3193 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3194 3195 Level: intermediate 3196 3197 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3198 @*/ 3199 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3200 { 3201 PetscErrorCode ierr; 3202 3203 PetscFunctionBegin; 3204 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3205 PetscValidType(A,1); 3206 PetscValidLogicalCollectiveInt(A,col,2); 3207 PetscValidPointer(v,3); 3208 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3209 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3210 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3211 PetscFunctionReturn(0); 3212 } 3213 3214 /*@C 3215 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3216 3217 Collective 3218 3219 Input Parameters: 3220 + mat - the Mat object 3221 . col - the column index 3222 - v - the Vec object 3223 3224 Level: intermediate 3225 3226 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite() 3227 @*/ 3228 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3229 { 3230 PetscErrorCode ierr; 3231 3232 PetscFunctionBegin; 3233 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3234 PetscValidType(A,1); 3235 PetscValidLogicalCollectiveInt(A,col,2); 3236 PetscValidPointer(v,3); 3237 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3238 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3239 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3240 PetscFunctionReturn(0); 3241 } 3242 3243 /*@C 3244 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3245 3246 Collective 3247 3248 Input Parameters: 3249 + mat - the Mat object 3250 - col - the column index 3251 3252 Output Parameter: 3253 . v - the vector 3254 3255 Notes: 3256 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3257 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3258 3259 Level: intermediate 3260 3261 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3262 @*/ 3263 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3264 { 3265 PetscErrorCode ierr; 3266 3267 PetscFunctionBegin; 3268 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3269 PetscValidType(A,1); 3270 PetscValidLogicalCollectiveInt(A,col,2); 3271 PetscValidPointer(v,3); 3272 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3273 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3274 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3275 PetscFunctionReturn(0); 3276 } 3277 3278 /*@C 3279 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3280 3281 Collective 3282 3283 Input Parameters: 3284 + mat - the Mat object 3285 . col - the column index 3286 - v - the Vec object 3287 3288 Level: intermediate 3289 3290 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead() 3291 @*/ 3292 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3293 { 3294 PetscErrorCode ierr; 3295 3296 PetscFunctionBegin; 3297 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3298 PetscValidType(A,1); 3299 PetscValidLogicalCollectiveInt(A,col,2); 3300 PetscValidPointer(v,3); 3301 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3302 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3303 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3304 PetscFunctionReturn(0); 3305 } 3306 3307 /*@C 3308 MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat. 3309 3310 Collective 3311 3312 Input Parameters: 3313 + mat - the Mat object 3314 . cbegin - the first index in the block 3315 - cend - the last index in the block 3316 3317 Output Parameter: 3318 . v - the matrix 3319 3320 Notes: 3321 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3322 3323 Level: intermediate 3324 3325 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix() 3326 @*/ 3327 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 3328 { 3329 PetscErrorCode ierr; 3330 3331 PetscFunctionBegin; 3332 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3333 PetscValidType(A,1); 3334 PetscValidLogicalCollectiveInt(A,cbegin,2); 3335 PetscValidLogicalCollectiveInt(A,cend,3); 3336 PetscValidPointer(v,4); 3337 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3338 if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N); 3339 if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N); 3340 ierr = PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));CHKERRQ(ierr); 3341 PetscFunctionReturn(0); 3342 } 3343 3344 /*@C 3345 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3346 3347 Collective 3348 3349 Input Parameters: 3350 + mat - the Mat object 3351 - v - the Mat object 3352 3353 Level: intermediate 3354 3355 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix() 3356 @*/ 3357 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3358 { 3359 PetscErrorCode ierr; 3360 3361 PetscFunctionBegin; 3362 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3363 PetscValidType(A,1); 3364 PetscValidPointer(v,2); 3365 ierr = PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));CHKERRQ(ierr); 3366 PetscFunctionReturn(0); 3367 } 3368