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