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