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