1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 11 PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian) 12 { 13 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 14 PetscInt j, k, n = A->rmap->n; 15 PetscScalar *v; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 21 if (!hermitian) { 22 for (k=0;k<n;k++) { 23 for (j=k;j<n;j++) { 24 v[j*mat->lda + k] = v[k*mat->lda + j]; 25 } 26 } 27 } else { 28 for (k=0;k<n;k++) { 29 for (j=k;j<n;j++) { 30 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 31 } 32 } 33 } 34 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 39 { 40 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 41 PetscErrorCode ierr; 42 PetscBLASInt info,n; 43 44 PetscFunctionBegin; 45 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 46 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 47 if (A->factortype == MAT_FACTOR_LU) { 48 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 49 if (!mat->fwork) { 50 mat->lfwork = n; 51 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 52 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 53 } 54 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 55 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 56 ierr = PetscFPTrapPop();CHKERRQ(ierr); 57 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 58 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 59 if (A->spd) { 60 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 61 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 62 ierr = PetscFPTrapPop();CHKERRQ(ierr); 63 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 64 #if defined(PETSC_USE_COMPLEX) 65 } else if (A->hermitian) { 66 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 67 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 70 ierr = PetscFPTrapPop();CHKERRQ(ierr); 71 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 72 #endif 73 } else { /* symmetric case */ 74 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 75 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 76 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 77 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 78 ierr = PetscFPTrapPop();CHKERRQ(ierr); 79 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr); 80 } 81 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1); 82 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 84 85 A->ops->solve = NULL; 86 A->ops->matsolve = NULL; 87 A->ops->solvetranspose = NULL; 88 A->ops->matsolvetranspose = NULL; 89 A->ops->solveadd = NULL; 90 A->ops->solvetransposeadd = NULL; 91 A->factortype = MAT_FACTOR_NONE; 92 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 97 { 98 PetscErrorCode ierr; 99 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 100 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 101 PetscScalar *slot,*bb,*v; 102 const PetscScalar *xx; 103 104 PetscFunctionBegin; 105 #if defined(PETSC_USE_DEBUG) 106 for (i=0; i<N; i++) { 107 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 108 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 109 if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n); 110 } 111 #endif 112 if (!N) PetscFunctionReturn(0); 113 114 /* fix right hand side if needed */ 115 if (x && b) { 116 Vec xt; 117 118 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 119 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr); 120 ierr = VecCopy(x,xt);CHKERRQ(ierr); 121 ierr = VecScale(xt,-1.0);CHKERRQ(ierr); 122 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr); 123 ierr = VecDestroy(&xt);CHKERRQ(ierr); 124 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 125 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 126 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 127 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 128 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 129 } 130 131 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 132 for (i=0; i<N; i++) { 133 slot = v + rows[i]*m; 134 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 135 } 136 for (i=0; i<N; i++) { 137 slot = v + rows[i]; 138 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 139 } 140 if (diag != 0.0) { 141 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 142 for (i=0; i<N; i++) { 143 slot = v + (m+1)*rows[i]; 144 *slot = diag; 145 } 146 } 147 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 152 { 153 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 if (c->ptapwork) { 158 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr); 159 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr); 160 } else 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 defined(PETSC_USE_DEBUG) 917 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 918 #endif 919 for (i=0; i<m; i++) { 920 if (indexm[i] < 0) {idx++; continue;} 921 #if defined(PETSC_USE_DEBUG) 922 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 923 #endif 924 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 925 } 926 } 927 } else { 928 for (j=0; j<n; j++) { 929 if (indexn[j] < 0) {idx += m; continue;} 930 #if defined(PETSC_USE_DEBUG) 931 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 932 #endif 933 for (i=0; i<m; i++) { 934 if (indexm[i] < 0) {idx++; continue;} 935 #if defined(PETSC_USE_DEBUG) 936 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 937 #endif 938 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 939 } 940 } 941 } 942 } else { 943 if (addv == INSERT_VALUES) { 944 for (i=0; i<m; i++) { 945 if (indexm[i] < 0) { idx += n; continue;} 946 #if defined(PETSC_USE_DEBUG) 947 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 948 #endif 949 for (j=0; j<n; j++) { 950 if (indexn[j] < 0) { idx++; continue;} 951 #if defined(PETSC_USE_DEBUG) 952 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 953 #endif 954 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 955 } 956 } 957 } else { 958 for (i=0; i<m; i++) { 959 if (indexm[i] < 0) { idx += n; continue;} 960 #if defined(PETSC_USE_DEBUG) 961 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 962 #endif 963 for (j=0; j<n; j++) { 964 if (indexn[j] < 0) { idx++; continue;} 965 #if defined(PETSC_USE_DEBUG) 966 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 967 #endif 968 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 969 } 970 } 971 } 972 } 973 /* hack to prevent unneeded copy to the GPU while returning the array */ 974 #if defined(PETSC_HAVE_CUDA) 975 oldf = A->offloadmask; 976 A->offloadmask = PETSC_OFFLOAD_GPU; 977 #endif 978 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr); 979 #if defined(PETSC_HAVE_CUDA) 980 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 981 #endif 982 PetscFunctionReturn(0); 983 } 984 985 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 986 { 987 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 988 const PetscScalar *vv; 989 PetscInt i,j; 990 PetscErrorCode ierr; 991 992 PetscFunctionBegin; 993 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 994 /* row-oriented output */ 995 for (i=0; i<m; i++) { 996 if (indexm[i] < 0) {v += n;continue;} 997 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); 998 for (j=0; j<n; j++) { 999 if (indexn[j] < 0) {v++; continue;} 1000 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); 1001 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1002 } 1003 } 1004 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1005 PetscFunctionReturn(0); 1006 } 1007 1008 /* -----------------------------------------------------------------*/ 1009 1010 PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer) 1011 { 1012 PetscErrorCode ierr; 1013 PetscBool skipHeader; 1014 PetscViewerFormat format; 1015 PetscInt header[4],M,N,m,lda,i,j,k; 1016 const PetscScalar *v; 1017 PetscScalar *vwork; 1018 1019 PetscFunctionBegin; 1020 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1021 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1022 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1023 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1024 1025 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1026 1027 /* write matrix header */ 1028 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1029 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1030 if (!skipHeader) {ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr);} 1031 1032 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1033 if (format != PETSC_VIEWER_NATIVE) { 1034 PetscInt nnz = m*N, *iwork; 1035 /* store row lengths for each row */ 1036 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr); 1037 for (i=0; i<m; i++) iwork[i] = N; 1038 ierr = PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1039 /* store column indices (zero start index) */ 1040 for (k=0, i=0; i<m; i++) 1041 for (j=0; j<N; j++, k++) 1042 iwork[k] = j; 1043 ierr = PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1044 ierr = PetscFree(iwork);CHKERRQ(ierr); 1045 } 1046 /* store matrix values as a dense matrix in row major order */ 1047 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr); 1048 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr); 1049 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1050 for (k=0, i=0; i<m; i++) 1051 for (j=0; j<N; j++, k++) 1052 vwork[k] = v[i+lda*j]; 1053 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr); 1054 ierr = PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1055 ierr = PetscFree(vwork);CHKERRQ(ierr); 1056 PetscFunctionReturn(0); 1057 } 1058 1059 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer) 1060 { 1061 PetscErrorCode ierr; 1062 PetscBool skipHeader; 1063 PetscInt header[4],M,N,m,nz,lda,i,j,k; 1064 PetscInt rows,cols; 1065 PetscScalar *v,*vwork; 1066 1067 PetscFunctionBegin; 1068 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1069 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1070 1071 if (!skipHeader) { 1072 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1073 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1074 M = header[1]; N = header[2]; 1075 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 1076 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 1077 nz = header[3]; 1078 if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz); 1079 } else { 1080 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1081 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"); 1082 nz = MATRIX_BINARY_FORMAT_DENSE; 1083 } 1084 1085 /* setup global sizes if not set */ 1086 if (mat->rmap->N < 0) mat->rmap->N = M; 1087 if (mat->cmap->N < 0) mat->cmap->N = N; 1088 ierr = MatSetUp(mat);CHKERRQ(ierr); 1089 /* check if global sizes are correct */ 1090 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1091 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); 1092 1093 ierr = MatGetSize(mat,NULL,&N);CHKERRQ(ierr); 1094 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1095 ierr = MatDenseGetArray(mat,&v);CHKERRQ(ierr); 1096 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1097 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1098 PetscInt nnz = m*N; 1099 /* read in matrix values */ 1100 ierr = PetscMalloc1(nnz,&vwork);CHKERRQ(ierr); 1101 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1102 /* store values in column major order */ 1103 for (j=0; j<N; j++) 1104 for (i=0; i<m; i++) 1105 v[i+lda*j] = vwork[i*N+j]; 1106 ierr = PetscFree(vwork);CHKERRQ(ierr); 1107 } else { /* matrix in file is sparse format */ 1108 PetscInt nnz = 0, *rlens, *icols; 1109 /* read in row lengths */ 1110 ierr = PetscMalloc1(m,&rlens);CHKERRQ(ierr); 1111 ierr = PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1112 for (i=0; i<m; i++) nnz += rlens[i]; 1113 /* read in column indices and values */ 1114 ierr = PetscMalloc2(nnz,&icols,nnz,&vwork);CHKERRQ(ierr); 1115 ierr = PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1116 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1117 /* store values in column major order */ 1118 for (k=0, i=0; i<m; i++) 1119 for (j=0; j<rlens[i]; j++, k++) 1120 v[i+lda*icols[k]] = vwork[k]; 1121 ierr = PetscFree(rlens);CHKERRQ(ierr); 1122 ierr = PetscFree2(icols,vwork);CHKERRQ(ierr); 1123 } 1124 ierr = MatDenseRestoreArray(mat,&v);CHKERRQ(ierr); 1125 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1126 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1131 { 1132 PetscBool isbinary, ishdf5; 1133 PetscErrorCode ierr; 1134 1135 PetscFunctionBegin; 1136 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1137 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1138 /* force binary viewer to load .info file if it has not yet done so */ 1139 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1140 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1141 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1142 if (isbinary) { 1143 ierr = MatLoad_Dense_Binary(newMat,viewer);CHKERRQ(ierr); 1144 } else if (ishdf5) { 1145 #if defined(PETSC_HAVE_HDF5) 1146 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1147 #else 1148 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1149 #endif 1150 } else { 1151 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); 1152 } 1153 PetscFunctionReturn(0); 1154 } 1155 1156 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1157 { 1158 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1159 PetscErrorCode ierr; 1160 PetscInt i,j; 1161 const char *name; 1162 PetscScalar *v,*av; 1163 PetscViewerFormat format; 1164 #if defined(PETSC_USE_COMPLEX) 1165 PetscBool allreal = PETSC_TRUE; 1166 #endif 1167 1168 PetscFunctionBegin; 1169 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1170 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1171 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1172 PetscFunctionReturn(0); /* do nothing for now */ 1173 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1174 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1175 for (i=0; i<A->rmap->n; i++) { 1176 v = av + i; 1177 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1178 for (j=0; j<A->cmap->n; j++) { 1179 #if defined(PETSC_USE_COMPLEX) 1180 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1181 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1182 } else if (PetscRealPart(*v)) { 1183 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1184 } 1185 #else 1186 if (*v) { 1187 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1188 } 1189 #endif 1190 v += a->lda; 1191 } 1192 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1193 } 1194 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1195 } else { 1196 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1197 #if defined(PETSC_USE_COMPLEX) 1198 /* determine if matrix has all real values */ 1199 v = av; 1200 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1201 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1202 } 1203 #endif 1204 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1205 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1206 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1207 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1208 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1209 } 1210 1211 for (i=0; i<A->rmap->n; i++) { 1212 v = av + i; 1213 for (j=0; j<A->cmap->n; j++) { 1214 #if defined(PETSC_USE_COMPLEX) 1215 if (allreal) { 1216 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1217 } else { 1218 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1219 } 1220 #else 1221 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1222 #endif 1223 v += a->lda; 1224 } 1225 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1226 } 1227 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1228 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1229 } 1230 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1231 } 1232 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1233 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1234 PetscFunctionReturn(0); 1235 } 1236 1237 static PetscErrorCode MatView_SeqDense_Binary(Mat mat,PetscViewer viewer) 1238 { 1239 PetscErrorCode ierr; 1240 PetscViewerFormat format; 1241 PetscInt header[4],M,N,m,lda,i,j,k; 1242 const PetscScalar *v; 1243 PetscScalar *vwork; 1244 1245 PetscFunctionBegin; 1246 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1247 1248 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1249 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1250 1251 /* write matrix header */ 1252 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1253 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1254 ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr); 1255 1256 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1257 if (format != PETSC_VIEWER_NATIVE) { 1258 PetscInt nnz = m*N, *iwork; 1259 /* store row lengths for each row */ 1260 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr); 1261 for (i=0; i<m; i++) iwork[i] = N; 1262 ierr = PetscViewerBinaryWrite(viewer,iwork,m,PETSC_INT);CHKERRQ(ierr); 1263 /* store column indices (zero start index) */ 1264 for (k=0, i=0; i<m; i++) 1265 for (j=0; j<N; j++, k++) 1266 iwork[k] = j; 1267 ierr = PetscViewerBinaryWrite(viewer,iwork,nnz,PETSC_INT);CHKERRQ(ierr); 1268 ierr = PetscFree(iwork);CHKERRQ(ierr); 1269 } 1270 /* store the matrix values as a dense matrix in row major order */ 1271 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr); 1272 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr); 1273 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1274 for (k=0, i=0; i<m; i++) 1275 for (j=0; j<N; j++, k++) 1276 vwork[k] = v[i+lda*j]; 1277 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr); 1278 ierr = PetscViewerBinaryWrite(viewer,vwork,m*N,PETSC_SCALAR);CHKERRQ(ierr); 1279 ierr = PetscFree(vwork);CHKERRQ(ierr); 1280 PetscFunctionReturn(0); 1281 } 1282 1283 #include <petscdraw.h> 1284 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1285 { 1286 Mat A = (Mat) Aa; 1287 PetscErrorCode ierr; 1288 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1289 int color = PETSC_DRAW_WHITE; 1290 const PetscScalar *v; 1291 PetscViewer viewer; 1292 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1293 PetscViewerFormat format; 1294 1295 PetscFunctionBegin; 1296 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1297 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1298 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1299 1300 /* Loop over matrix elements drawing boxes */ 1301 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 1302 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1303 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1304 /* Blue for negative and Red for positive */ 1305 for (j = 0; j < n; j++) { 1306 x_l = j; x_r = x_l + 1.0; 1307 for (i = 0; i < m; i++) { 1308 y_l = m - i - 1.0; 1309 y_r = y_l + 1.0; 1310 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1311 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1312 else continue; 1313 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1314 } 1315 } 1316 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1317 } else { 1318 /* use contour shading to indicate magnitude of values */ 1319 /* first determine max of all nonzero values */ 1320 PetscReal minv = 0.0, maxv = 0.0; 1321 PetscDraw popup; 1322 1323 for (i=0; i < m*n; i++) { 1324 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1325 } 1326 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1327 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1328 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1329 1330 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1331 for (j=0; j<n; j++) { 1332 x_l = j; 1333 x_r = x_l + 1.0; 1334 for (i=0; i<m; i++) { 1335 y_l = m - i - 1.0; 1336 y_r = y_l + 1.0; 1337 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1338 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1339 } 1340 } 1341 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1342 } 1343 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 1344 PetscFunctionReturn(0); 1345 } 1346 1347 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1348 { 1349 PetscDraw draw; 1350 PetscBool isnull; 1351 PetscReal xr,yr,xl,yl,h,w; 1352 PetscErrorCode ierr; 1353 1354 PetscFunctionBegin; 1355 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1356 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1357 if (isnull) PetscFunctionReturn(0); 1358 1359 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1360 xr += w; yr += h; xl = -w; yl = -h; 1361 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1362 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1363 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1364 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1365 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1366 PetscFunctionReturn(0); 1367 } 1368 1369 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1370 { 1371 PetscErrorCode ierr; 1372 PetscBool iascii,isbinary,isdraw; 1373 1374 PetscFunctionBegin; 1375 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1376 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1377 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1378 1379 if (iascii) { 1380 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1381 } else if (isbinary) { 1382 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1383 } else if (isdraw) { 1384 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1385 } 1386 PetscFunctionReturn(0); 1387 } 1388 1389 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[]) 1390 { 1391 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1392 1393 PetscFunctionBegin; 1394 a->unplacedarray = a->v; 1395 a->unplaced_user_alloc = a->user_alloc; 1396 a->v = (PetscScalar*) array; 1397 #if defined(PETSC_HAVE_CUDA) 1398 A->offloadmask = PETSC_OFFLOAD_CPU; 1399 #endif 1400 PetscFunctionReturn(0); 1401 } 1402 1403 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1404 { 1405 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1406 1407 PetscFunctionBegin; 1408 a->v = a->unplacedarray; 1409 a->user_alloc = a->unplaced_user_alloc; 1410 a->unplacedarray = NULL; 1411 #if defined(PETSC_HAVE_CUDA) 1412 A->offloadmask = PETSC_OFFLOAD_CPU; 1413 #endif 1414 PetscFunctionReturn(0); 1415 } 1416 1417 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1418 { 1419 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1420 PetscErrorCode ierr; 1421 1422 PetscFunctionBegin; 1423 #if defined(PETSC_USE_LOG) 1424 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1425 #endif 1426 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1427 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1428 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1429 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1430 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1431 1432 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1433 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1434 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1435 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1436 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1437 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr); 1440 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1441 #if defined(PETSC_HAVE_ELEMENTAL) 1442 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1443 #endif 1444 #if defined(PETSC_HAVE_CUDA) 1445 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1446 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr); 1447 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr); 1448 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",NULL);CHKERRQ(ierr); 1449 #endif 1450 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1451 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1452 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr); 1453 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1454 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1455 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1456 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1457 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1458 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1459 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1460 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1461 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1464 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_nest_seqdense_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_nest_seqdense_C",NULL);CHKERRQ(ierr); 1469 1470 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1471 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1472 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",NULL);CHKERRQ(ierr); 1474 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1475 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",NULL);CHKERRQ(ierr); 1476 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1477 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",NULL);CHKERRQ(ierr); 1478 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1479 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1484 { 1485 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1486 PetscErrorCode ierr; 1487 PetscInt k,j,m,n,M; 1488 PetscScalar *v,tmp; 1489 1490 PetscFunctionBegin; 1491 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1492 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1493 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1494 for (j=0; j<m; j++) { 1495 for (k=0; k<j; k++) { 1496 tmp = v[j + k*M]; 1497 v[j + k*M] = v[k + j*M]; 1498 v[k + j*M] = tmp; 1499 } 1500 } 1501 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1502 } else { /* out-of-place transpose */ 1503 Mat tmat; 1504 Mat_SeqDense *tmatd; 1505 PetscScalar *v2; 1506 PetscInt M2; 1507 1508 if (reuse != MAT_REUSE_MATRIX) { 1509 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1510 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1511 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1512 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1513 } else tmat = *matout; 1514 1515 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1516 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1517 tmatd = (Mat_SeqDense*)tmat->data; 1518 M2 = tmatd->lda; 1519 for (j=0; j<n; j++) { 1520 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1521 } 1522 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1523 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1524 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1525 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1526 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1527 else { 1528 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1529 } 1530 } 1531 PetscFunctionReturn(0); 1532 } 1533 1534 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1535 { 1536 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1537 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1538 PetscInt i; 1539 const PetscScalar *v1,*v2; 1540 PetscErrorCode ierr; 1541 1542 PetscFunctionBegin; 1543 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1544 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1545 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1546 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1547 for (i=0; i<A1->cmap->n; i++) { 1548 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1549 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1550 v1 += mat1->lda; 1551 v2 += mat2->lda; 1552 } 1553 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1554 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1555 *flg = PETSC_TRUE; 1556 PetscFunctionReturn(0); 1557 } 1558 1559 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1560 { 1561 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1562 PetscInt i,n,len; 1563 PetscScalar *x; 1564 const PetscScalar *vv; 1565 PetscErrorCode ierr; 1566 1567 PetscFunctionBegin; 1568 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1569 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1570 len = PetscMin(A->rmap->n,A->cmap->n); 1571 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1572 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1573 for (i=0; i<len; i++) { 1574 x[i] = vv[i*mat->lda + i]; 1575 } 1576 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1577 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1578 PetscFunctionReturn(0); 1579 } 1580 1581 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1582 { 1583 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1584 const PetscScalar *l,*r; 1585 PetscScalar x,*v,*vv; 1586 PetscErrorCode ierr; 1587 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1588 1589 PetscFunctionBegin; 1590 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1591 if (ll) { 1592 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1593 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1594 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1595 for (i=0; i<m; i++) { 1596 x = l[i]; 1597 v = vv + i; 1598 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1599 } 1600 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1601 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1602 } 1603 if (rr) { 1604 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1605 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1606 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1607 for (i=0; i<n; i++) { 1608 x = r[i]; 1609 v = vv + i*mat->lda; 1610 for (j=0; j<m; j++) (*v++) *= x; 1611 } 1612 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1613 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1614 } 1615 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1616 PetscFunctionReturn(0); 1617 } 1618 1619 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1620 { 1621 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1622 PetscScalar *v,*vv; 1623 PetscReal sum = 0.0; 1624 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1625 PetscErrorCode ierr; 1626 1627 PetscFunctionBegin; 1628 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1629 v = vv; 1630 if (type == NORM_FROBENIUS) { 1631 if (lda>m) { 1632 for (j=0; j<A->cmap->n; j++) { 1633 v = vv+j*lda; 1634 for (i=0; i<m; i++) { 1635 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1636 } 1637 } 1638 } else { 1639 #if defined(PETSC_USE_REAL___FP16) 1640 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1641 *nrm = BLASnrm2_(&cnt,v,&one); 1642 } 1643 #else 1644 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1645 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1646 } 1647 } 1648 *nrm = PetscSqrtReal(sum); 1649 #endif 1650 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1651 } else if (type == NORM_1) { 1652 *nrm = 0.0; 1653 for (j=0; j<A->cmap->n; j++) { 1654 v = vv + j*mat->lda; 1655 sum = 0.0; 1656 for (i=0; i<A->rmap->n; i++) { 1657 sum += PetscAbsScalar(*v); v++; 1658 } 1659 if (sum > *nrm) *nrm = sum; 1660 } 1661 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1662 } else if (type == NORM_INFINITY) { 1663 *nrm = 0.0; 1664 for (j=0; j<A->rmap->n; j++) { 1665 v = vv + j; 1666 sum = 0.0; 1667 for (i=0; i<A->cmap->n; i++) { 1668 sum += PetscAbsScalar(*v); v += mat->lda; 1669 } 1670 if (sum > *nrm) *nrm = sum; 1671 } 1672 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1673 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1674 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1675 PetscFunctionReturn(0); 1676 } 1677 1678 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1679 { 1680 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1681 PetscErrorCode ierr; 1682 1683 PetscFunctionBegin; 1684 switch (op) { 1685 case MAT_ROW_ORIENTED: 1686 aij->roworiented = flg; 1687 break; 1688 case MAT_NEW_NONZERO_LOCATIONS: 1689 case MAT_NEW_NONZERO_LOCATION_ERR: 1690 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1691 case MAT_NEW_DIAGONALS: 1692 case MAT_KEEP_NONZERO_PATTERN: 1693 case MAT_IGNORE_OFF_PROC_ENTRIES: 1694 case MAT_USE_HASH_TABLE: 1695 case MAT_IGNORE_ZERO_ENTRIES: 1696 case MAT_IGNORE_LOWER_TRIANGULAR: 1697 case MAT_SORTED_FULL: 1698 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1699 break; 1700 case MAT_SPD: 1701 case MAT_SYMMETRIC: 1702 case MAT_STRUCTURALLY_SYMMETRIC: 1703 case MAT_HERMITIAN: 1704 case MAT_SYMMETRY_ETERNAL: 1705 /* These options are handled directly by MatSetOption() */ 1706 break; 1707 default: 1708 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1709 } 1710 PetscFunctionReturn(0); 1711 } 1712 1713 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1714 { 1715 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1716 PetscErrorCode ierr; 1717 PetscInt lda=l->lda,m=A->rmap->n,j; 1718 PetscScalar *v; 1719 1720 PetscFunctionBegin; 1721 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1722 if (lda>m) { 1723 for (j=0; j<A->cmap->n; j++) { 1724 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1725 } 1726 } else { 1727 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1728 } 1729 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1730 PetscFunctionReturn(0); 1731 } 1732 1733 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1734 { 1735 PetscErrorCode ierr; 1736 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1737 PetscInt m = l->lda, n = A->cmap->n, i,j; 1738 PetscScalar *slot,*bb,*v; 1739 const PetscScalar *xx; 1740 1741 PetscFunctionBegin; 1742 #if defined(PETSC_USE_DEBUG) 1743 for (i=0; i<N; i++) { 1744 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1745 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); 1746 } 1747 #endif 1748 if (!N) PetscFunctionReturn(0); 1749 1750 /* fix right hand side if needed */ 1751 if (x && b) { 1752 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1753 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1754 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1755 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1756 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1757 } 1758 1759 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1760 for (i=0; i<N; i++) { 1761 slot = v + rows[i]; 1762 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1763 } 1764 if (diag != 0.0) { 1765 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1766 for (i=0; i<N; i++) { 1767 slot = v + (m+1)*rows[i]; 1768 *slot = diag; 1769 } 1770 } 1771 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1772 PetscFunctionReturn(0); 1773 } 1774 1775 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1776 { 1777 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1778 1779 PetscFunctionBegin; 1780 *lda = mat->lda; 1781 PetscFunctionReturn(0); 1782 } 1783 1784 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1785 { 1786 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1787 1788 PetscFunctionBegin; 1789 *array = mat->v; 1790 PetscFunctionReturn(0); 1791 } 1792 1793 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1794 { 1795 PetscFunctionBegin; 1796 PetscFunctionReturn(0); 1797 } 1798 1799 /*@C 1800 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1801 1802 Logically Collective on Mat 1803 1804 Input Parameter: 1805 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1806 1807 Output Parameter: 1808 . lda - the leading dimension 1809 1810 Level: intermediate 1811 1812 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1813 @*/ 1814 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1815 { 1816 PetscErrorCode ierr; 1817 1818 PetscFunctionBegin; 1819 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1820 PetscFunctionReturn(0); 1821 } 1822 1823 /*@C 1824 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1825 1826 Logically Collective on Mat 1827 1828 Input Parameter: 1829 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1830 1831 Output Parameter: 1832 . array - pointer to the data 1833 1834 Level: intermediate 1835 1836 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1837 @*/ 1838 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1839 { 1840 PetscErrorCode ierr; 1841 1842 PetscFunctionBegin; 1843 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1844 PetscFunctionReturn(0); 1845 } 1846 1847 /*@C 1848 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1849 1850 Logically Collective on Mat 1851 1852 Input Parameters: 1853 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1854 - array - pointer to the data 1855 1856 Level: intermediate 1857 1858 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1859 @*/ 1860 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1861 { 1862 PetscErrorCode ierr; 1863 1864 PetscFunctionBegin; 1865 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1866 if (array) *array = NULL; 1867 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1868 PetscFunctionReturn(0); 1869 } 1870 1871 /*@C 1872 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1873 1874 Not Collective 1875 1876 Input Parameter: 1877 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1878 1879 Output Parameter: 1880 . array - pointer to the data 1881 1882 Level: intermediate 1883 1884 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1885 @*/ 1886 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1887 { 1888 PetscErrorCode ierr; 1889 1890 PetscFunctionBegin; 1891 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1892 PetscFunctionReturn(0); 1893 } 1894 1895 /*@C 1896 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1897 1898 Not Collective 1899 1900 Input Parameters: 1901 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1902 - array - pointer to the data 1903 1904 Level: intermediate 1905 1906 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1907 @*/ 1908 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1909 { 1910 PetscErrorCode ierr; 1911 1912 PetscFunctionBegin; 1913 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1914 if (array) *array = NULL; 1915 PetscFunctionReturn(0); 1916 } 1917 1918 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1919 { 1920 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1921 PetscErrorCode ierr; 1922 PetscInt i,j,nrows,ncols,blda; 1923 const PetscInt *irow,*icol; 1924 PetscScalar *av,*bv,*v = mat->v; 1925 Mat newmat; 1926 1927 PetscFunctionBegin; 1928 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1929 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1930 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1931 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1932 1933 /* Check submatrixcall */ 1934 if (scall == MAT_REUSE_MATRIX) { 1935 PetscInt n_cols,n_rows; 1936 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1937 if (n_rows != nrows || n_cols != ncols) { 1938 /* resize the result matrix to match number of requested rows/columns */ 1939 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1940 } 1941 newmat = *B; 1942 } else { 1943 /* Create and fill new matrix */ 1944 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1945 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1946 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1947 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1948 } 1949 1950 /* Now extract the data pointers and do the copy,column at a time */ 1951 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 1952 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 1953 for (i=0; i<ncols; i++) { 1954 av = v + mat->lda*icol[i]; 1955 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 1956 bv += blda; 1957 } 1958 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 1959 1960 /* Assemble the matrices so that the correct flags are set */ 1961 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1962 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1963 1964 /* Free work space */ 1965 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1966 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1967 *B = newmat; 1968 PetscFunctionReturn(0); 1969 } 1970 1971 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1972 { 1973 PetscErrorCode ierr; 1974 PetscInt i; 1975 1976 PetscFunctionBegin; 1977 if (scall == MAT_INITIAL_MATRIX) { 1978 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1979 } 1980 1981 for (i=0; i<n; i++) { 1982 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1983 } 1984 PetscFunctionReturn(0); 1985 } 1986 1987 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1988 { 1989 PetscFunctionBegin; 1990 PetscFunctionReturn(0); 1991 } 1992 1993 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1994 { 1995 PetscFunctionBegin; 1996 PetscFunctionReturn(0); 1997 } 1998 1999 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2000 { 2001 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2002 PetscErrorCode ierr; 2003 const PetscScalar *va; 2004 PetscScalar *vb; 2005 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2006 2007 PetscFunctionBegin; 2008 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2009 if (A->ops->copy != B->ops->copy) { 2010 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2011 PetscFunctionReturn(0); 2012 } 2013 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2014 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2015 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2016 if (lda1>m || lda2>m) { 2017 for (j=0; j<n; j++) { 2018 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2019 } 2020 } else { 2021 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2022 } 2023 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2024 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2025 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2026 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2027 PetscFunctionReturn(0); 2028 } 2029 2030 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2031 { 2032 PetscErrorCode ierr; 2033 2034 PetscFunctionBegin; 2035 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2036 PetscFunctionReturn(0); 2037 } 2038 2039 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2040 { 2041 PetscInt i,nz = A->rmap->n*A->cmap->n; 2042 PetscScalar *aa; 2043 PetscErrorCode ierr; 2044 2045 PetscFunctionBegin; 2046 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2047 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2048 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2049 PetscFunctionReturn(0); 2050 } 2051 2052 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2053 { 2054 PetscInt i,nz = A->rmap->n*A->cmap->n; 2055 PetscScalar *aa; 2056 PetscErrorCode ierr; 2057 2058 PetscFunctionBegin; 2059 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2060 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2061 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2062 PetscFunctionReturn(0); 2063 } 2064 2065 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2066 { 2067 PetscInt i,nz = A->rmap->n*A->cmap->n; 2068 PetscScalar *aa; 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2073 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2074 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2075 PetscFunctionReturn(0); 2076 } 2077 2078 /* ----------------------------------------------------------------*/ 2079 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2080 { 2081 PetscErrorCode ierr; 2082 PetscInt m=A->rmap->n,n=B->cmap->n; 2083 PetscBool flg; 2084 2085 PetscFunctionBegin; 2086 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2087 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2088 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2089 ierr = MatSeqDenseSetPreallocation(C,NULL);CHKERRQ(ierr); 2090 PetscFunctionReturn(0); 2091 } 2092 2093 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2094 { 2095 Mat_SeqDense *a,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2096 PetscBLASInt m,n,k; 2097 const PetscScalar *av,*bv; 2098 PetscScalar *cv; 2099 PetscScalar _DOne=1.0,_DZero=0.0; 2100 PetscBool flg; 2101 PetscErrorCode (*numeric)(Mat,Mat,Mat)=NULL; 2102 PetscErrorCode ierr; 2103 2104 PetscFunctionBegin; 2105 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2106 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2107 if (flg) numeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2108 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&flg);CHKERRQ(ierr); 2109 if (flg) numeric = MatMatMultNumeric_SeqBAIJ_SeqDense; 2110 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&flg);CHKERRQ(ierr); 2111 if (flg) numeric = MatMatMultNumeric_SeqSBAIJ_SeqDense; 2112 ierr = PetscObjectTypeCompare((PetscObject)A,MATNEST,&flg);CHKERRQ(ierr); 2113 if (flg) numeric = MatMatMultNumeric_Nest_Dense; 2114 if (numeric) { 2115 C->ops->matmultnumeric = numeric; 2116 ierr = (*numeric)(A,B,C);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 a = (Mat_SeqDense*)A->data; 2120 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2121 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2122 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2123 if (!m || !n || !k) PetscFunctionReturn(0); 2124 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2125 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2126 ierr = MatDenseGetArray(C,&cv);CHKERRQ(ierr); 2127 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2128 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2129 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2130 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2131 ierr = MatDenseRestoreArray(C,&cv);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2136 { 2137 PetscErrorCode ierr; 2138 PetscInt m=A->rmap->n,n=B->rmap->n; 2139 PetscBool flg; 2140 2141 PetscFunctionBegin; 2142 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2143 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2144 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2145 ierr = MatSeqDenseSetPreallocation(C,NULL);CHKERRQ(ierr); 2146 C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqDense_SeqDense; 2147 PetscFunctionReturn(0); 2148 } 2149 2150 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2151 { 2152 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2153 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2154 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2155 PetscBLASInt m,n,k; 2156 PetscScalar _DOne=1.0,_DZero=0.0; 2157 PetscErrorCode ierr; 2158 2159 PetscFunctionBegin; 2160 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2161 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2162 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2163 if (!m || !n || !k) PetscFunctionReturn(0); 2164 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2165 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2166 PetscFunctionReturn(0); 2167 } 2168 2169 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2170 { 2171 PetscErrorCode ierr; 2172 PetscInt m=A->cmap->n,n=B->cmap->n; 2173 PetscBool flg; 2174 2175 PetscFunctionBegin; 2176 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2177 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2178 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2179 ierr = MatSeqDenseSetPreallocation(C,NULL);CHKERRQ(ierr); 2180 PetscFunctionReturn(0); 2181 } 2182 2183 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2184 { 2185 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2186 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2187 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2188 PetscBLASInt m,n,k; 2189 PetscScalar _DOne=1.0,_DZero=0.0; 2190 PetscErrorCode ierr; 2191 2192 PetscFunctionBegin; 2193 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2194 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2195 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2196 if (!m || !n || !k) PetscFunctionReturn(0); 2197 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2198 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2199 PetscFunctionReturn(0); 2200 } 2201 2202 /* ----------------------------------------------- */ 2203 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2204 { 2205 PetscFunctionBegin; 2206 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2207 C->ops->productsymbolic = MatProductSymbolic_AB; 2208 /* dense mat may not call MatProductSymbolic(), thus set C->ops->productnumeric here */ 2209 C->ops->productnumeric = MatProductNumeric_AB; 2210 PetscFunctionReturn(0); 2211 } 2212 2213 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2214 { 2215 PetscFunctionBegin; 2216 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2217 C->ops->productsymbolic = MatProductSymbolic_AtB; 2218 C->ops->productnumeric = MatProductNumeric_AtB; 2219 PetscFunctionReturn(0); 2220 } 2221 2222 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2223 { 2224 PetscFunctionBegin; 2225 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2226 C->ops->productsymbolic = MatProductSymbolic_ABt; 2227 C->ops->productnumeric = MatProductNumeric_ABt; 2228 PetscFunctionReturn(0); 2229 } 2230 2231 static PetscErrorCode MatProductSetFromOptions_SeqDense_PtAP(Mat C) 2232 { 2233 PetscFunctionBegin; 2234 C->ops->productsymbolic = MatProductSymbolic_Basic; 2235 PetscFunctionReturn(0); 2236 } 2237 2238 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2239 { 2240 PetscErrorCode ierr; 2241 Mat_Product *product = C->product; 2242 2243 PetscFunctionBegin; 2244 switch (product->type) { 2245 case MATPRODUCT_AB: 2246 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2247 break; 2248 case MATPRODUCT_AtB: 2249 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2250 break; 2251 case MATPRODUCT_ABt: 2252 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2253 break; 2254 case MATPRODUCT_PtAP: 2255 ierr = MatProductSetFromOptions_SeqDense_PtAP(C);CHKERRQ(ierr); 2256 break; 2257 default: SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type is not supported"); 2258 } 2259 PetscFunctionReturn(0); 2260 } 2261 /* ----------------------------------------------- */ 2262 2263 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2264 { 2265 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2266 PetscErrorCode ierr; 2267 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2268 PetscScalar *x; 2269 const PetscScalar *aa; 2270 2271 PetscFunctionBegin; 2272 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2273 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2274 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2275 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2276 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2277 for (i=0; i<m; i++) { 2278 x[i] = aa[i]; if (idx) idx[i] = 0; 2279 for (j=1; j<n; j++) { 2280 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2281 } 2282 } 2283 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2284 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2285 PetscFunctionReturn(0); 2286 } 2287 2288 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2289 { 2290 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2291 PetscErrorCode ierr; 2292 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2293 PetscScalar *x; 2294 PetscReal atmp; 2295 const PetscScalar *aa; 2296 2297 PetscFunctionBegin; 2298 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2299 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2300 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2301 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2302 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2303 for (i=0; i<m; i++) { 2304 x[i] = PetscAbsScalar(aa[i]); 2305 for (j=1; j<n; j++) { 2306 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2307 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2308 } 2309 } 2310 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2311 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 2315 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2316 { 2317 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2318 PetscErrorCode ierr; 2319 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2320 PetscScalar *x; 2321 const PetscScalar *aa; 2322 2323 PetscFunctionBegin; 2324 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2325 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2326 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2327 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2328 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2329 for (i=0; i<m; i++) { 2330 x[i] = aa[i]; if (idx) idx[i] = 0; 2331 for (j=1; j<n; j++) { 2332 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2333 } 2334 } 2335 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2336 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2337 PetscFunctionReturn(0); 2338 } 2339 2340 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2341 { 2342 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2343 PetscErrorCode ierr; 2344 PetscScalar *x; 2345 const PetscScalar *aa; 2346 2347 PetscFunctionBegin; 2348 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2349 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2350 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2351 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2352 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2353 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2358 { 2359 PetscErrorCode ierr; 2360 PetscInt i,j,m,n; 2361 const PetscScalar *a; 2362 2363 PetscFunctionBegin; 2364 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2365 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2366 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2367 if (type == NORM_2) { 2368 for (i=0; i<n; i++) { 2369 for (j=0; j<m; j++) { 2370 norms[i] += PetscAbsScalar(a[j]*a[j]); 2371 } 2372 a += m; 2373 } 2374 } else if (type == NORM_1) { 2375 for (i=0; i<n; i++) { 2376 for (j=0; j<m; j++) { 2377 norms[i] += PetscAbsScalar(a[j]); 2378 } 2379 a += m; 2380 } 2381 } else if (type == NORM_INFINITY) { 2382 for (i=0; i<n; i++) { 2383 for (j=0; j<m; j++) { 2384 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2385 } 2386 a += m; 2387 } 2388 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2389 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2390 if (type == NORM_2) { 2391 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2392 } 2393 PetscFunctionReturn(0); 2394 } 2395 2396 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2397 { 2398 PetscErrorCode ierr; 2399 PetscScalar *a; 2400 PetscInt m,n,i; 2401 2402 PetscFunctionBegin; 2403 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2404 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2405 for (i=0; i<m*n; i++) { 2406 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2407 } 2408 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2409 PetscFunctionReturn(0); 2410 } 2411 2412 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2413 { 2414 PetscFunctionBegin; 2415 *missing = PETSC_FALSE; 2416 PetscFunctionReturn(0); 2417 } 2418 2419 /* vals is not const */ 2420 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2421 { 2422 PetscErrorCode ierr; 2423 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2424 PetscScalar *v; 2425 2426 PetscFunctionBegin; 2427 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2428 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2429 *vals = v+col*a->lda; 2430 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2431 PetscFunctionReturn(0); 2432 } 2433 2434 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2435 { 2436 PetscFunctionBegin; 2437 *vals = 0; /* user cannot accidently use the array later */ 2438 PetscFunctionReturn(0); 2439 } 2440 2441 /* -------------------------------------------------------------------*/ 2442 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2443 MatGetRow_SeqDense, 2444 MatRestoreRow_SeqDense, 2445 MatMult_SeqDense, 2446 /* 4*/ MatMultAdd_SeqDense, 2447 MatMultTranspose_SeqDense, 2448 MatMultTransposeAdd_SeqDense, 2449 0, 2450 0, 2451 0, 2452 /* 10*/ 0, 2453 MatLUFactor_SeqDense, 2454 MatCholeskyFactor_SeqDense, 2455 MatSOR_SeqDense, 2456 MatTranspose_SeqDense, 2457 /* 15*/ MatGetInfo_SeqDense, 2458 MatEqual_SeqDense, 2459 MatGetDiagonal_SeqDense, 2460 MatDiagonalScale_SeqDense, 2461 MatNorm_SeqDense, 2462 /* 20*/ MatAssemblyBegin_SeqDense, 2463 MatAssemblyEnd_SeqDense, 2464 MatSetOption_SeqDense, 2465 MatZeroEntries_SeqDense, 2466 /* 24*/ MatZeroRows_SeqDense, 2467 0, 2468 0, 2469 0, 2470 0, 2471 /* 29*/ MatSetUp_SeqDense, 2472 0, 2473 0, 2474 0, 2475 0, 2476 /* 34*/ MatDuplicate_SeqDense, 2477 0, 2478 0, 2479 0, 2480 0, 2481 /* 39*/ MatAXPY_SeqDense, 2482 MatCreateSubMatrices_SeqDense, 2483 0, 2484 MatGetValues_SeqDense, 2485 MatCopy_SeqDense, 2486 /* 44*/ MatGetRowMax_SeqDense, 2487 MatScale_SeqDense, 2488 MatShift_Basic, 2489 0, 2490 MatZeroRowsColumns_SeqDense, 2491 /* 49*/ MatSetRandom_SeqDense, 2492 0, 2493 0, 2494 0, 2495 0, 2496 /* 54*/ 0, 2497 0, 2498 0, 2499 0, 2500 0, 2501 /* 59*/ 0, 2502 MatDestroy_SeqDense, 2503 MatView_SeqDense, 2504 0, 2505 0, 2506 /* 64*/ 0, 2507 0, 2508 0, 2509 0, 2510 0, 2511 /* 69*/ MatGetRowMaxAbs_SeqDense, 2512 0, 2513 0, 2514 0, 2515 0, 2516 /* 74*/ 0, 2517 0, 2518 0, 2519 0, 2520 0, 2521 /* 79*/ 0, 2522 0, 2523 0, 2524 0, 2525 /* 83*/ MatLoad_SeqDense, 2526 0, 2527 MatIsHermitian_SeqDense, 2528 0, 2529 0, 2530 0, 2531 /* 89*/ 0, 2532 0, 2533 MatMatMultNumeric_SeqDense_SeqDense, 2534 0, 2535 0, 2536 /* 94*/ 0, 2537 0, 2538 0, 2539 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2540 0, 2541 /* 99*/ MatProductSetFromOptions_SeqDense, 2542 0, 2543 0, 2544 MatConjugate_SeqDense, 2545 0, 2546 /*104*/ 0, 2547 MatRealPart_SeqDense, 2548 MatImaginaryPart_SeqDense, 2549 0, 2550 0, 2551 /*109*/ 0, 2552 0, 2553 MatGetRowMin_SeqDense, 2554 MatGetColumnVector_SeqDense, 2555 MatMissingDiagonal_SeqDense, 2556 /*114*/ 0, 2557 0, 2558 0, 2559 0, 2560 0, 2561 /*119*/ 0, 2562 0, 2563 0, 2564 0, 2565 0, 2566 /*124*/ 0, 2567 MatGetColumnNorms_SeqDense, 2568 0, 2569 0, 2570 0, 2571 /*129*/ 0, 2572 0, 2573 0, 2574 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2575 0, 2576 /*134*/ 0, 2577 0, 2578 0, 2579 0, 2580 0, 2581 /*139*/ 0, 2582 0, 2583 0, 2584 0, 2585 0, 2586 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2587 /*145*/ 0, 2588 0, 2589 0 2590 }; 2591 2592 /*@C 2593 MatCreateSeqDense - Creates a sequential dense matrix that 2594 is stored in column major order (the usual Fortran 77 manner). Many 2595 of the matrix operations use the BLAS and LAPACK routines. 2596 2597 Collective 2598 2599 Input Parameters: 2600 + comm - MPI communicator, set to PETSC_COMM_SELF 2601 . m - number of rows 2602 . n - number of columns 2603 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2604 to control all matrix memory allocation. 2605 2606 Output Parameter: 2607 . A - the matrix 2608 2609 Notes: 2610 The data input variable is intended primarily for Fortran programmers 2611 who wish to allocate their own matrix memory space. Most users should 2612 set data=NULL. 2613 2614 Level: intermediate 2615 2616 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2617 @*/ 2618 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2619 { 2620 PetscErrorCode ierr; 2621 2622 PetscFunctionBegin; 2623 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2624 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2625 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2626 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2627 PetscFunctionReturn(0); 2628 } 2629 2630 /*@C 2631 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2632 2633 Collective 2634 2635 Input Parameters: 2636 + B - the matrix 2637 - data - the array (or NULL) 2638 2639 Notes: 2640 The data input variable is intended primarily for Fortran programmers 2641 who wish to allocate their own matrix memory space. Most users should 2642 need not call this routine. 2643 2644 Level: intermediate 2645 2646 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2647 2648 @*/ 2649 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2650 { 2651 PetscErrorCode ierr; 2652 2653 PetscFunctionBegin; 2654 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2655 PetscFunctionReturn(0); 2656 } 2657 2658 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2659 { 2660 Mat_SeqDense *b; 2661 PetscErrorCode ierr; 2662 2663 PetscFunctionBegin; 2664 B->preallocated = PETSC_TRUE; 2665 2666 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2667 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2668 2669 b = (Mat_SeqDense*)B->data; 2670 b->Mmax = B->rmap->n; 2671 b->Nmax = B->cmap->n; 2672 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2673 2674 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2675 if (!data) { /* petsc-allocated storage */ 2676 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2677 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2678 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2679 2680 b->user_alloc = PETSC_FALSE; 2681 } else { /* user-allocated storage */ 2682 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2683 b->v = data; 2684 b->user_alloc = PETSC_TRUE; 2685 } 2686 B->assembled = PETSC_TRUE; 2687 PetscFunctionReturn(0); 2688 } 2689 2690 #if defined(PETSC_HAVE_ELEMENTAL) 2691 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2692 { 2693 Mat mat_elemental; 2694 PetscErrorCode ierr; 2695 const PetscScalar *array; 2696 PetscScalar *v_colwise; 2697 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2698 2699 PetscFunctionBegin; 2700 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2701 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2702 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2703 k = 0; 2704 for (j=0; j<N; j++) { 2705 cols[j] = j; 2706 for (i=0; i<M; i++) { 2707 v_colwise[j*M+i] = array[k++]; 2708 } 2709 } 2710 for (i=0; i<M; i++) { 2711 rows[i] = i; 2712 } 2713 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2714 2715 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2716 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2717 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2718 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2719 2720 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2721 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2722 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2723 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2724 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2725 2726 if (reuse == MAT_INPLACE_MATRIX) { 2727 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2728 } else { 2729 *newmat = mat_elemental; 2730 } 2731 PetscFunctionReturn(0); 2732 } 2733 #endif 2734 2735 /*@C 2736 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2737 2738 Input parameter: 2739 + A - the matrix 2740 - lda - the leading dimension 2741 2742 Notes: 2743 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2744 it asserts that the preallocation has a leading dimension (the LDA parameter 2745 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2746 2747 Level: intermediate 2748 2749 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2750 2751 @*/ 2752 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2753 { 2754 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2755 2756 PetscFunctionBegin; 2757 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); 2758 b->lda = lda; 2759 b->changelda = PETSC_FALSE; 2760 b->Mmax = PetscMax(b->Mmax,lda); 2761 PetscFunctionReturn(0); 2762 } 2763 2764 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2765 { 2766 PetscErrorCode ierr; 2767 PetscMPIInt size; 2768 2769 PetscFunctionBegin; 2770 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2771 if (size == 1) { 2772 if (scall == MAT_INITIAL_MATRIX) { 2773 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2774 } else { 2775 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2776 } 2777 } else { 2778 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2779 } 2780 PetscFunctionReturn(0); 2781 } 2782 2783 /*MC 2784 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2785 2786 Options Database Keys: 2787 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2788 2789 Level: beginner 2790 2791 .seealso: MatCreateSeqDense() 2792 2793 M*/ 2794 PetscErrorCode MatCreate_SeqDense(Mat B) 2795 { 2796 Mat_SeqDense *b; 2797 PetscErrorCode ierr; 2798 PetscMPIInt size; 2799 2800 PetscFunctionBegin; 2801 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2802 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2803 2804 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2805 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2806 B->data = (void*)b; 2807 2808 b->roworiented = PETSC_TRUE; 2809 2810 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2811 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2812 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2813 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2814 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2815 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2816 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2817 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2818 #if defined(PETSC_HAVE_ELEMENTAL) 2819 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2820 #endif 2821 #if defined(PETSC_HAVE_CUDA) 2822 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 2823 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2824 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2825 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2826 #endif 2827 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2828 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 2829 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 2830 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2831 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2832 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2833 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqbaij_seqdense_C",MatMatMultSymbolic_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2834 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqbaij_seqdense_C",MatMatMultNumeric_SeqBAIJ_SeqDense);CHKERRQ(ierr); 2835 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 2836 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqsbaij_seqdense_C",MatMatMultSymbolic_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2837 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqsbaij_seqdense_C",MatMatMultNumeric_SeqSBAIJ_SeqDense);CHKERRQ(ierr); 2838 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2839 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2840 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2841 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2842 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2843 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2844 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_nest_seqdense_C",MatMatMultSymbolic_Nest_Dense);CHKERRQ(ierr); 2845 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_nest_seqdense_C",MatMatMultNumeric_Nest_Dense);CHKERRQ(ierr); 2846 2847 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2848 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2849 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2850 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2851 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2852 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2853 2854 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2855 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2856 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2857 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2858 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2859 PetscFunctionReturn(0); 2860 } 2861 2862 /*@C 2863 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. 2864 2865 Not Collective 2866 2867 Input Parameter: 2868 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2869 - col - column index 2870 2871 Output Parameter: 2872 . vals - pointer to the data 2873 2874 Level: intermediate 2875 2876 .seealso: MatDenseRestoreColumn() 2877 @*/ 2878 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2879 { 2880 PetscErrorCode ierr; 2881 2882 PetscFunctionBegin; 2883 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2884 PetscFunctionReturn(0); 2885 } 2886 2887 /*@C 2888 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2889 2890 Not Collective 2891 2892 Input Parameter: 2893 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2894 2895 Output Parameter: 2896 . vals - pointer to the data 2897 2898 Level: intermediate 2899 2900 .seealso: MatDenseGetColumn() 2901 @*/ 2902 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2903 { 2904 PetscErrorCode ierr; 2905 2906 PetscFunctionBegin; 2907 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2908 PetscFunctionReturn(0); 2909 } 2910