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