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 = 0,lday = 0,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.0*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 = 0; 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,NULL,NULL,&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,NULL,&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 || ftype == MAT_FACTOR_ILU) { 749 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 750 (*fact)->ops->ilufactorsymbolic = MatLUFactorSymbolic_SeqDense; 751 } else { 752 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 753 } 754 (*fact)->factortype = ftype; 755 756 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 757 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 758 PetscFunctionReturn(0); 759 } 760 761 /* ------------------------------------------------------------------*/ 762 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 763 { 764 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 765 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 766 const PetscScalar *b; 767 PetscErrorCode ierr; 768 PetscInt m = A->rmap->n,i; 769 PetscBLASInt o = 1,bm = 0; 770 771 PetscFunctionBegin; 772 #if defined(PETSC_HAVE_CUDA) 773 if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented"); 774 #endif 775 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 776 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 777 if (flag & SOR_ZERO_INITIAL_GUESS) { 778 /* this is a hack fix, should have another version without the second BLASdotu */ 779 ierr = VecSet(xx,zero);CHKERRQ(ierr); 780 } 781 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 782 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 783 its = its*lits; 784 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 785 while (its--) { 786 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 787 for (i=0; i<m; i++) { 788 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 789 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 790 } 791 } 792 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 793 for (i=m-1; i>=0; i--) { 794 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 795 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 796 } 797 } 798 } 799 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 800 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 /* -----------------------------------------------------------------*/ 805 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 806 { 807 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 808 const PetscScalar *v = mat->v,*x; 809 PetscScalar *y; 810 PetscErrorCode ierr; 811 PetscBLASInt m, n,_One=1; 812 PetscScalar _DOne=1.0,_DZero=0.0; 813 814 PetscFunctionBegin; 815 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 816 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 817 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 818 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 819 if (!A->rmap->n || !A->cmap->n) { 820 PetscBLASInt i; 821 for (i=0; i<n; i++) y[i] = 0.0; 822 } else { 823 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 824 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 825 } 826 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 827 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 828 PetscFunctionReturn(0); 829 } 830 831 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 832 { 833 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 834 PetscScalar *y,_DOne=1.0,_DZero=0.0; 835 PetscErrorCode ierr; 836 PetscBLASInt m, n, _One=1; 837 const PetscScalar *v = mat->v,*x; 838 839 PetscFunctionBegin; 840 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 841 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 842 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 843 ierr = VecGetArrayWrite(yy,&y);CHKERRQ(ierr); 844 if (!A->rmap->n || !A->cmap->n) { 845 PetscBLASInt i; 846 for (i=0; i<m; i++) y[i] = 0.0; 847 } else { 848 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 849 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 850 } 851 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 852 ierr = VecRestoreArrayWrite(yy,&y);CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 857 { 858 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 859 const PetscScalar *v = mat->v,*x; 860 PetscScalar *y,_DOne=1.0; 861 PetscErrorCode ierr; 862 PetscBLASInt m, n, _One=1; 863 864 PetscFunctionBegin; 865 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 866 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 867 ierr = VecCopy(zz,yy);CHKERRQ(ierr); 868 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 869 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 870 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 871 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 872 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 873 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 874 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 875 PetscFunctionReturn(0); 876 } 877 878 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 879 { 880 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 881 const PetscScalar *v = mat->v,*x; 882 PetscScalar *y; 883 PetscErrorCode ierr; 884 PetscBLASInt m, n, _One=1; 885 PetscScalar _DOne=1.0; 886 887 PetscFunctionBegin; 888 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 889 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 890 ierr = VecCopy(zz,yy);CHKERRQ(ierr); 891 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 892 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 893 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 894 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 895 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 896 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 897 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 898 PetscFunctionReturn(0); 899 } 900 901 /* -----------------------------------------------------------------*/ 902 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 903 { 904 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 905 PetscErrorCode ierr; 906 PetscInt i; 907 908 PetscFunctionBegin; 909 *ncols = A->cmap->n; 910 if (cols) { 911 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 912 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 913 } 914 if (vals) { 915 const PetscScalar *v; 916 917 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 918 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 919 v += row; 920 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 921 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 922 } 923 PetscFunctionReturn(0); 924 } 925 926 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 927 { 928 PetscErrorCode ierr; 929 930 PetscFunctionBegin; 931 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 932 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 933 PetscFunctionReturn(0); 934 } 935 /* ----------------------------------------------------------------*/ 936 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 937 { 938 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 939 PetscScalar *av; 940 PetscInt i,j,idx=0; 941 #if defined(PETSC_HAVE_CUDA) 942 PetscOffloadMask oldf; 943 #endif 944 PetscErrorCode ierr; 945 946 PetscFunctionBegin; 947 ierr = MatDenseGetArray(A,&av);CHKERRQ(ierr); 948 if (!mat->roworiented) { 949 if (addv == INSERT_VALUES) { 950 for (j=0; j<n; j++) { 951 if (indexn[j] < 0) {idx += m; continue;} 952 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); 953 for (i=0; i<m; i++) { 954 if (indexm[i] < 0) {idx++; continue;} 955 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); 956 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 957 } 958 } 959 } else { 960 for (j=0; j<n; j++) { 961 if (indexn[j] < 0) {idx += m; continue;} 962 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); 963 for (i=0; i<m; i++) { 964 if (indexm[i] < 0) {idx++; continue;} 965 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); 966 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 967 } 968 } 969 } 970 } else { 971 if (addv == INSERT_VALUES) { 972 for (i=0; i<m; i++) { 973 if (indexm[i] < 0) { idx += n; continue;} 974 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); 975 for (j=0; j<n; j++) { 976 if (indexn[j] < 0) { idx++; continue;} 977 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); 978 av[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 979 } 980 } 981 } else { 982 for (i=0; i<m; i++) { 983 if (indexm[i] < 0) { idx += n; continue;} 984 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); 985 for (j=0; j<n; j++) { 986 if (indexn[j] < 0) { idx++; continue;} 987 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); 988 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 989 } 990 } 991 } 992 } 993 /* hack to prevent unneeded copy to the GPU while returning the array */ 994 #if defined(PETSC_HAVE_CUDA) 995 oldf = A->offloadmask; 996 A->offloadmask = PETSC_OFFLOAD_GPU; 997 #endif 998 ierr = MatDenseRestoreArray(A,&av);CHKERRQ(ierr); 999 #if defined(PETSC_HAVE_CUDA) 1000 A->offloadmask = (oldf == PETSC_OFFLOAD_UNALLOCATED ? PETSC_OFFLOAD_UNALLOCATED : PETSC_OFFLOAD_CPU); 1001 #endif 1002 PetscFunctionReturn(0); 1003 } 1004 1005 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 1006 { 1007 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1008 const PetscScalar *vv; 1009 PetscInt i,j; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1014 /* row-oriented output */ 1015 for (i=0; i<m; i++) { 1016 if (indexm[i] < 0) {v += n;continue;} 1017 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); 1018 for (j=0; j<n; j++) { 1019 if (indexn[j] < 0) {v++; continue;} 1020 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); 1021 *v++ = vv[indexn[j]*mat->lda + indexm[i]]; 1022 } 1023 } 1024 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 1028 /* -----------------------------------------------------------------*/ 1029 1030 PetscErrorCode MatView_Dense_Binary(Mat mat,PetscViewer viewer) 1031 { 1032 PetscErrorCode ierr; 1033 PetscBool skipHeader; 1034 PetscViewerFormat format; 1035 PetscInt header[4],M,N,m,lda,i,j,k; 1036 const PetscScalar *v; 1037 PetscScalar *vwork; 1038 1039 PetscFunctionBegin; 1040 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1041 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1042 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1043 if (skipHeader) format = PETSC_VIEWER_NATIVE; 1044 1045 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1046 1047 /* write matrix header */ 1048 header[0] = MAT_FILE_CLASSID; header[1] = M; header[2] = N; 1049 header[3] = (format == PETSC_VIEWER_NATIVE) ? MATRIX_BINARY_FORMAT_DENSE : M*N; 1050 if (!skipHeader) {ierr = PetscViewerBinaryWrite(viewer,header,4,PETSC_INT);CHKERRQ(ierr);} 1051 1052 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1053 if (format != PETSC_VIEWER_NATIVE) { 1054 PetscInt nnz = m*N, *iwork; 1055 /* store row lengths for each row */ 1056 ierr = PetscMalloc1(nnz,&iwork);CHKERRQ(ierr); 1057 for (i=0; i<m; i++) iwork[i] = N; 1058 ierr = PetscViewerBinaryWriteAll(viewer,iwork,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1059 /* store column indices (zero start index) */ 1060 for (k=0, i=0; i<m; i++) 1061 for (j=0; j<N; j++, k++) 1062 iwork[k] = j; 1063 ierr = PetscViewerBinaryWriteAll(viewer,iwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1064 ierr = PetscFree(iwork);CHKERRQ(ierr); 1065 } 1066 /* store matrix values as a dense matrix in row major order */ 1067 ierr = PetscMalloc1(m*N,&vwork);CHKERRQ(ierr); 1068 ierr = MatDenseGetArrayRead(mat,&v);CHKERRQ(ierr); 1069 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1070 for (k=0, i=0; i<m; i++) 1071 for (j=0; j<N; j++, k++) 1072 vwork[k] = v[i+lda*j]; 1073 ierr = MatDenseRestoreArrayRead(mat,&v);CHKERRQ(ierr); 1074 ierr = PetscViewerBinaryWriteAll(viewer,vwork,m*N,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1075 ierr = PetscFree(vwork);CHKERRQ(ierr); 1076 PetscFunctionReturn(0); 1077 } 1078 1079 PetscErrorCode MatLoad_Dense_Binary(Mat mat,PetscViewer viewer) 1080 { 1081 PetscErrorCode ierr; 1082 PetscBool skipHeader; 1083 PetscInt header[4],M,N,m,nz,lda,i,j,k; 1084 PetscInt rows,cols; 1085 PetscScalar *v,*vwork; 1086 1087 PetscFunctionBegin; 1088 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1089 ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);CHKERRQ(ierr); 1090 1091 if (!skipHeader) { 1092 ierr = PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1093 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1094 M = header[1]; N = header[2]; 1095 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 1096 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 1097 nz = header[3]; 1098 if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz); 1099 } else { 1100 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1101 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"); 1102 nz = MATRIX_BINARY_FORMAT_DENSE; 1103 } 1104 1105 /* setup global sizes if not set */ 1106 if (mat->rmap->N < 0) mat->rmap->N = M; 1107 if (mat->cmap->N < 0) mat->cmap->N = N; 1108 ierr = MatSetUp(mat);CHKERRQ(ierr); 1109 /* check if global sizes are correct */ 1110 ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr); 1111 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); 1112 1113 ierr = MatGetSize(mat,NULL,&N);CHKERRQ(ierr); 1114 ierr = MatGetLocalSize(mat,&m,NULL);CHKERRQ(ierr); 1115 ierr = MatDenseGetArray(mat,&v);CHKERRQ(ierr); 1116 ierr = MatDenseGetLDA(mat,&lda);CHKERRQ(ierr); 1117 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense format */ 1118 PetscInt nnz = m*N; 1119 /* read in matrix values */ 1120 ierr = PetscMalloc1(nnz,&vwork);CHKERRQ(ierr); 1121 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1122 /* store values in column major order */ 1123 for (j=0; j<N; j++) 1124 for (i=0; i<m; i++) 1125 v[i+lda*j] = vwork[i*N+j]; 1126 ierr = PetscFree(vwork);CHKERRQ(ierr); 1127 } else { /* matrix in file is sparse format */ 1128 PetscInt nnz = 0, *rlens, *icols; 1129 /* read in row lengths */ 1130 ierr = PetscMalloc1(m,&rlens);CHKERRQ(ierr); 1131 ierr = PetscViewerBinaryReadAll(viewer,rlens,m,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1132 for (i=0; i<m; i++) nnz += rlens[i]; 1133 /* read in column indices and values */ 1134 ierr = PetscMalloc2(nnz,&icols,nnz,&vwork);CHKERRQ(ierr); 1135 ierr = PetscViewerBinaryReadAll(viewer,icols,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_INT);CHKERRQ(ierr); 1136 ierr = PetscViewerBinaryReadAll(viewer,vwork,nnz,PETSC_DETERMINE,PETSC_DETERMINE,PETSC_SCALAR);CHKERRQ(ierr); 1137 /* store values in column major order */ 1138 for (k=0, i=0; i<m; i++) 1139 for (j=0; j<rlens[i]; j++, k++) 1140 v[i+lda*icols[k]] = vwork[k]; 1141 ierr = PetscFree(rlens);CHKERRQ(ierr); 1142 ierr = PetscFree2(icols,vwork);CHKERRQ(ierr); 1143 } 1144 ierr = MatDenseRestoreArray(mat,&v);CHKERRQ(ierr); 1145 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1146 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1147 PetscFunctionReturn(0); 1148 } 1149 1150 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1151 { 1152 PetscBool isbinary, ishdf5; 1153 PetscErrorCode ierr; 1154 1155 PetscFunctionBegin; 1156 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1157 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1158 /* force binary viewer to load .info file if it has not yet done so */ 1159 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1160 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1161 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1162 if (isbinary) { 1163 ierr = MatLoad_Dense_Binary(newMat,viewer);CHKERRQ(ierr); 1164 } else if (ishdf5) { 1165 #if defined(PETSC_HAVE_HDF5) 1166 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1167 #else 1168 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1169 #endif 1170 } else { 1171 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); 1172 } 1173 PetscFunctionReturn(0); 1174 } 1175 1176 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1177 { 1178 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1179 PetscErrorCode ierr; 1180 PetscInt i,j; 1181 const char *name; 1182 PetscScalar *v,*av; 1183 PetscViewerFormat format; 1184 #if defined(PETSC_USE_COMPLEX) 1185 PetscBool allreal = PETSC_TRUE; 1186 #endif 1187 1188 PetscFunctionBegin; 1189 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1190 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1191 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1192 PetscFunctionReturn(0); /* do nothing for now */ 1193 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1194 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1195 for (i=0; i<A->rmap->n; i++) { 1196 v = av + i; 1197 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1198 for (j=0; j<A->cmap->n; j++) { 1199 #if defined(PETSC_USE_COMPLEX) 1200 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1201 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1202 } else if (PetscRealPart(*v)) { 1203 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1204 } 1205 #else 1206 if (*v) { 1207 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1208 } 1209 #endif 1210 v += a->lda; 1211 } 1212 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1213 } 1214 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1215 } else { 1216 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1217 #if defined(PETSC_USE_COMPLEX) 1218 /* determine if matrix has all real values */ 1219 v = av; 1220 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1221 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1222 } 1223 #endif 1224 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1225 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1226 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1228 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1229 } 1230 1231 for (i=0; i<A->rmap->n; i++) { 1232 v = av + i; 1233 for (j=0; j<A->cmap->n; j++) { 1234 #if defined(PETSC_USE_COMPLEX) 1235 if (allreal) { 1236 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1237 } else { 1238 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1239 } 1240 #else 1241 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1242 #endif 1243 v += a->lda; 1244 } 1245 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1246 } 1247 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1248 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1249 } 1250 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1251 } 1252 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&av);CHKERRQ(ierr); 1253 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1254 PetscFunctionReturn(0); 1255 } 1256 1257 #include <petscdraw.h> 1258 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1259 { 1260 Mat A = (Mat) Aa; 1261 PetscErrorCode ierr; 1262 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1263 int color = PETSC_DRAW_WHITE; 1264 const PetscScalar *v; 1265 PetscViewer viewer; 1266 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1267 PetscViewerFormat format; 1268 1269 PetscFunctionBegin; 1270 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1271 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1272 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1273 1274 /* Loop over matrix elements drawing boxes */ 1275 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 1276 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1277 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1278 /* Blue for negative and Red for positive */ 1279 for (j = 0; j < n; j++) { 1280 x_l = j; x_r = x_l + 1.0; 1281 for (i = 0; i < m; i++) { 1282 y_l = m - i - 1.0; 1283 y_r = y_l + 1.0; 1284 if (PetscRealPart(v[j*m+i]) > 0.) color = PETSC_DRAW_RED; 1285 else if (PetscRealPart(v[j*m+i]) < 0.) color = PETSC_DRAW_BLUE; 1286 else continue; 1287 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1288 } 1289 } 1290 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1291 } else { 1292 /* use contour shading to indicate magnitude of values */ 1293 /* first determine max of all nonzero values */ 1294 PetscReal minv = 0.0, maxv = 0.0; 1295 PetscDraw popup; 1296 1297 for (i=0; i < m*n; i++) { 1298 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1299 } 1300 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1301 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1302 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1303 1304 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1305 for (j=0; j<n; j++) { 1306 x_l = j; 1307 x_r = x_l + 1.0; 1308 for (i=0; i<m; i++) { 1309 y_l = m - i - 1.0; 1310 y_r = y_l + 1.0; 1311 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1312 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1313 } 1314 } 1315 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1316 } 1317 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 1318 PetscFunctionReturn(0); 1319 } 1320 1321 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1322 { 1323 PetscDraw draw; 1324 PetscBool isnull; 1325 PetscReal xr,yr,xl,yl,h,w; 1326 PetscErrorCode ierr; 1327 1328 PetscFunctionBegin; 1329 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1330 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1331 if (isnull) PetscFunctionReturn(0); 1332 1333 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1334 xr += w; yr += h; xl = -w; yl = -h; 1335 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1336 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1337 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1338 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1339 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1344 { 1345 PetscErrorCode ierr; 1346 PetscBool iascii,isbinary,isdraw; 1347 1348 PetscFunctionBegin; 1349 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1350 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1351 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1352 1353 if (iascii) { 1354 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1355 } else if (isbinary) { 1356 ierr = MatView_Dense_Binary(A,viewer);CHKERRQ(ierr); 1357 } else if (isdraw) { 1358 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1359 } 1360 PetscFunctionReturn(0); 1361 } 1362 1363 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array) 1364 { 1365 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1366 1367 PetscFunctionBegin; 1368 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1369 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1370 if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first"); 1371 a->unplacedarray = a->v; 1372 a->unplaced_user_alloc = a->user_alloc; 1373 a->v = (PetscScalar*) array; 1374 a->user_alloc = PETSC_TRUE; 1375 #if defined(PETSC_HAVE_CUDA) 1376 A->offloadmask = PETSC_OFFLOAD_CPU; 1377 #endif 1378 PetscFunctionReturn(0); 1379 } 1380 1381 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1382 { 1383 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1384 1385 PetscFunctionBegin; 1386 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1387 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1388 a->v = a->unplacedarray; 1389 a->user_alloc = a->unplaced_user_alloc; 1390 a->unplacedarray = NULL; 1391 #if defined(PETSC_HAVE_CUDA) 1392 A->offloadmask = PETSC_OFFLOAD_CPU; 1393 #endif 1394 PetscFunctionReturn(0); 1395 } 1396 1397 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array) 1398 { 1399 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1400 PetscErrorCode ierr; 1401 1402 PetscFunctionBegin; 1403 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1404 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1405 if (!a->user_alloc) { ierr = PetscFree(a->v);CHKERRQ(ierr); } 1406 a->v = (PetscScalar*) array; 1407 a->user_alloc = PETSC_FALSE; 1408 #if defined(PETSC_HAVE_CUDA) 1409 A->offloadmask = PETSC_OFFLOAD_CPU; 1410 #endif 1411 PetscFunctionReturn(0); 1412 } 1413 1414 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1415 { 1416 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1417 PetscErrorCode ierr; 1418 1419 PetscFunctionBegin; 1420 #if defined(PETSC_USE_LOG) 1421 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1422 #endif 1423 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1424 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1425 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1426 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1427 if (!l->unplaced_user_alloc) {ierr = PetscFree(l->unplacedarray);CHKERRQ(ierr);} 1428 if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1429 if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1430 ierr = VecDestroy(&l->cvec);CHKERRQ(ierr); 1431 ierr = MatDestroy(&l->cmat);CHKERRQ(ierr); 1432 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1433 1434 ierr = PetscObjectChangeTypeName((PetscObject)mat,NULL);CHKERRQ(ierr); 1435 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1436 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);CHKERRQ(ierr); 1437 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1438 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1439 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1440 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1441 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);CHKERRQ(ierr); 1442 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr); 1443 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr); 1444 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);CHKERRQ(ierr); 1445 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);CHKERRQ(ierr); 1446 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1447 #if defined(PETSC_HAVE_ELEMENTAL) 1448 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1449 #endif 1450 #if defined(PETSC_HAVE_SCALAPACK) 1451 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);CHKERRQ(ierr); 1452 #endif 1453 #if defined(PETSC_HAVE_CUDA) 1454 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1455 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr); 1456 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr); 1457 #endif 1458 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1459 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1460 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr); 1461 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1463 1464 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1466 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);CHKERRQ(ierr); 1467 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);CHKERRQ(ierr); 1468 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);CHKERRQ(ierr); 1469 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);CHKERRQ(ierr); 1470 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);CHKERRQ(ierr); 1471 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);CHKERRQ(ierr); 1472 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);CHKERRQ(ierr); 1473 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);CHKERRQ(ierr); 1474 PetscFunctionReturn(0); 1475 } 1476 1477 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1478 { 1479 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1480 PetscErrorCode ierr; 1481 PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n; 1482 PetscScalar *v,tmp; 1483 1484 PetscFunctionBegin; 1485 if (reuse == MAT_INPLACE_MATRIX) { 1486 if (m == n) { /* in place transpose */ 1487 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1488 for (j=0; j<m; j++) { 1489 for (k=0; k<j; k++) { 1490 tmp = v[j + k*M]; 1491 v[j + k*M] = v[k + j*M]; 1492 v[k + j*M] = tmp; 1493 } 1494 } 1495 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1496 } else { /* reuse memory, temporary allocates new memory */ 1497 PetscScalar *v2; 1498 PetscLayout tmplayout; 1499 1500 ierr = PetscMalloc1((size_t)m*n,&v2);CHKERRQ(ierr); 1501 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1502 for (j=0; j<n; j++) { 1503 for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M]; 1504 } 1505 ierr = PetscArraycpy(v,v2,(size_t)m*n);CHKERRQ(ierr); 1506 ierr = PetscFree(v2);CHKERRQ(ierr); 1507 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1508 /* cleanup size dependent quantities */ 1509 ierr = VecDestroy(&mat->cvec);CHKERRQ(ierr); 1510 ierr = MatDestroy(&mat->cmat);CHKERRQ(ierr); 1511 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 1512 ierr = PetscFree(mat->fwork);CHKERRQ(ierr); 1513 ierr = MatDestroy(&mat->ptapwork);CHKERRQ(ierr); 1514 /* swap row/col layouts */ 1515 mat->lda = n; 1516 tmplayout = A->rmap; 1517 A->rmap = A->cmap; 1518 A->cmap = tmplayout; 1519 } 1520 } else { /* out-of-place transpose */ 1521 Mat tmat; 1522 Mat_SeqDense *tmatd; 1523 PetscScalar *v2; 1524 PetscInt M2; 1525 1526 if (reuse == MAT_INITIAL_MATRIX) { 1527 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1528 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1529 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1530 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1531 } else tmat = *matout; 1532 1533 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1534 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1535 tmatd = (Mat_SeqDense*)tmat->data; 1536 M2 = tmatd->lda; 1537 for (j=0; j<n; j++) { 1538 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1539 } 1540 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1541 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1542 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1543 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 *matout = tmat; 1545 } 1546 PetscFunctionReturn(0); 1547 } 1548 1549 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1550 { 1551 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1552 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1553 PetscInt i; 1554 const PetscScalar *v1,*v2; 1555 PetscErrorCode ierr; 1556 1557 PetscFunctionBegin; 1558 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1559 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1560 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1561 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1562 for (i=0; i<A1->cmap->n; i++) { 1563 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1564 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1565 v1 += mat1->lda; 1566 v2 += mat2->lda; 1567 } 1568 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1569 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1570 *flg = PETSC_TRUE; 1571 PetscFunctionReturn(0); 1572 } 1573 1574 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1575 { 1576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1577 PetscInt i,n,len; 1578 PetscScalar *x; 1579 const PetscScalar *vv; 1580 PetscErrorCode ierr; 1581 1582 PetscFunctionBegin; 1583 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1584 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1585 len = PetscMin(A->rmap->n,A->cmap->n); 1586 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1587 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1588 for (i=0; i<len; i++) { 1589 x[i] = vv[i*mat->lda + i]; 1590 } 1591 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1592 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1593 PetscFunctionReturn(0); 1594 } 1595 1596 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1597 { 1598 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1599 const PetscScalar *l,*r; 1600 PetscScalar x,*v,*vv; 1601 PetscErrorCode ierr; 1602 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1603 1604 PetscFunctionBegin; 1605 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1606 if (ll) { 1607 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1608 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1609 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1610 for (i=0; i<m; i++) { 1611 x = l[i]; 1612 v = vv + i; 1613 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1614 } 1615 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1616 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1617 } 1618 if (rr) { 1619 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1620 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1621 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1622 for (i=0; i<n; i++) { 1623 x = r[i]; 1624 v = vv + i*mat->lda; 1625 for (j=0; j<m; j++) (*v++) *= x; 1626 } 1627 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1628 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1629 } 1630 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1631 PetscFunctionReturn(0); 1632 } 1633 1634 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1635 { 1636 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1637 PetscScalar *v,*vv; 1638 PetscReal sum = 0.0; 1639 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1640 PetscErrorCode ierr; 1641 1642 PetscFunctionBegin; 1643 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1644 v = vv; 1645 if (type == NORM_FROBENIUS) { 1646 if (lda>m) { 1647 for (j=0; j<A->cmap->n; j++) { 1648 v = vv+j*lda; 1649 for (i=0; i<m; i++) { 1650 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1651 } 1652 } 1653 } else { 1654 #if defined(PETSC_USE_REAL___FP16) 1655 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1656 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one)); 1657 } 1658 #else 1659 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1660 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1661 } 1662 } 1663 *nrm = PetscSqrtReal(sum); 1664 #endif 1665 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1666 } else if (type == NORM_1) { 1667 *nrm = 0.0; 1668 for (j=0; j<A->cmap->n; j++) { 1669 v = vv + j*mat->lda; 1670 sum = 0.0; 1671 for (i=0; i<A->rmap->n; i++) { 1672 sum += PetscAbsScalar(*v); v++; 1673 } 1674 if (sum > *nrm) *nrm = sum; 1675 } 1676 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1677 } else if (type == NORM_INFINITY) { 1678 *nrm = 0.0; 1679 for (j=0; j<A->rmap->n; j++) { 1680 v = vv + j; 1681 sum = 0.0; 1682 for (i=0; i<A->cmap->n; i++) { 1683 sum += PetscAbsScalar(*v); v += mat->lda; 1684 } 1685 if (sum > *nrm) *nrm = sum; 1686 } 1687 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1688 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1689 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1690 PetscFunctionReturn(0); 1691 } 1692 1693 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1694 { 1695 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1696 PetscErrorCode ierr; 1697 1698 PetscFunctionBegin; 1699 switch (op) { 1700 case MAT_ROW_ORIENTED: 1701 aij->roworiented = flg; 1702 break; 1703 case MAT_NEW_NONZERO_LOCATIONS: 1704 case MAT_NEW_NONZERO_LOCATION_ERR: 1705 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1706 case MAT_FORCE_DIAGONAL_ENTRIES: 1707 case MAT_KEEP_NONZERO_PATTERN: 1708 case MAT_IGNORE_OFF_PROC_ENTRIES: 1709 case MAT_USE_HASH_TABLE: 1710 case MAT_IGNORE_ZERO_ENTRIES: 1711 case MAT_IGNORE_LOWER_TRIANGULAR: 1712 case MAT_SORTED_FULL: 1713 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1714 break; 1715 case MAT_SPD: 1716 case MAT_SYMMETRIC: 1717 case MAT_STRUCTURALLY_SYMMETRIC: 1718 case MAT_HERMITIAN: 1719 case MAT_SYMMETRY_ETERNAL: 1720 /* These options are handled directly by MatSetOption() */ 1721 break; 1722 default: 1723 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1724 } 1725 PetscFunctionReturn(0); 1726 } 1727 1728 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1729 { 1730 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1731 PetscErrorCode ierr; 1732 PetscInt lda=l->lda,m=A->rmap->n,j; 1733 PetscScalar *v; 1734 1735 PetscFunctionBegin; 1736 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1737 if (lda>m) { 1738 for (j=0; j<A->cmap->n; j++) { 1739 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1740 } 1741 } else { 1742 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1743 } 1744 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1745 PetscFunctionReturn(0); 1746 } 1747 1748 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1749 { 1750 PetscErrorCode ierr; 1751 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1752 PetscInt m = l->lda, n = A->cmap->n, i,j; 1753 PetscScalar *slot,*bb,*v; 1754 const PetscScalar *xx; 1755 1756 PetscFunctionBegin; 1757 if (PetscDefined(USE_DEBUG)) { 1758 for (i=0; i<N; i++) { 1759 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1760 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); 1761 } 1762 } 1763 if (!N) PetscFunctionReturn(0); 1764 1765 /* fix right hand side if needed */ 1766 if (x && b) { 1767 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1768 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1769 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1770 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1771 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1772 } 1773 1774 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1775 for (i=0; i<N; i++) { 1776 slot = v + rows[i]; 1777 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1778 } 1779 if (diag != 0.0) { 1780 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1781 for (i=0; i<N; i++) { 1782 slot = v + (m+1)*rows[i]; 1783 *slot = diag; 1784 } 1785 } 1786 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1787 PetscFunctionReturn(0); 1788 } 1789 1790 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1791 { 1792 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1793 1794 PetscFunctionBegin; 1795 *lda = mat->lda; 1796 PetscFunctionReturn(0); 1797 } 1798 1799 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 1800 { 1801 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1802 1803 PetscFunctionBegin; 1804 if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1805 *array = mat->v; 1806 PetscFunctionReturn(0); 1807 } 1808 1809 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 1810 { 1811 PetscFunctionBegin; 1812 *array = NULL; 1813 PetscFunctionReturn(0); 1814 } 1815 1816 /*@C 1817 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1818 1819 Not collective 1820 1821 Input Parameter: 1822 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1823 1824 Output Parameter: 1825 . lda - the leading dimension 1826 1827 Level: intermediate 1828 1829 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA() 1830 @*/ 1831 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1832 { 1833 PetscErrorCode ierr; 1834 1835 PetscFunctionBegin; 1836 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1837 PetscValidPointer(lda,2); 1838 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1839 PetscFunctionReturn(0); 1840 } 1841 1842 /*@C 1843 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 1844 1845 Not collective 1846 1847 Input Parameter: 1848 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1849 - lda - the leading dimension 1850 1851 Level: intermediate 1852 1853 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA() 1854 @*/ 1855 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 1856 { 1857 PetscErrorCode ierr; 1858 1859 PetscFunctionBegin; 1860 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1861 ierr = PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));CHKERRQ(ierr); 1862 PetscFunctionReturn(0); 1863 } 1864 1865 /*@C 1866 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 1867 1868 Logically Collective on Mat 1869 1870 Input Parameter: 1871 . mat - a dense matrix 1872 1873 Output Parameter: 1874 . array - pointer to the data 1875 1876 Level: intermediate 1877 1878 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1879 @*/ 1880 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1881 { 1882 PetscErrorCode ierr; 1883 1884 PetscFunctionBegin; 1885 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1886 PetscValidPointer(array,2); 1887 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 /*@C 1892 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1893 1894 Logically Collective on Mat 1895 1896 Input Parameters: 1897 + mat - a dense matrix 1898 - array - pointer to the data 1899 1900 Level: intermediate 1901 1902 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1903 @*/ 1904 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1905 { 1906 PetscErrorCode ierr; 1907 1908 PetscFunctionBegin; 1909 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1910 PetscValidPointer(array,2); 1911 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1912 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1913 #if defined(PETSC_HAVE_CUDA) 1914 A->offloadmask = PETSC_OFFLOAD_CPU; 1915 #endif 1916 PetscFunctionReturn(0); 1917 } 1918 1919 /*@C 1920 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 1921 1922 Not Collective 1923 1924 Input Parameter: 1925 . mat - a dense matrix 1926 1927 Output Parameter: 1928 . array - pointer to the data 1929 1930 Level: intermediate 1931 1932 .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1933 @*/ 1934 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1935 { 1936 PetscErrorCode ierr; 1937 1938 PetscFunctionBegin; 1939 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1940 PetscValidPointer(array,2); 1941 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1942 PetscFunctionReturn(0); 1943 } 1944 1945 /*@C 1946 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 1947 1948 Not Collective 1949 1950 Input Parameters: 1951 + mat - a dense matrix 1952 - array - pointer to the data 1953 1954 Level: intermediate 1955 1956 .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1957 @*/ 1958 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1959 { 1960 PetscErrorCode ierr; 1961 1962 PetscFunctionBegin; 1963 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1964 PetscValidPointer(array,2); 1965 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1966 PetscFunctionReturn(0); 1967 } 1968 1969 /*@C 1970 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 1971 1972 Not Collective 1973 1974 Input Parameter: 1975 . mat - a dense matrix 1976 1977 Output Parameter: 1978 . array - pointer to the data 1979 1980 Level: intermediate 1981 1982 .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1983 @*/ 1984 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 1985 { 1986 PetscErrorCode ierr; 1987 1988 PetscFunctionBegin; 1989 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1990 PetscValidPointer(array,2); 1991 ierr = PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1992 PetscFunctionReturn(0); 1993 } 1994 1995 /*@C 1996 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 1997 1998 Not Collective 1999 2000 Input Parameters: 2001 + mat - a dense matrix 2002 - array - pointer to the data 2003 2004 Level: intermediate 2005 2006 .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 2007 @*/ 2008 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 2009 { 2010 PetscErrorCode ierr; 2011 2012 PetscFunctionBegin; 2013 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2014 PetscValidPointer(array,2); 2015 ierr = PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 2016 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 2017 #if defined(PETSC_HAVE_CUDA) 2018 A->offloadmask = PETSC_OFFLOAD_CPU; 2019 #endif 2020 PetscFunctionReturn(0); 2021 } 2022 2023 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 2024 { 2025 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2026 PetscErrorCode ierr; 2027 PetscInt i,j,nrows,ncols,blda; 2028 const PetscInt *irow,*icol; 2029 PetscScalar *av,*bv,*v = mat->v; 2030 Mat newmat; 2031 2032 PetscFunctionBegin; 2033 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2034 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2035 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2036 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2037 2038 /* Check submatrixcall */ 2039 if (scall == MAT_REUSE_MATRIX) { 2040 PetscInt n_cols,n_rows; 2041 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2042 if (n_rows != nrows || n_cols != ncols) { 2043 /* resize the result matrix to match number of requested rows/columns */ 2044 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2045 } 2046 newmat = *B; 2047 } else { 2048 /* Create and fill new matrix */ 2049 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 2050 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2051 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 2052 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 2053 } 2054 2055 /* Now extract the data pointers and do the copy,column at a time */ 2056 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 2057 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 2058 for (i=0; i<ncols; i++) { 2059 av = v + mat->lda*icol[i]; 2060 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2061 bv += blda; 2062 } 2063 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 2064 2065 /* Assemble the matrices so that the correct flags are set */ 2066 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2067 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2068 2069 /* Free work space */ 2070 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2071 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2072 *B = newmat; 2073 PetscFunctionReturn(0); 2074 } 2075 2076 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2077 { 2078 PetscErrorCode ierr; 2079 PetscInt i; 2080 2081 PetscFunctionBegin; 2082 if (scall == MAT_INITIAL_MATRIX) { 2083 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2084 } 2085 2086 for (i=0; i<n; i++) { 2087 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2088 } 2089 PetscFunctionReturn(0); 2090 } 2091 2092 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2093 { 2094 PetscFunctionBegin; 2095 PetscFunctionReturn(0); 2096 } 2097 2098 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2099 { 2100 PetscFunctionBegin; 2101 PetscFunctionReturn(0); 2102 } 2103 2104 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2105 { 2106 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2107 PetscErrorCode ierr; 2108 const PetscScalar *va; 2109 PetscScalar *vb; 2110 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2111 2112 PetscFunctionBegin; 2113 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2114 if (A->ops->copy != B->ops->copy) { 2115 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2116 PetscFunctionReturn(0); 2117 } 2118 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2119 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2120 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2121 if (lda1>m || lda2>m) { 2122 for (j=0; j<n; j++) { 2123 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2124 } 2125 } else { 2126 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2127 } 2128 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2129 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2130 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2131 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2136 { 2137 PetscErrorCode ierr; 2138 2139 PetscFunctionBegin; 2140 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2141 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2142 if (!A->preallocated) { 2143 ierr = MatSeqDenseSetPreallocation(A,NULL);CHKERRQ(ierr); 2144 } 2145 PetscFunctionReturn(0); 2146 } 2147 2148 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2149 { 2150 PetscInt i,nz = A->rmap->n*A->cmap->n; 2151 PetscScalar *aa; 2152 PetscErrorCode ierr; 2153 2154 PetscFunctionBegin; 2155 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2156 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2157 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2158 PetscFunctionReturn(0); 2159 } 2160 2161 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2162 { 2163 PetscInt i,nz = A->rmap->n*A->cmap->n; 2164 PetscScalar *aa; 2165 PetscErrorCode ierr; 2166 2167 PetscFunctionBegin; 2168 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2169 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2170 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2171 PetscFunctionReturn(0); 2172 } 2173 2174 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2175 { 2176 PetscInt i,nz = A->rmap->n*A->cmap->n; 2177 PetscScalar *aa; 2178 PetscErrorCode ierr; 2179 2180 PetscFunctionBegin; 2181 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2182 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2183 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2184 PetscFunctionReturn(0); 2185 } 2186 2187 /* ----------------------------------------------------------------*/ 2188 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2189 { 2190 PetscErrorCode ierr; 2191 PetscInt m=A->rmap->n,n=B->cmap->n; 2192 PetscBool cisdense; 2193 2194 PetscFunctionBegin; 2195 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2196 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2197 if (!cisdense) { 2198 PetscBool flg; 2199 2200 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2201 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2202 } 2203 ierr = MatSetUp(C);CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2208 { 2209 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2210 PetscBLASInt m,n,k; 2211 const PetscScalar *av,*bv; 2212 PetscScalar *cv; 2213 PetscScalar _DOne=1.0,_DZero=0.0; 2214 PetscErrorCode ierr; 2215 2216 PetscFunctionBegin; 2217 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2218 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2219 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2220 if (!m || !n || !k) PetscFunctionReturn(0); 2221 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2222 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2223 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2224 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2225 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2226 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2227 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2228 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2229 PetscFunctionReturn(0); 2230 } 2231 2232 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2233 { 2234 PetscErrorCode ierr; 2235 PetscInt m=A->rmap->n,n=B->rmap->n; 2236 PetscBool cisdense; 2237 2238 PetscFunctionBegin; 2239 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2240 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2241 if (!cisdense) { 2242 PetscBool flg; 2243 2244 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2245 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2246 } 2247 ierr = MatSetUp(C);CHKERRQ(ierr); 2248 PetscFunctionReturn(0); 2249 } 2250 2251 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2252 { 2253 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2254 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2255 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2256 const PetscScalar *av,*bv; 2257 PetscScalar *cv; 2258 PetscBLASInt m,n,k; 2259 PetscScalar _DOne=1.0,_DZero=0.0; 2260 PetscErrorCode ierr; 2261 2262 PetscFunctionBegin; 2263 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2264 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2265 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2266 if (!m || !n || !k) PetscFunctionReturn(0); 2267 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2268 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2269 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2270 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2271 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2272 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2273 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2274 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2275 PetscFunctionReturn(0); 2276 } 2277 2278 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2279 { 2280 PetscErrorCode ierr; 2281 PetscInt m=A->cmap->n,n=B->cmap->n; 2282 PetscBool cisdense; 2283 2284 PetscFunctionBegin; 2285 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2286 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2287 if (!cisdense) { 2288 PetscBool flg; 2289 2290 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2291 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2292 } 2293 ierr = MatSetUp(C);CHKERRQ(ierr); 2294 PetscFunctionReturn(0); 2295 } 2296 2297 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2298 { 2299 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2300 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2301 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2302 const PetscScalar *av,*bv; 2303 PetscScalar *cv; 2304 PetscBLASInt m,n,k; 2305 PetscScalar _DOne=1.0,_DZero=0.0; 2306 PetscErrorCode ierr; 2307 2308 PetscFunctionBegin; 2309 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2310 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2311 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2312 if (!m || !n || !k) PetscFunctionReturn(0); 2313 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2314 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2315 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2316 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2317 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2318 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2319 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2320 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2321 PetscFunctionReturn(0); 2322 } 2323 2324 /* ----------------------------------------------- */ 2325 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2326 { 2327 PetscFunctionBegin; 2328 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2329 C->ops->productsymbolic = MatProductSymbolic_AB; 2330 PetscFunctionReturn(0); 2331 } 2332 2333 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2334 { 2335 PetscFunctionBegin; 2336 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2337 C->ops->productsymbolic = MatProductSymbolic_AtB; 2338 PetscFunctionReturn(0); 2339 } 2340 2341 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2342 { 2343 PetscFunctionBegin; 2344 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2345 C->ops->productsymbolic = MatProductSymbolic_ABt; 2346 PetscFunctionReturn(0); 2347 } 2348 2349 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2350 { 2351 PetscErrorCode ierr; 2352 Mat_Product *product = C->product; 2353 2354 PetscFunctionBegin; 2355 switch (product->type) { 2356 case MATPRODUCT_AB: 2357 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2358 break; 2359 case MATPRODUCT_AtB: 2360 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2361 break; 2362 case MATPRODUCT_ABt: 2363 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2364 break; 2365 default: 2366 break; 2367 } 2368 PetscFunctionReturn(0); 2369 } 2370 /* ----------------------------------------------- */ 2371 2372 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2373 { 2374 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2375 PetscErrorCode ierr; 2376 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2377 PetscScalar *x; 2378 const PetscScalar *aa; 2379 2380 PetscFunctionBegin; 2381 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2382 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2383 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2384 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2385 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2386 for (i=0; i<m; i++) { 2387 x[i] = aa[i]; if (idx) idx[i] = 0; 2388 for (j=1; j<n; j++) { 2389 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2390 } 2391 } 2392 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2393 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2394 PetscFunctionReturn(0); 2395 } 2396 2397 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2398 { 2399 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2400 PetscErrorCode ierr; 2401 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2402 PetscScalar *x; 2403 PetscReal atmp; 2404 const PetscScalar *aa; 2405 2406 PetscFunctionBegin; 2407 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2408 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2409 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2410 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2411 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2412 for (i=0; i<m; i++) { 2413 x[i] = PetscAbsScalar(aa[i]); 2414 for (j=1; j<n; j++) { 2415 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2416 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2417 } 2418 } 2419 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2420 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2421 PetscFunctionReturn(0); 2422 } 2423 2424 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2425 { 2426 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2427 PetscErrorCode ierr; 2428 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2429 PetscScalar *x; 2430 const PetscScalar *aa; 2431 2432 PetscFunctionBegin; 2433 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2434 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2435 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2436 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2437 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2438 for (i=0; i<m; i++) { 2439 x[i] = aa[i]; if (idx) idx[i] = 0; 2440 for (j=1; j<n; j++) { 2441 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2442 } 2443 } 2444 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2445 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2446 PetscFunctionReturn(0); 2447 } 2448 2449 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2450 { 2451 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2452 PetscErrorCode ierr; 2453 PetscScalar *x; 2454 const PetscScalar *aa; 2455 2456 PetscFunctionBegin; 2457 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2458 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2459 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2460 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2461 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2462 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2463 PetscFunctionReturn(0); 2464 } 2465 2466 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2467 { 2468 PetscErrorCode ierr; 2469 PetscInt i,j,m,n; 2470 const PetscScalar *a; 2471 2472 PetscFunctionBegin; 2473 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2474 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2475 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2476 if (type == NORM_2) { 2477 for (i=0; i<n; i++) { 2478 for (j=0; j<m; j++) { 2479 norms[i] += PetscAbsScalar(a[j]*a[j]); 2480 } 2481 a += m; 2482 } 2483 } else if (type == NORM_1) { 2484 for (i=0; i<n; i++) { 2485 for (j=0; j<m; j++) { 2486 norms[i] += PetscAbsScalar(a[j]); 2487 } 2488 a += m; 2489 } 2490 } else if (type == NORM_INFINITY) { 2491 for (i=0; i<n; i++) { 2492 for (j=0; j<m; j++) { 2493 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2494 } 2495 a += m; 2496 } 2497 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2498 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2499 if (type == NORM_2) { 2500 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2501 } 2502 PetscFunctionReturn(0); 2503 } 2504 2505 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2506 { 2507 PetscErrorCode ierr; 2508 PetscScalar *a; 2509 PetscInt lda,m,n,i,j; 2510 2511 PetscFunctionBegin; 2512 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2513 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2514 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2515 for (j=0; j<n; j++) { 2516 for (i=0; i<m; i++) { 2517 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2518 } 2519 } 2520 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2521 PetscFunctionReturn(0); 2522 } 2523 2524 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2525 { 2526 PetscFunctionBegin; 2527 *missing = PETSC_FALSE; 2528 PetscFunctionReturn(0); 2529 } 2530 2531 /* vals is not const */ 2532 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2533 { 2534 PetscErrorCode ierr; 2535 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2536 PetscScalar *v; 2537 2538 PetscFunctionBegin; 2539 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2540 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2541 *vals = v+col*a->lda; 2542 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2543 PetscFunctionReturn(0); 2544 } 2545 2546 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2547 { 2548 PetscFunctionBegin; 2549 *vals = NULL; /* user cannot accidently use the array later */ 2550 PetscFunctionReturn(0); 2551 } 2552 2553 /* -------------------------------------------------------------------*/ 2554 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2555 MatGetRow_SeqDense, 2556 MatRestoreRow_SeqDense, 2557 MatMult_SeqDense, 2558 /* 4*/ MatMultAdd_SeqDense, 2559 MatMultTranspose_SeqDense, 2560 MatMultTransposeAdd_SeqDense, 2561 NULL, 2562 NULL, 2563 NULL, 2564 /* 10*/ NULL, 2565 MatLUFactor_SeqDense, 2566 MatCholeskyFactor_SeqDense, 2567 MatSOR_SeqDense, 2568 MatTranspose_SeqDense, 2569 /* 15*/ MatGetInfo_SeqDense, 2570 MatEqual_SeqDense, 2571 MatGetDiagonal_SeqDense, 2572 MatDiagonalScale_SeqDense, 2573 MatNorm_SeqDense, 2574 /* 20*/ MatAssemblyBegin_SeqDense, 2575 MatAssemblyEnd_SeqDense, 2576 MatSetOption_SeqDense, 2577 MatZeroEntries_SeqDense, 2578 /* 24*/ MatZeroRows_SeqDense, 2579 NULL, 2580 NULL, 2581 NULL, 2582 NULL, 2583 /* 29*/ MatSetUp_SeqDense, 2584 NULL, 2585 NULL, 2586 NULL, 2587 NULL, 2588 /* 34*/ MatDuplicate_SeqDense, 2589 NULL, 2590 NULL, 2591 NULL, 2592 NULL, 2593 /* 39*/ MatAXPY_SeqDense, 2594 MatCreateSubMatrices_SeqDense, 2595 NULL, 2596 MatGetValues_SeqDense, 2597 MatCopy_SeqDense, 2598 /* 44*/ MatGetRowMax_SeqDense, 2599 MatScale_SeqDense, 2600 MatShift_Basic, 2601 NULL, 2602 MatZeroRowsColumns_SeqDense, 2603 /* 49*/ MatSetRandom_SeqDense, 2604 NULL, 2605 NULL, 2606 NULL, 2607 NULL, 2608 /* 54*/ NULL, 2609 NULL, 2610 NULL, 2611 NULL, 2612 NULL, 2613 /* 59*/ NULL, 2614 MatDestroy_SeqDense, 2615 MatView_SeqDense, 2616 NULL, 2617 NULL, 2618 /* 64*/ NULL, 2619 NULL, 2620 NULL, 2621 NULL, 2622 NULL, 2623 /* 69*/ MatGetRowMaxAbs_SeqDense, 2624 NULL, 2625 NULL, 2626 NULL, 2627 NULL, 2628 /* 74*/ NULL, 2629 NULL, 2630 NULL, 2631 NULL, 2632 NULL, 2633 /* 79*/ NULL, 2634 NULL, 2635 NULL, 2636 NULL, 2637 /* 83*/ MatLoad_SeqDense, 2638 MatIsSymmetric_SeqDense, 2639 MatIsHermitian_SeqDense, 2640 NULL, 2641 NULL, 2642 NULL, 2643 /* 89*/ NULL, 2644 NULL, 2645 MatMatMultNumeric_SeqDense_SeqDense, 2646 NULL, 2647 NULL, 2648 /* 94*/ NULL, 2649 NULL, 2650 NULL, 2651 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2652 NULL, 2653 /* 99*/ MatProductSetFromOptions_SeqDense, 2654 NULL, 2655 NULL, 2656 MatConjugate_SeqDense, 2657 NULL, 2658 /*104*/ NULL, 2659 MatRealPart_SeqDense, 2660 MatImaginaryPart_SeqDense, 2661 NULL, 2662 NULL, 2663 /*109*/ NULL, 2664 NULL, 2665 MatGetRowMin_SeqDense, 2666 MatGetColumnVector_SeqDense, 2667 MatMissingDiagonal_SeqDense, 2668 /*114*/ NULL, 2669 NULL, 2670 NULL, 2671 NULL, 2672 NULL, 2673 /*119*/ NULL, 2674 NULL, 2675 NULL, 2676 NULL, 2677 NULL, 2678 /*124*/ NULL, 2679 MatGetColumnNorms_SeqDense, 2680 NULL, 2681 NULL, 2682 NULL, 2683 /*129*/ NULL, 2684 NULL, 2685 NULL, 2686 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2687 NULL, 2688 /*134*/ NULL, 2689 NULL, 2690 NULL, 2691 NULL, 2692 NULL, 2693 /*139*/ NULL, 2694 NULL, 2695 NULL, 2696 NULL, 2697 NULL, 2698 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2699 /*145*/ NULL, 2700 NULL, 2701 NULL 2702 }; 2703 2704 /*@C 2705 MatCreateSeqDense - Creates a sequential dense matrix that 2706 is stored in column major order (the usual Fortran 77 manner). Many 2707 of the matrix operations use the BLAS and LAPACK routines. 2708 2709 Collective 2710 2711 Input Parameters: 2712 + comm - MPI communicator, set to PETSC_COMM_SELF 2713 . m - number of rows 2714 . n - number of columns 2715 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2716 to control all matrix memory allocation. 2717 2718 Output Parameter: 2719 . A - the matrix 2720 2721 Notes: 2722 The data input variable is intended primarily for Fortran programmers 2723 who wish to allocate their own matrix memory space. Most users should 2724 set data=NULL. 2725 2726 Level: intermediate 2727 2728 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2729 @*/ 2730 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2731 { 2732 PetscErrorCode ierr; 2733 2734 PetscFunctionBegin; 2735 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2736 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2737 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2738 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2739 PetscFunctionReturn(0); 2740 } 2741 2742 /*@C 2743 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2744 2745 Collective 2746 2747 Input Parameters: 2748 + B - the matrix 2749 - data - the array (or NULL) 2750 2751 Notes: 2752 The data input variable is intended primarily for Fortran programmers 2753 who wish to allocate their own matrix memory space. Most users should 2754 need not call this routine. 2755 2756 Level: intermediate 2757 2758 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA() 2759 2760 @*/ 2761 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2762 { 2763 PetscErrorCode ierr; 2764 2765 PetscFunctionBegin; 2766 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2767 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2768 PetscFunctionReturn(0); 2769 } 2770 2771 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2772 { 2773 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2774 PetscErrorCode ierr; 2775 2776 PetscFunctionBegin; 2777 if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2778 B->preallocated = PETSC_TRUE; 2779 2780 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2781 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2782 2783 if (b->lda <= 0) b->lda = B->rmap->n; 2784 2785 ierr = PetscIntMultError(b->lda,B->cmap->n,NULL);CHKERRQ(ierr); 2786 if (!data) { /* petsc-allocated storage */ 2787 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2788 ierr = PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);CHKERRQ(ierr); 2789 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2790 2791 b->user_alloc = PETSC_FALSE; 2792 } else { /* user-allocated storage */ 2793 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2794 b->v = data; 2795 b->user_alloc = PETSC_TRUE; 2796 } 2797 B->assembled = PETSC_TRUE; 2798 PetscFunctionReturn(0); 2799 } 2800 2801 #if defined(PETSC_HAVE_ELEMENTAL) 2802 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2803 { 2804 Mat mat_elemental; 2805 PetscErrorCode ierr; 2806 const PetscScalar *array; 2807 PetscScalar *v_colwise; 2808 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2809 2810 PetscFunctionBegin; 2811 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2812 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2813 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2814 k = 0; 2815 for (j=0; j<N; j++) { 2816 cols[j] = j; 2817 for (i=0; i<M; i++) { 2818 v_colwise[j*M+i] = array[k++]; 2819 } 2820 } 2821 for (i=0; i<M; i++) { 2822 rows[i] = i; 2823 } 2824 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2825 2826 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2827 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2828 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2829 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2830 2831 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2832 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2833 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2834 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2835 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2836 2837 if (reuse == MAT_INPLACE_MATRIX) { 2838 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2839 } else { 2840 *newmat = mat_elemental; 2841 } 2842 PetscFunctionReturn(0); 2843 } 2844 #endif 2845 2846 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 2847 { 2848 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2849 PetscBool data; 2850 2851 PetscFunctionBegin; 2852 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 2853 if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage"); 2854 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); 2855 b->lda = lda; 2856 PetscFunctionReturn(0); 2857 } 2858 2859 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2860 { 2861 PetscErrorCode ierr; 2862 PetscMPIInt size; 2863 2864 PetscFunctionBegin; 2865 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2866 if (size == 1) { 2867 if (scall == MAT_INITIAL_MATRIX) { 2868 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2869 } else { 2870 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2871 } 2872 } else { 2873 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2874 } 2875 PetscFunctionReturn(0); 2876 } 2877 2878 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2879 { 2880 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2881 PetscErrorCode ierr; 2882 2883 PetscFunctionBegin; 2884 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2885 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2886 if (!a->cvec) { 2887 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2888 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2889 } 2890 a->vecinuse = col + 1; 2891 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2892 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2893 *v = a->cvec; 2894 PetscFunctionReturn(0); 2895 } 2896 2897 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2898 { 2899 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2900 PetscErrorCode ierr; 2901 2902 PetscFunctionBegin; 2903 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2904 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2905 a->vecinuse = 0; 2906 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2907 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2908 *v = NULL; 2909 PetscFunctionReturn(0); 2910 } 2911 2912 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2913 { 2914 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2915 PetscErrorCode ierr; 2916 2917 PetscFunctionBegin; 2918 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2919 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2920 if (!a->cvec) { 2921 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2922 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2923 } 2924 a->vecinuse = col + 1; 2925 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2926 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2927 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr); 2928 *v = a->cvec; 2929 PetscFunctionReturn(0); 2930 } 2931 2932 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2933 { 2934 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2935 PetscErrorCode ierr; 2936 2937 PetscFunctionBegin; 2938 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2939 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2940 a->vecinuse = 0; 2941 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2942 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr); 2943 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2944 *v = NULL; 2945 PetscFunctionReturn(0); 2946 } 2947 2948 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2949 { 2950 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2951 PetscErrorCode ierr; 2952 2953 PetscFunctionBegin; 2954 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2955 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2956 if (!a->cvec) { 2957 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2958 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2959 } 2960 a->vecinuse = col + 1; 2961 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2962 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2963 *v = a->cvec; 2964 PetscFunctionReturn(0); 2965 } 2966 2967 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2968 { 2969 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2970 PetscErrorCode ierr; 2971 2972 PetscFunctionBegin; 2973 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2974 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2975 a->vecinuse = 0; 2976 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2977 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2978 *v = NULL; 2979 PetscFunctionReturn(0); 2980 } 2981 2982 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 2983 { 2984 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2985 PetscErrorCode ierr; 2986 2987 PetscFunctionBegin; 2988 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2989 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2990 if (a->cmat && cend-cbegin != a->cmat->cmap->N) { 2991 ierr = MatDestroy(&a->cmat);CHKERRQ(ierr); 2992 } 2993 if (!a->cmat) { 2994 ierr = MatCreateDense(PetscObjectComm((PetscObject)A),A->rmap->n,PETSC_DECIDE,A->rmap->N,cend-cbegin,(PetscScalar*)a->v + (size_t)cbegin * (size_t)a->lda,&a->cmat);CHKERRQ(ierr); 2995 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);CHKERRQ(ierr); 2996 } else { 2997 ierr = MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);CHKERRQ(ierr); 2998 } 2999 ierr = MatDenseSetLDA(a->cmat,a->lda);CHKERRQ(ierr); 3000 a->matinuse = cbegin + 1; 3001 *v = a->cmat; 3002 PetscFunctionReturn(0); 3003 } 3004 3005 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 3006 { 3007 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3008 PetscErrorCode ierr; 3009 3010 PetscFunctionBegin; 3011 if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 3012 if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 3013 if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 3014 a->matinuse = 0; 3015 ierr = MatDenseResetArray(a->cmat);CHKERRQ(ierr); 3016 *v = NULL; 3017 PetscFunctionReturn(0); 3018 } 3019 3020 /*MC 3021 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3022 3023 Options Database Keys: 3024 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 3025 3026 Level: beginner 3027 3028 .seealso: MatCreateSeqDense() 3029 3030 M*/ 3031 PetscErrorCode MatCreate_SeqDense(Mat B) 3032 { 3033 Mat_SeqDense *b; 3034 PetscErrorCode ierr; 3035 PetscMPIInt size; 3036 3037 PetscFunctionBegin; 3038 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3039 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3040 3041 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3042 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3043 B->data = (void*)b; 3044 3045 b->roworiented = PETSC_TRUE; 3046 3047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 3048 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);CHKERRQ(ierr); 3049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 3052 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 3053 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);CHKERRQ(ierr); 3054 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3055 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3056 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3057 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3058 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 3059 #if defined(PETSC_HAVE_ELEMENTAL) 3060 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 3061 #endif 3062 #if defined(PETSC_HAVE_SCALAPACK) 3063 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);CHKERRQ(ierr); 3064 #endif 3065 #if defined(PETSC_HAVE_CUDA) 3066 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 3067 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3068 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3069 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3070 #endif 3071 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 3072 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 3073 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3074 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3075 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3076 3077 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 3078 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 3079 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr); 3080 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr); 3081 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr); 3082 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr); 3083 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr); 3084 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr); 3085 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);CHKERRQ(ierr); 3086 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);CHKERRQ(ierr); 3087 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 3088 PetscFunctionReturn(0); 3089 } 3090 3091 /*@C 3092 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. 3093 3094 Not Collective 3095 3096 Input Parameters: 3097 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3098 - col - column index 3099 3100 Output Parameter: 3101 . vals - pointer to the data 3102 3103 Level: intermediate 3104 3105 .seealso: MatDenseRestoreColumn() 3106 @*/ 3107 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3108 { 3109 PetscErrorCode ierr; 3110 3111 PetscFunctionBegin; 3112 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3113 PetscValidLogicalCollectiveInt(A,col,2); 3114 PetscValidPointer(vals,3); 3115 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 3116 PetscFunctionReturn(0); 3117 } 3118 3119 /*@C 3120 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3121 3122 Not Collective 3123 3124 Input Parameter: 3125 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3126 3127 Output Parameter: 3128 . vals - pointer to the data 3129 3130 Level: intermediate 3131 3132 .seealso: MatDenseGetColumn() 3133 @*/ 3134 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3135 { 3136 PetscErrorCode ierr; 3137 3138 PetscFunctionBegin; 3139 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3140 PetscValidPointer(vals,2); 3141 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 3142 PetscFunctionReturn(0); 3143 } 3144 3145 /*@C 3146 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3147 3148 Collective 3149 3150 Input Parameters: 3151 + mat - the Mat object 3152 - col - the column index 3153 3154 Output Parameter: 3155 . v - the vector 3156 3157 Notes: 3158 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3159 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3160 3161 Level: intermediate 3162 3163 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3164 @*/ 3165 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3166 { 3167 PetscErrorCode ierr; 3168 3169 PetscFunctionBegin; 3170 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3171 PetscValidType(A,1); 3172 PetscValidLogicalCollectiveInt(A,col,2); 3173 PetscValidPointer(v,3); 3174 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3175 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); 3176 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3177 PetscFunctionReturn(0); 3178 } 3179 3180 /*@C 3181 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3182 3183 Collective 3184 3185 Input Parameters: 3186 + mat - the Mat object 3187 . col - the column index 3188 - v - the Vec object 3189 3190 Level: intermediate 3191 3192 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3193 @*/ 3194 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3195 { 3196 PetscErrorCode ierr; 3197 3198 PetscFunctionBegin; 3199 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3200 PetscValidType(A,1); 3201 PetscValidLogicalCollectiveInt(A,col,2); 3202 PetscValidPointer(v,3); 3203 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3204 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); 3205 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3206 PetscFunctionReturn(0); 3207 } 3208 3209 /*@C 3210 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3211 3212 Collective 3213 3214 Input Parameters: 3215 + mat - the Mat object 3216 - col - the column index 3217 3218 Output Parameter: 3219 . v - the vector 3220 3221 Notes: 3222 The vector is owned by PETSc and users cannot modify it. 3223 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3224 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3225 3226 Level: intermediate 3227 3228 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3229 @*/ 3230 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3231 { 3232 PetscErrorCode ierr; 3233 3234 PetscFunctionBegin; 3235 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3236 PetscValidType(A,1); 3237 PetscValidLogicalCollectiveInt(A,col,2); 3238 PetscValidPointer(v,3); 3239 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3240 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); 3241 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3242 PetscFunctionReturn(0); 3243 } 3244 3245 /*@C 3246 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3247 3248 Collective 3249 3250 Input Parameters: 3251 + mat - the Mat object 3252 . col - the column index 3253 - v - the Vec object 3254 3255 Level: intermediate 3256 3257 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite() 3258 @*/ 3259 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3260 { 3261 PetscErrorCode ierr; 3262 3263 PetscFunctionBegin; 3264 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3265 PetscValidType(A,1); 3266 PetscValidLogicalCollectiveInt(A,col,2); 3267 PetscValidPointer(v,3); 3268 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3269 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); 3270 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3271 PetscFunctionReturn(0); 3272 } 3273 3274 /*@C 3275 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3276 3277 Collective 3278 3279 Input Parameters: 3280 + mat - the Mat object 3281 - col - the column index 3282 3283 Output Parameter: 3284 . v - the vector 3285 3286 Notes: 3287 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3288 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3289 3290 Level: intermediate 3291 3292 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3293 @*/ 3294 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3295 { 3296 PetscErrorCode ierr; 3297 3298 PetscFunctionBegin; 3299 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3300 PetscValidType(A,1); 3301 PetscValidLogicalCollectiveInt(A,col,2); 3302 PetscValidPointer(v,3); 3303 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3304 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); 3305 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3306 PetscFunctionReturn(0); 3307 } 3308 3309 /*@C 3310 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3311 3312 Collective 3313 3314 Input Parameters: 3315 + mat - the Mat object 3316 . col - the column index 3317 - v - the Vec object 3318 3319 Level: intermediate 3320 3321 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead() 3322 @*/ 3323 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3324 { 3325 PetscErrorCode ierr; 3326 3327 PetscFunctionBegin; 3328 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3329 PetscValidType(A,1); 3330 PetscValidLogicalCollectiveInt(A,col,2); 3331 PetscValidPointer(v,3); 3332 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3333 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); 3334 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3335 PetscFunctionReturn(0); 3336 } 3337 3338 /*@C 3339 MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat. 3340 3341 Collective 3342 3343 Input Parameters: 3344 + mat - the Mat object 3345 . cbegin - the first index in the block 3346 - cend - the last index in the block 3347 3348 Output Parameter: 3349 . v - the matrix 3350 3351 Notes: 3352 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3353 3354 Level: intermediate 3355 3356 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix() 3357 @*/ 3358 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 3359 { 3360 PetscErrorCode ierr; 3361 3362 PetscFunctionBegin; 3363 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3364 PetscValidType(A,1); 3365 PetscValidLogicalCollectiveInt(A,cbegin,2); 3366 PetscValidLogicalCollectiveInt(A,cend,3); 3367 PetscValidPointer(v,4); 3368 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3369 if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N); 3370 if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N); 3371 ierr = PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));CHKERRQ(ierr); 3372 PetscFunctionReturn(0); 3373 } 3374 3375 /*@C 3376 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3377 3378 Collective 3379 3380 Input Parameters: 3381 + mat - the Mat object 3382 - v - the Mat object 3383 3384 Level: intermediate 3385 3386 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix() 3387 @*/ 3388 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3389 { 3390 PetscErrorCode ierr; 3391 3392 PetscFunctionBegin; 3393 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3394 PetscValidType(A,1); 3395 PetscValidPointer(v,2); 3396 ierr = PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));CHKERRQ(ierr); 3397 PetscFunctionReturn(0); 3398 } 3399