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,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1499 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1500 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1501 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1502 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1503 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1504 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1505 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1506 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1507 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1508 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1509 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1510 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_nest_seqdense_C",NULL);CHKERRQ(ierr); 1511 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_nest_seqdense_C",NULL);CHKERRQ(ierr); 1512 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_nest_seqdense_C",NULL);CHKERRQ(ierr); 1513 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1514 1515 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1516 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1517 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1518 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1519 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1520 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1521 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1522 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1523 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1524 1525 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1526 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1527 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1528 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1529 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1530 PetscFunctionReturn(0); 1531 } 1532 1533 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1534 { 1535 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1536 PetscErrorCode ierr; 1537 PetscInt k,j,m,n,M; 1538 PetscScalar *v,tmp; 1539 1540 PetscFunctionBegin; 1541 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1542 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1543 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1544 for (j=0; j<m; j++) { 1545 for (k=0; k<j; k++) { 1546 tmp = v[j + k*M]; 1547 v[j + k*M] = v[k + j*M]; 1548 v[k + j*M] = tmp; 1549 } 1550 } 1551 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1552 } else { /* out-of-place transpose */ 1553 Mat tmat; 1554 Mat_SeqDense *tmatd; 1555 PetscScalar *v2; 1556 PetscInt M2; 1557 1558 if (reuse != MAT_REUSE_MATRIX) { 1559 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1560 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1561 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1562 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1563 } else tmat = *matout; 1564 1565 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1566 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1567 tmatd = (Mat_SeqDense*)tmat->data; 1568 M2 = tmatd->lda; 1569 for (j=0; j<n; j++) { 1570 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1571 } 1572 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1573 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1574 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1575 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1576 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1577 else { 1578 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1579 } 1580 } 1581 PetscFunctionReturn(0); 1582 } 1583 1584 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1585 { 1586 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1587 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1588 PetscInt i; 1589 const PetscScalar *v1,*v2; 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1594 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1595 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1596 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1597 for (i=0; i<A1->cmap->n; i++) { 1598 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1599 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1600 v1 += mat1->lda; 1601 v2 += mat2->lda; 1602 } 1603 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1604 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1605 *flg = PETSC_TRUE; 1606 PetscFunctionReturn(0); 1607 } 1608 1609 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1610 { 1611 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1612 PetscInt i,n,len; 1613 PetscScalar *x; 1614 const PetscScalar *vv; 1615 PetscErrorCode ierr; 1616 1617 PetscFunctionBegin; 1618 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1619 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1620 len = PetscMin(A->rmap->n,A->cmap->n); 1621 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1622 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1623 for (i=0; i<len; i++) { 1624 x[i] = vv[i*mat->lda + i]; 1625 } 1626 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1627 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1632 { 1633 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1634 const PetscScalar *l,*r; 1635 PetscScalar x,*v,*vv; 1636 PetscErrorCode ierr; 1637 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1638 1639 PetscFunctionBegin; 1640 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1641 if (ll) { 1642 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1643 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1644 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1645 for (i=0; i<m; i++) { 1646 x = l[i]; 1647 v = vv + i; 1648 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1649 } 1650 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1651 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1652 } 1653 if (rr) { 1654 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1655 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1656 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1657 for (i=0; i<n; i++) { 1658 x = r[i]; 1659 v = vv + i*mat->lda; 1660 for (j=0; j<m; j++) (*v++) *= x; 1661 } 1662 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1663 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1664 } 1665 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 1669 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1670 { 1671 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1672 PetscScalar *v,*vv; 1673 PetscReal sum = 0.0; 1674 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1675 PetscErrorCode ierr; 1676 1677 PetscFunctionBegin; 1678 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1679 v = vv; 1680 if (type == NORM_FROBENIUS) { 1681 if (lda>m) { 1682 for (j=0; j<A->cmap->n; j++) { 1683 v = vv+j*lda; 1684 for (i=0; i<m; i++) { 1685 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1686 } 1687 } 1688 } else { 1689 #if defined(PETSC_USE_REAL___FP16) 1690 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1691 *nrm = BLASnrm2_(&cnt,v,&one); 1692 } 1693 #else 1694 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1695 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1696 } 1697 } 1698 *nrm = PetscSqrtReal(sum); 1699 #endif 1700 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1701 } else if (type == NORM_1) { 1702 *nrm = 0.0; 1703 for (j=0; j<A->cmap->n; j++) { 1704 v = vv + j*mat->lda; 1705 sum = 0.0; 1706 for (i=0; i<A->rmap->n; i++) { 1707 sum += PetscAbsScalar(*v); v++; 1708 } 1709 if (sum > *nrm) *nrm = sum; 1710 } 1711 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1712 } else if (type == NORM_INFINITY) { 1713 *nrm = 0.0; 1714 for (j=0; j<A->rmap->n; j++) { 1715 v = vv + j; 1716 sum = 0.0; 1717 for (i=0; i<A->cmap->n; i++) { 1718 sum += PetscAbsScalar(*v); v += mat->lda; 1719 } 1720 if (sum > *nrm) *nrm = sum; 1721 } 1722 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1723 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1724 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1725 PetscFunctionReturn(0); 1726 } 1727 1728 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1729 { 1730 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1731 PetscErrorCode ierr; 1732 1733 PetscFunctionBegin; 1734 switch (op) { 1735 case MAT_ROW_ORIENTED: 1736 aij->roworiented = flg; 1737 break; 1738 case MAT_NEW_NONZERO_LOCATIONS: 1739 case MAT_NEW_NONZERO_LOCATION_ERR: 1740 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1741 case MAT_NEW_DIAGONALS: 1742 case MAT_KEEP_NONZERO_PATTERN: 1743 case MAT_IGNORE_OFF_PROC_ENTRIES: 1744 case MAT_USE_HASH_TABLE: 1745 case MAT_IGNORE_ZERO_ENTRIES: 1746 case MAT_IGNORE_LOWER_TRIANGULAR: 1747 case MAT_SORTED_FULL: 1748 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1749 break; 1750 case MAT_SPD: 1751 case MAT_SYMMETRIC: 1752 case MAT_STRUCTURALLY_SYMMETRIC: 1753 case MAT_HERMITIAN: 1754 case MAT_SYMMETRY_ETERNAL: 1755 /* These options are handled directly by MatSetOption() */ 1756 break; 1757 default: 1758 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1759 } 1760 PetscFunctionReturn(0); 1761 } 1762 1763 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1764 { 1765 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1766 PetscErrorCode ierr; 1767 PetscInt lda=l->lda,m=A->rmap->n,j; 1768 PetscScalar *v; 1769 1770 PetscFunctionBegin; 1771 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1772 if (lda>m) { 1773 for (j=0; j<A->cmap->n; j++) { 1774 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1775 } 1776 } else { 1777 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1778 } 1779 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1780 PetscFunctionReturn(0); 1781 } 1782 1783 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1784 { 1785 PetscErrorCode ierr; 1786 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1787 PetscInt m = l->lda, n = A->cmap->n, i,j; 1788 PetscScalar *slot,*bb,*v; 1789 const PetscScalar *xx; 1790 1791 PetscFunctionBegin; 1792 #if defined(PETSC_USE_DEBUG) 1793 for (i=0; i<N; i++) { 1794 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1795 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); 1796 } 1797 #endif 1798 if (!N) PetscFunctionReturn(0); 1799 1800 /* fix right hand side if needed */ 1801 if (x && b) { 1802 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1803 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1804 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1805 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1806 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1807 } 1808 1809 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1810 for (i=0; i<N; i++) { 1811 slot = v + rows[i]; 1812 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1813 } 1814 if (diag != 0.0) { 1815 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1816 for (i=0; i<N; i++) { 1817 slot = v + (m+1)*rows[i]; 1818 *slot = diag; 1819 } 1820 } 1821 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1822 PetscFunctionReturn(0); 1823 } 1824 1825 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1826 { 1827 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1828 1829 PetscFunctionBegin; 1830 *lda = mat->lda; 1831 PetscFunctionReturn(0); 1832 } 1833 1834 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1835 { 1836 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1837 1838 PetscFunctionBegin; 1839 *array = mat->v; 1840 PetscFunctionReturn(0); 1841 } 1842 1843 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1844 { 1845 PetscFunctionBegin; 1846 PetscFunctionReturn(0); 1847 } 1848 1849 /*@C 1850 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1851 1852 Logically Collective on Mat 1853 1854 Input Parameter: 1855 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1856 1857 Output Parameter: 1858 . lda - the leading dimension 1859 1860 Level: intermediate 1861 1862 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1863 @*/ 1864 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1865 { 1866 PetscErrorCode ierr; 1867 1868 PetscFunctionBegin; 1869 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /*@C 1874 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1875 1876 Logically Collective on Mat 1877 1878 Input Parameter: 1879 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1880 1881 Output Parameter: 1882 . array - pointer to the data 1883 1884 Level: intermediate 1885 1886 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1887 @*/ 1888 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1889 { 1890 PetscErrorCode ierr; 1891 1892 PetscFunctionBegin; 1893 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /*@C 1898 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1899 1900 Logically Collective on Mat 1901 1902 Input Parameters: 1903 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1904 - array - pointer to the data 1905 1906 Level: intermediate 1907 1908 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1909 @*/ 1910 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1911 { 1912 PetscErrorCode ierr; 1913 1914 PetscFunctionBegin; 1915 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1916 if (array) *array = NULL; 1917 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 1921 /*@C 1922 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1923 1924 Not Collective 1925 1926 Input Parameter: 1927 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1928 1929 Output Parameter: 1930 . array - pointer to the data 1931 1932 Level: intermediate 1933 1934 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1935 @*/ 1936 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1937 { 1938 PetscErrorCode ierr; 1939 1940 PetscFunctionBegin; 1941 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 /*@C 1946 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1947 1948 Not Collective 1949 1950 Input Parameters: 1951 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1952 - array - pointer to the data 1953 1954 Level: intermediate 1955 1956 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1957 @*/ 1958 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1959 { 1960 PetscErrorCode ierr; 1961 1962 PetscFunctionBegin; 1963 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1964 if (array) *array = NULL; 1965 PetscFunctionReturn(0); 1966 } 1967 1968 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1969 { 1970 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1971 PetscErrorCode ierr; 1972 PetscInt i,j,nrows,ncols,blda; 1973 const PetscInt *irow,*icol; 1974 PetscScalar *av,*bv,*v = mat->v; 1975 Mat newmat; 1976 1977 PetscFunctionBegin; 1978 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1979 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1980 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1981 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1982 1983 /* Check submatrixcall */ 1984 if (scall == MAT_REUSE_MATRIX) { 1985 PetscInt n_cols,n_rows; 1986 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1987 if (n_rows != nrows || n_cols != ncols) { 1988 /* resize the result matrix to match number of requested rows/columns */ 1989 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1990 } 1991 newmat = *B; 1992 } else { 1993 /* Create and fill new matrix */ 1994 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1995 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1996 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1997 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1998 } 1999 2000 /* Now extract the data pointers and do the copy,column at a time */ 2001 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 2002 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 2003 for (i=0; i<ncols; i++) { 2004 av = v + mat->lda*icol[i]; 2005 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2006 bv += blda; 2007 } 2008 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 2009 2010 /* Assemble the matrices so that the correct flags are set */ 2011 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2012 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2013 2014 /* Free work space */ 2015 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2016 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2017 *B = newmat; 2018 PetscFunctionReturn(0); 2019 } 2020 2021 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2022 { 2023 PetscErrorCode ierr; 2024 PetscInt i; 2025 2026 PetscFunctionBegin; 2027 if (scall == MAT_INITIAL_MATRIX) { 2028 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2029 } 2030 2031 for (i=0; i<n; i++) { 2032 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2033 } 2034 PetscFunctionReturn(0); 2035 } 2036 2037 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2038 { 2039 PetscFunctionBegin; 2040 PetscFunctionReturn(0); 2041 } 2042 2043 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2044 { 2045 PetscFunctionBegin; 2046 PetscFunctionReturn(0); 2047 } 2048 2049 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2050 { 2051 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2052 PetscErrorCode ierr; 2053 const PetscScalar *va; 2054 PetscScalar *vb; 2055 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2056 2057 PetscFunctionBegin; 2058 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2059 if (A->ops->copy != B->ops->copy) { 2060 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2061 PetscFunctionReturn(0); 2062 } 2063 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2064 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2065 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2066 if (lda1>m || lda2>m) { 2067 for (j=0; j<n; j++) { 2068 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2069 } 2070 } else { 2071 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2072 } 2073 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2074 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2075 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2076 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2077 PetscFunctionReturn(0); 2078 } 2079 2080 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2081 { 2082 PetscErrorCode ierr; 2083 2084 PetscFunctionBegin; 2085 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2086 PetscFunctionReturn(0); 2087 } 2088 2089 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2090 { 2091 PetscInt i,nz = A->rmap->n*A->cmap->n; 2092 PetscScalar *aa; 2093 PetscErrorCode ierr; 2094 2095 PetscFunctionBegin; 2096 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2097 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2098 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2099 PetscFunctionReturn(0); 2100 } 2101 2102 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2103 { 2104 PetscInt i,nz = A->rmap->n*A->cmap->n; 2105 PetscScalar *aa; 2106 PetscErrorCode ierr; 2107 2108 PetscFunctionBegin; 2109 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2110 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2111 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2112 PetscFunctionReturn(0); 2113 } 2114 2115 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2116 { 2117 PetscInt i,nz = A->rmap->n*A->cmap->n; 2118 PetscScalar *aa; 2119 PetscErrorCode ierr; 2120 2121 PetscFunctionBegin; 2122 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2123 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2124 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2125 PetscFunctionReturn(0); 2126 } 2127 2128 /* ----------------------------------------------------------------*/ 2129 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2130 { 2131 PetscErrorCode ierr; 2132 2133 PetscFunctionBegin; 2134 if (scall == MAT_INITIAL_MATRIX) { 2135 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2136 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2137 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2138 } 2139 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2140 if ((*C)->ops->matmultnumeric) { 2141 ierr = (*(*C)->ops->matmultnumeric)(A,B,*C);CHKERRQ(ierr); 2142 } else { 2143 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2144 } 2145 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2146 PetscFunctionReturn(0); 2147 } 2148 2149 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2150 { 2151 PetscErrorCode ierr; 2152 PetscInt m=A->rmap->n,n=B->cmap->n; 2153 Mat Cmat; 2154 PetscBool flg; 2155 2156 PetscFunctionBegin; 2157 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2158 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2159 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2160 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2161 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2162 *C = Cmat; 2163 PetscFunctionReturn(0); 2164 } 2165 2166 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2167 { 2168 Mat_SeqDense *a,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2169 PetscBLASInt m,n,k; 2170 const PetscScalar *av,*bv; 2171 PetscScalar *cv; 2172 PetscScalar _DOne=1.0,_DZero=0.0; 2173 PetscErrorCode ierr; 2174 PetscBool flg; 2175 2176 PetscFunctionBegin; 2177 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2178 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2179 if (flg) { 2180 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2181 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2182 PetscFunctionReturn(0); 2183 } 2184 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);CHKERRQ(ierr); 2185 if (flg) { 2186 C->ops->matmultnumeric = MatMatMultNumeric_SeqBAIJ_SeqDense; 2187 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2188 PetscFunctionReturn(0); 2189 } 2190 ierr = PetscObjectTypeCompare((PetscObject)A,MATNEST,&flg);CHKERRQ(ierr); 2191 if (flg) { 2192 C->ops->matmultnumeric = MatMatMultNumeric_Nest_Dense; 2193 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2194 PetscFunctionReturn(0); 2195 } 2196 a = (Mat_SeqDense*)A->data; 2197 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2198 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2199 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2200 if (!m || !n || !k) PetscFunctionReturn(0); 2201 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2202 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2203 ierr = MatDenseGetArray(C,&cv);CHKERRQ(ierr); 2204 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2205 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2206 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2207 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2208 ierr = MatDenseRestoreArray(C,&cv);CHKERRQ(ierr); 2209 PetscFunctionReturn(0); 2210 } 2211 2212 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2213 { 2214 PetscErrorCode ierr; 2215 2216 PetscFunctionBegin; 2217 if (scall == MAT_INITIAL_MATRIX) { 2218 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2219 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2220 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2221 } 2222 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2223 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2224 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2225 PetscFunctionReturn(0); 2226 } 2227 2228 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2229 { 2230 PetscErrorCode ierr; 2231 PetscInt m=A->rmap->n,n=B->rmap->n; 2232 Mat Cmat; 2233 PetscBool flg; 2234 2235 PetscFunctionBegin; 2236 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2237 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2238 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2239 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2240 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2241 *C = Cmat; 2242 PetscFunctionReturn(0); 2243 } 2244 2245 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2246 { 2247 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2248 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2249 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2250 PetscBLASInt m,n,k; 2251 PetscScalar _DOne=1.0,_DZero=0.0; 2252 PetscErrorCode ierr; 2253 2254 PetscFunctionBegin; 2255 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2256 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2257 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2258 if (!m || !n || !k) PetscFunctionReturn(0); 2259 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2260 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2261 PetscFunctionReturn(0); 2262 } 2263 2264 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2265 { 2266 PetscErrorCode ierr; 2267 2268 PetscFunctionBegin; 2269 if (scall == MAT_INITIAL_MATRIX) { 2270 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2271 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2272 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2273 } 2274 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2275 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2276 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2277 PetscFunctionReturn(0); 2278 } 2279 2280 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2281 { 2282 PetscErrorCode ierr; 2283 PetscInt m=A->cmap->n,n=B->cmap->n; 2284 Mat Cmat; 2285 PetscBool flg; 2286 2287 PetscFunctionBegin; 2288 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2289 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2290 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2291 ierr = MatSetType(Cmat,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2292 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2293 *C = Cmat; 2294 PetscFunctionReturn(0); 2295 } 2296 2297 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2298 { 2299 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2300 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2301 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2302 PetscBLASInt m,n,k; 2303 PetscScalar _DOne=1.0,_DZero=0.0; 2304 PetscErrorCode ierr; 2305 2306 PetscFunctionBegin; 2307 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2308 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2309 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2310 if (!m || !n || !k) PetscFunctionReturn(0); 2311 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2312 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2313 PetscFunctionReturn(0); 2314 } 2315 2316 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2317 { 2318 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2319 PetscErrorCode ierr; 2320 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2321 PetscScalar *x; 2322 const PetscScalar *aa; 2323 2324 PetscFunctionBegin; 2325 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2326 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2327 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2328 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2329 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2330 for (i=0; i<m; i++) { 2331 x[i] = aa[i]; if (idx) idx[i] = 0; 2332 for (j=1; j<n; j++) { 2333 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2334 } 2335 } 2336 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2337 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2338 PetscFunctionReturn(0); 2339 } 2340 2341 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2342 { 2343 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2344 PetscErrorCode ierr; 2345 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2346 PetscScalar *x; 2347 PetscReal atmp; 2348 const PetscScalar *aa; 2349 2350 PetscFunctionBegin; 2351 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2352 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2353 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2354 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2355 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2356 for (i=0; i<m; i++) { 2357 x[i] = PetscAbsScalar(aa[i]); 2358 for (j=1; j<n; j++) { 2359 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2360 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2361 } 2362 } 2363 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2364 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2365 PetscFunctionReturn(0); 2366 } 2367 2368 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2369 { 2370 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2371 PetscErrorCode ierr; 2372 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2373 PetscScalar *x; 2374 const PetscScalar *aa; 2375 2376 PetscFunctionBegin; 2377 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2378 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2379 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2380 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2381 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2382 for (i=0; i<m; i++) { 2383 x[i] = aa[i]; if (idx) idx[i] = 0; 2384 for (j=1; j<n; j++) { 2385 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2386 } 2387 } 2388 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2389 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2390 PetscFunctionReturn(0); 2391 } 2392 2393 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2394 { 2395 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2396 PetscErrorCode ierr; 2397 PetscScalar *x; 2398 const PetscScalar *aa; 2399 2400 PetscFunctionBegin; 2401 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2402 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2403 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2404 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2405 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2406 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2407 PetscFunctionReturn(0); 2408 } 2409 2410 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2411 { 2412 PetscErrorCode ierr; 2413 PetscInt i,j,m,n; 2414 const PetscScalar *a; 2415 2416 PetscFunctionBegin; 2417 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2418 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2419 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2420 if (type == NORM_2) { 2421 for (i=0; i<n; i++) { 2422 for (j=0; j<m; j++) { 2423 norms[i] += PetscAbsScalar(a[j]*a[j]); 2424 } 2425 a += m; 2426 } 2427 } else if (type == NORM_1) { 2428 for (i=0; i<n; i++) { 2429 for (j=0; j<m; j++) { 2430 norms[i] += PetscAbsScalar(a[j]); 2431 } 2432 a += m; 2433 } 2434 } else if (type == NORM_INFINITY) { 2435 for (i=0; i<n; i++) { 2436 for (j=0; j<m; j++) { 2437 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2438 } 2439 a += m; 2440 } 2441 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2442 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2443 if (type == NORM_2) { 2444 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2445 } 2446 PetscFunctionReturn(0); 2447 } 2448 2449 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2450 { 2451 PetscErrorCode ierr; 2452 PetscScalar *a; 2453 PetscInt m,n,i; 2454 2455 PetscFunctionBegin; 2456 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2457 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2458 for (i=0; i<m*n; i++) { 2459 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2460 } 2461 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2462 PetscFunctionReturn(0); 2463 } 2464 2465 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2466 { 2467 PetscFunctionBegin; 2468 *missing = PETSC_FALSE; 2469 PetscFunctionReturn(0); 2470 } 2471 2472 /* vals is not const */ 2473 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2474 { 2475 PetscErrorCode ierr; 2476 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2477 PetscScalar *v; 2478 2479 PetscFunctionBegin; 2480 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2481 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2482 *vals = v+col*a->lda; 2483 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2484 PetscFunctionReturn(0); 2485 } 2486 2487 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2488 { 2489 PetscFunctionBegin; 2490 *vals = 0; /* user cannot accidently use the array later */ 2491 PetscFunctionReturn(0); 2492 } 2493 2494 /* -------------------------------------------------------------------*/ 2495 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2496 MatGetRow_SeqDense, 2497 MatRestoreRow_SeqDense, 2498 MatMult_SeqDense, 2499 /* 4*/ MatMultAdd_SeqDense, 2500 MatMultTranspose_SeqDense, 2501 MatMultTransposeAdd_SeqDense, 2502 0, 2503 0, 2504 0, 2505 /* 10*/ 0, 2506 MatLUFactor_SeqDense, 2507 MatCholeskyFactor_SeqDense, 2508 MatSOR_SeqDense, 2509 MatTranspose_SeqDense, 2510 /* 15*/ MatGetInfo_SeqDense, 2511 MatEqual_SeqDense, 2512 MatGetDiagonal_SeqDense, 2513 MatDiagonalScale_SeqDense, 2514 MatNorm_SeqDense, 2515 /* 20*/ MatAssemblyBegin_SeqDense, 2516 MatAssemblyEnd_SeqDense, 2517 MatSetOption_SeqDense, 2518 MatZeroEntries_SeqDense, 2519 /* 24*/ MatZeroRows_SeqDense, 2520 0, 2521 0, 2522 0, 2523 0, 2524 /* 29*/ MatSetUp_SeqDense, 2525 0, 2526 0, 2527 0, 2528 0, 2529 /* 34*/ MatDuplicate_SeqDense, 2530 0, 2531 0, 2532 0, 2533 0, 2534 /* 39*/ MatAXPY_SeqDense, 2535 MatCreateSubMatrices_SeqDense, 2536 0, 2537 MatGetValues_SeqDense, 2538 MatCopy_SeqDense, 2539 /* 44*/ MatGetRowMax_SeqDense, 2540 MatScale_SeqDense, 2541 MatShift_Basic, 2542 0, 2543 MatZeroRowsColumns_SeqDense, 2544 /* 49*/ MatSetRandom_SeqDense, 2545 0, 2546 0, 2547 0, 2548 0, 2549 /* 54*/ 0, 2550 0, 2551 0, 2552 0, 2553 0, 2554 /* 59*/ 0, 2555 MatDestroy_SeqDense, 2556 MatView_SeqDense, 2557 0, 2558 0, 2559 /* 64*/ 0, 2560 0, 2561 0, 2562 0, 2563 0, 2564 /* 69*/ MatGetRowMaxAbs_SeqDense, 2565 0, 2566 0, 2567 0, 2568 0, 2569 /* 74*/ 0, 2570 0, 2571 0, 2572 0, 2573 0, 2574 /* 79*/ 0, 2575 0, 2576 0, 2577 0, 2578 /* 83*/ MatLoad_SeqDense, 2579 0, 2580 MatIsHermitian_SeqDense, 2581 0, 2582 0, 2583 0, 2584 /* 89*/ MatMatMult_SeqDense_SeqDense, 2585 MatMatMultSymbolic_SeqDense_SeqDense, 2586 MatMatMultNumeric_SeqDense_SeqDense, 2587 MatPtAP_SeqDense_SeqDense, 2588 MatPtAPSymbolic_SeqDense_SeqDense, 2589 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2590 MatMatTransposeMult_SeqDense_SeqDense, 2591 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2592 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2593 0, 2594 /* 99*/ 0, 2595 0, 2596 0, 2597 MatConjugate_SeqDense, 2598 0, 2599 /*104*/ 0, 2600 MatRealPart_SeqDense, 2601 MatImaginaryPart_SeqDense, 2602 0, 2603 0, 2604 /*109*/ 0, 2605 0, 2606 MatGetRowMin_SeqDense, 2607 MatGetColumnVector_SeqDense, 2608 MatMissingDiagonal_SeqDense, 2609 /*114*/ 0, 2610 0, 2611 0, 2612 0, 2613 0, 2614 /*119*/ 0, 2615 0, 2616 0, 2617 0, 2618 0, 2619 /*124*/ 0, 2620 MatGetColumnNorms_SeqDense, 2621 0, 2622 0, 2623 0, 2624 /*129*/ 0, 2625 MatTransposeMatMult_SeqDense_SeqDense, 2626 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2627 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2628 0, 2629 /*134*/ 0, 2630 0, 2631 0, 2632 0, 2633 0, 2634 /*139*/ 0, 2635 0, 2636 0, 2637 0, 2638 0, 2639 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2640 }; 2641 2642 /*@C 2643 MatCreateSeqDense - Creates a sequential dense matrix that 2644 is stored in column major order (the usual Fortran 77 manner). Many 2645 of the matrix operations use the BLAS and LAPACK routines. 2646 2647 Collective 2648 2649 Input Parameters: 2650 + comm - MPI communicator, set to PETSC_COMM_SELF 2651 . m - number of rows 2652 . n - number of columns 2653 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2654 to control all matrix memory allocation. 2655 2656 Output Parameter: 2657 . A - the matrix 2658 2659 Notes: 2660 The data input variable is intended primarily for Fortran programmers 2661 who wish to allocate their own matrix memory space. Most users should 2662 set data=NULL. 2663 2664 Level: intermediate 2665 2666 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2667 @*/ 2668 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2669 { 2670 PetscErrorCode ierr; 2671 2672 PetscFunctionBegin; 2673 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2674 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2675 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2676 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2677 PetscFunctionReturn(0); 2678 } 2679 2680 /*@C 2681 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2682 2683 Collective 2684 2685 Input Parameters: 2686 + B - the matrix 2687 - data - the array (or NULL) 2688 2689 Notes: 2690 The data input variable is intended primarily for Fortran programmers 2691 who wish to allocate their own matrix memory space. Most users should 2692 need not call this routine. 2693 2694 Level: intermediate 2695 2696 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2697 2698 @*/ 2699 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2700 { 2701 PetscErrorCode ierr; 2702 2703 PetscFunctionBegin; 2704 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2705 PetscFunctionReturn(0); 2706 } 2707 2708 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2709 { 2710 Mat_SeqDense *b; 2711 PetscErrorCode ierr; 2712 2713 PetscFunctionBegin; 2714 B->preallocated = PETSC_TRUE; 2715 2716 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2717 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2718 2719 b = (Mat_SeqDense*)B->data; 2720 b->Mmax = B->rmap->n; 2721 b->Nmax = B->cmap->n; 2722 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2723 2724 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2725 if (!data) { /* petsc-allocated storage */ 2726 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2727 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2728 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2729 2730 b->user_alloc = PETSC_FALSE; 2731 } else { /* user-allocated storage */ 2732 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2733 b->v = data; 2734 b->user_alloc = PETSC_TRUE; 2735 } 2736 B->assembled = PETSC_TRUE; 2737 PetscFunctionReturn(0); 2738 } 2739 2740 #if defined(PETSC_HAVE_ELEMENTAL) 2741 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2742 { 2743 Mat mat_elemental; 2744 PetscErrorCode ierr; 2745 const PetscScalar *array; 2746 PetscScalar *v_colwise; 2747 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2748 2749 PetscFunctionBegin; 2750 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2751 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2752 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2753 k = 0; 2754 for (j=0; j<N; j++) { 2755 cols[j] = j; 2756 for (i=0; i<M; i++) { 2757 v_colwise[j*M+i] = array[k++]; 2758 } 2759 } 2760 for (i=0; i<M; i++) { 2761 rows[i] = i; 2762 } 2763 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2764 2765 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2766 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2767 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2768 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2769 2770 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2771 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2772 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2773 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2774 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2775 2776 if (reuse == MAT_INPLACE_MATRIX) { 2777 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2778 } else { 2779 *newmat = mat_elemental; 2780 } 2781 PetscFunctionReturn(0); 2782 } 2783 #endif 2784 2785 /*@C 2786 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2787 2788 Input parameter: 2789 + A - the matrix 2790 - lda - the leading dimension 2791 2792 Notes: 2793 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2794 it asserts that the preallocation has a leading dimension (the LDA parameter 2795 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2796 2797 Level: intermediate 2798 2799 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2800 2801 @*/ 2802 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2803 { 2804 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2805 2806 PetscFunctionBegin; 2807 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); 2808 b->lda = lda; 2809 b->changelda = PETSC_FALSE; 2810 b->Mmax = PetscMax(b->Mmax,lda); 2811 PetscFunctionReturn(0); 2812 } 2813 2814 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2815 { 2816 PetscErrorCode ierr; 2817 PetscMPIInt size; 2818 2819 PetscFunctionBegin; 2820 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2821 if (size == 1) { 2822 if (scall == MAT_INITIAL_MATRIX) { 2823 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2824 } else { 2825 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2826 } 2827 } else { 2828 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2829 } 2830 PetscFunctionReturn(0); 2831 } 2832 2833 /*MC 2834 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2835 2836 Options Database Keys: 2837 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2838 2839 Level: beginner 2840 2841 .seealso: MatCreateSeqDense() 2842 2843 M*/ 2844 PetscErrorCode MatCreate_SeqDense(Mat B) 2845 { 2846 Mat_SeqDense *b; 2847 PetscErrorCode ierr; 2848 PetscMPIInt size; 2849 2850 PetscFunctionBegin; 2851 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2852 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2853 2854 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2855 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2856 B->data = (void*)b; 2857 2858 b->roworiented = PETSC_TRUE; 2859 2860 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2861 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2862 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2863 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2864 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2865 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2866 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2867 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2868 #if defined(PETSC_HAVE_ELEMENTAL) 2869 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2870 #endif 2871 #if defined(PETSC_HAVE_CUDA) 2872 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 2873 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijcusparse_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2874 #endif 2875 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2876 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2877 #if defined(PETSC_HAVE_VIENNACL) 2878 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijviennacl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2879 #endif 2880 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2881 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2882 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqbaij_seqdense_C",MatMatMult_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2883 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2884 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2885 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2886 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2887 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2888 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2889 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2890 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2891 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2892 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2893 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2894 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2895 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2896 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2897 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_nest_seqdense_C",MatMatMult_Nest_Dense);CHKERRQ(ierr); 2898 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);CHKERRQ(ierr); 2899 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);CHKERRQ(ierr); 2900 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2901 2902 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2903 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2904 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2905 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2906 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2907 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2908 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2909 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2910 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2911 2912 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2914 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2915 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2916 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2917 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2918 PetscFunctionReturn(0); 2919 } 2920 2921 /*@C 2922 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. 2923 2924 Not Collective 2925 2926 Input Parameter: 2927 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2928 - col - column index 2929 2930 Output Parameter: 2931 . vals - pointer to the data 2932 2933 Level: intermediate 2934 2935 .seealso: MatDenseRestoreColumn() 2936 @*/ 2937 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2938 { 2939 PetscErrorCode ierr; 2940 2941 PetscFunctionBegin; 2942 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2943 PetscFunctionReturn(0); 2944 } 2945 2946 /*@C 2947 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2948 2949 Not Collective 2950 2951 Input Parameter: 2952 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2953 2954 Output Parameter: 2955 . vals - pointer to the data 2956 2957 Level: intermediate 2958 2959 .seealso: MatDenseGetColumn() 2960 @*/ 2961 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2962 { 2963 PetscErrorCode ierr; 2964 2965 PetscFunctionBegin; 2966 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2967 PetscFunctionReturn(0); 2968 } 2969