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