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