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