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