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