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,n,M; 1481 PetscScalar *v,tmp; 1482 1483 PetscFunctionBegin; 1484 m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1485 if (reuse == MAT_INPLACE_MATRIX && 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 { /* out-of-place transpose */ 1496 Mat tmat; 1497 Mat_SeqDense *tmatd; 1498 PetscScalar *v2; 1499 PetscInt M2; 1500 1501 if (reuse != MAT_REUSE_MATRIX) { 1502 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1503 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1504 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1505 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1506 } else tmat = *matout; 1507 1508 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1509 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1510 tmatd = (Mat_SeqDense*)tmat->data; 1511 M2 = tmatd->lda; 1512 for (j=0; j<n; j++) { 1513 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1514 } 1515 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1516 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1517 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1518 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1519 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1520 else { 1521 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1522 } 1523 } 1524 PetscFunctionReturn(0); 1525 } 1526 1527 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1528 { 1529 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1530 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1531 PetscInt i; 1532 const PetscScalar *v1,*v2; 1533 PetscErrorCode ierr; 1534 1535 PetscFunctionBegin; 1536 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1537 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1538 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1539 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1540 for (i=0; i<A1->cmap->n; i++) { 1541 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1542 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1543 v1 += mat1->lda; 1544 v2 += mat2->lda; 1545 } 1546 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1547 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1548 *flg = PETSC_TRUE; 1549 PetscFunctionReturn(0); 1550 } 1551 1552 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1553 { 1554 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1555 PetscInt i,n,len; 1556 PetscScalar *x; 1557 const PetscScalar *vv; 1558 PetscErrorCode ierr; 1559 1560 PetscFunctionBegin; 1561 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1562 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1563 len = PetscMin(A->rmap->n,A->cmap->n); 1564 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1565 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1566 for (i=0; i<len; i++) { 1567 x[i] = vv[i*mat->lda + i]; 1568 } 1569 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1570 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1571 PetscFunctionReturn(0); 1572 } 1573 1574 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1575 { 1576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1577 const PetscScalar *l,*r; 1578 PetscScalar x,*v,*vv; 1579 PetscErrorCode ierr; 1580 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1581 1582 PetscFunctionBegin; 1583 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1584 if (ll) { 1585 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1586 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1587 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1588 for (i=0; i<m; i++) { 1589 x = l[i]; 1590 v = vv + i; 1591 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1592 } 1593 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1594 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1595 } 1596 if (rr) { 1597 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1598 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1599 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1600 for (i=0; i<n; i++) { 1601 x = r[i]; 1602 v = vv + i*mat->lda; 1603 for (j=0; j<m; j++) (*v++) *= x; 1604 } 1605 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1606 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1607 } 1608 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1609 PetscFunctionReturn(0); 1610 } 1611 1612 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1613 { 1614 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1615 PetscScalar *v,*vv; 1616 PetscReal sum = 0.0; 1617 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1618 PetscErrorCode ierr; 1619 1620 PetscFunctionBegin; 1621 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1622 v = vv; 1623 if (type == NORM_FROBENIUS) { 1624 if (lda>m) { 1625 for (j=0; j<A->cmap->n; j++) { 1626 v = vv+j*lda; 1627 for (i=0; i<m; i++) { 1628 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1629 } 1630 } 1631 } else { 1632 #if defined(PETSC_USE_REAL___FP16) 1633 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1634 *nrm = BLASnrm2_(&cnt,v,&one); 1635 } 1636 #else 1637 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1638 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1639 } 1640 } 1641 *nrm = PetscSqrtReal(sum); 1642 #endif 1643 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1644 } else if (type == NORM_1) { 1645 *nrm = 0.0; 1646 for (j=0; j<A->cmap->n; j++) { 1647 v = vv + j*mat->lda; 1648 sum = 0.0; 1649 for (i=0; i<A->rmap->n; i++) { 1650 sum += PetscAbsScalar(*v); v++; 1651 } 1652 if (sum > *nrm) *nrm = sum; 1653 } 1654 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1655 } else if (type == NORM_INFINITY) { 1656 *nrm = 0.0; 1657 for (j=0; j<A->rmap->n; j++) { 1658 v = vv + j; 1659 sum = 0.0; 1660 for (i=0; i<A->cmap->n; i++) { 1661 sum += PetscAbsScalar(*v); v += mat->lda; 1662 } 1663 if (sum > *nrm) *nrm = sum; 1664 } 1665 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1666 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1667 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1672 { 1673 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1674 PetscErrorCode ierr; 1675 1676 PetscFunctionBegin; 1677 switch (op) { 1678 case MAT_ROW_ORIENTED: 1679 aij->roworiented = flg; 1680 break; 1681 case MAT_NEW_NONZERO_LOCATIONS: 1682 case MAT_NEW_NONZERO_LOCATION_ERR: 1683 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1684 case MAT_NEW_DIAGONALS: 1685 case MAT_KEEP_NONZERO_PATTERN: 1686 case MAT_IGNORE_OFF_PROC_ENTRIES: 1687 case MAT_USE_HASH_TABLE: 1688 case MAT_IGNORE_ZERO_ENTRIES: 1689 case MAT_IGNORE_LOWER_TRIANGULAR: 1690 case MAT_SORTED_FULL: 1691 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1692 break; 1693 case MAT_SPD: 1694 case MAT_SYMMETRIC: 1695 case MAT_STRUCTURALLY_SYMMETRIC: 1696 case MAT_HERMITIAN: 1697 case MAT_SYMMETRY_ETERNAL: 1698 /* These options are handled directly by MatSetOption() */ 1699 break; 1700 default: 1701 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1702 } 1703 PetscFunctionReturn(0); 1704 } 1705 1706 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1707 { 1708 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1709 PetscErrorCode ierr; 1710 PetscInt lda=l->lda,m=A->rmap->n,j; 1711 PetscScalar *v; 1712 1713 PetscFunctionBegin; 1714 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1715 if (lda>m) { 1716 for (j=0; j<A->cmap->n; j++) { 1717 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 1718 } 1719 } else { 1720 ierr = PetscArrayzero(v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1721 } 1722 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1723 PetscFunctionReturn(0); 1724 } 1725 1726 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1727 { 1728 PetscErrorCode ierr; 1729 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1730 PetscInt m = l->lda, n = A->cmap->n, i,j; 1731 PetscScalar *slot,*bb,*v; 1732 const PetscScalar *xx; 1733 1734 PetscFunctionBegin; 1735 if (PetscDefined(USE_DEBUG)) { 1736 for (i=0; i<N; i++) { 1737 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1738 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); 1739 } 1740 } 1741 if (!N) PetscFunctionReturn(0); 1742 1743 /* fix right hand side if needed */ 1744 if (x && b) { 1745 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1746 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1747 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1748 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1749 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1750 } 1751 1752 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1753 for (i=0; i<N; i++) { 1754 slot = v + rows[i]; 1755 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1756 } 1757 if (diag != 0.0) { 1758 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1759 for (i=0; i<N; i++) { 1760 slot = v + (m+1)*rows[i]; 1761 *slot = diag; 1762 } 1763 } 1764 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1765 PetscFunctionReturn(0); 1766 } 1767 1768 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1769 { 1770 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1771 1772 PetscFunctionBegin; 1773 *lda = mat->lda; 1774 PetscFunctionReturn(0); 1775 } 1776 1777 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 1778 { 1779 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1780 1781 PetscFunctionBegin; 1782 if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1783 *array = mat->v; 1784 PetscFunctionReturn(0); 1785 } 1786 1787 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 1788 { 1789 PetscFunctionBegin; 1790 *array = NULL; 1791 PetscFunctionReturn(0); 1792 } 1793 1794 /*@C 1795 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1796 1797 Not collective 1798 1799 Input Parameter: 1800 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1801 1802 Output Parameter: 1803 . lda - the leading dimension 1804 1805 Level: intermediate 1806 1807 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA() 1808 @*/ 1809 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1810 { 1811 PetscErrorCode ierr; 1812 1813 PetscFunctionBegin; 1814 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1815 PetscValidPointer(lda,2); 1816 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 /*@C 1821 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 1822 1823 Not collective 1824 1825 Input Parameter: 1826 + mat - a MATSEQDENSE or MATMPIDENSE matrix 1827 - lda - the leading dimension 1828 1829 Level: intermediate 1830 1831 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA() 1832 @*/ 1833 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 1834 { 1835 PetscErrorCode ierr; 1836 1837 PetscFunctionBegin; 1838 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1839 ierr = PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));CHKERRQ(ierr); 1840 PetscFunctionReturn(0); 1841 } 1842 1843 /*@C 1844 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 1845 1846 Logically Collective on Mat 1847 1848 Input Parameter: 1849 . mat - a dense matrix 1850 1851 Output Parameter: 1852 . array - pointer to the data 1853 1854 Level: intermediate 1855 1856 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1857 @*/ 1858 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1859 { 1860 PetscErrorCode ierr; 1861 1862 PetscFunctionBegin; 1863 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1864 PetscValidPointer(array,2); 1865 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1866 PetscFunctionReturn(0); 1867 } 1868 1869 /*@C 1870 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1871 1872 Logically Collective on Mat 1873 1874 Input Parameters: 1875 + mat - a dense matrix 1876 - array - pointer to the data 1877 1878 Level: intermediate 1879 1880 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1881 @*/ 1882 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1883 { 1884 PetscErrorCode ierr; 1885 1886 PetscFunctionBegin; 1887 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1888 PetscValidPointer(array,2); 1889 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1890 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1891 #if defined(PETSC_HAVE_CUDA) 1892 A->offloadmask = PETSC_OFFLOAD_CPU; 1893 #endif 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /*@C 1898 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 1899 1900 Not Collective 1901 1902 Input Parameter: 1903 . mat - a dense matrix 1904 1905 Output Parameter: 1906 . array - pointer to the data 1907 1908 Level: intermediate 1909 1910 .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1911 @*/ 1912 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1913 { 1914 PetscErrorCode ierr; 1915 1916 PetscFunctionBegin; 1917 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1918 PetscValidPointer(array,2); 1919 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1920 PetscFunctionReturn(0); 1921 } 1922 1923 /*@C 1924 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 1925 1926 Not Collective 1927 1928 Input Parameters: 1929 + mat - a dense matrix 1930 - array - pointer to the data 1931 1932 Level: intermediate 1933 1934 .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 1935 @*/ 1936 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1937 { 1938 PetscErrorCode ierr; 1939 1940 PetscFunctionBegin; 1941 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1942 PetscValidPointer(array,2); 1943 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1944 PetscFunctionReturn(0); 1945 } 1946 1947 /*@C 1948 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 1949 1950 Not Collective 1951 1952 Input Parameter: 1953 . mat - a dense matrix 1954 1955 Output Parameter: 1956 . array - pointer to the data 1957 1958 Level: intermediate 1959 1960 .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1961 @*/ 1962 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 1963 { 1964 PetscErrorCode ierr; 1965 1966 PetscFunctionBegin; 1967 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1968 PetscValidPointer(array,2); 1969 ierr = PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1970 PetscFunctionReturn(0); 1971 } 1972 1973 /*@C 1974 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 1975 1976 Not Collective 1977 1978 Input Parameters: 1979 + mat - a dense matrix 1980 - array - pointer to the data 1981 1982 Level: intermediate 1983 1984 .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1985 @*/ 1986 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 1987 { 1988 PetscErrorCode ierr; 1989 1990 PetscFunctionBegin; 1991 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1992 PetscValidPointer(array,2); 1993 ierr = PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1994 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1995 #if defined(PETSC_HAVE_CUDA) 1996 A->offloadmask = PETSC_OFFLOAD_CPU; 1997 #endif 1998 PetscFunctionReturn(0); 1999 } 2000 2001 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 2002 { 2003 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2004 PetscErrorCode ierr; 2005 PetscInt i,j,nrows,ncols,blda; 2006 const PetscInt *irow,*icol; 2007 PetscScalar *av,*bv,*v = mat->v; 2008 Mat newmat; 2009 2010 PetscFunctionBegin; 2011 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2012 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2013 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2014 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2015 2016 /* Check submatrixcall */ 2017 if (scall == MAT_REUSE_MATRIX) { 2018 PetscInt n_cols,n_rows; 2019 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2020 if (n_rows != nrows || n_cols != ncols) { 2021 /* resize the result matrix to match number of requested rows/columns */ 2022 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2023 } 2024 newmat = *B; 2025 } else { 2026 /* Create and fill new matrix */ 2027 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 2028 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2029 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 2030 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 2031 } 2032 2033 /* Now extract the data pointers and do the copy,column at a time */ 2034 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 2035 ierr = MatDenseGetLDA(newmat,&blda);CHKERRQ(ierr); 2036 for (i=0; i<ncols; i++) { 2037 av = v + mat->lda*icol[i]; 2038 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2039 bv += blda; 2040 } 2041 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 2042 2043 /* Assemble the matrices so that the correct flags are set */ 2044 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2045 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2046 2047 /* Free work space */ 2048 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2049 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2050 *B = newmat; 2051 PetscFunctionReturn(0); 2052 } 2053 2054 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2055 { 2056 PetscErrorCode ierr; 2057 PetscInt i; 2058 2059 PetscFunctionBegin; 2060 if (scall == MAT_INITIAL_MATRIX) { 2061 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2062 } 2063 2064 for (i=0; i<n; i++) { 2065 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2066 } 2067 PetscFunctionReturn(0); 2068 } 2069 2070 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2071 { 2072 PetscFunctionBegin; 2073 PetscFunctionReturn(0); 2074 } 2075 2076 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2077 { 2078 PetscFunctionBegin; 2079 PetscFunctionReturn(0); 2080 } 2081 2082 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2083 { 2084 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2085 PetscErrorCode ierr; 2086 const PetscScalar *va; 2087 PetscScalar *vb; 2088 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2089 2090 PetscFunctionBegin; 2091 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2092 if (A->ops->copy != B->ops->copy) { 2093 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2094 PetscFunctionReturn(0); 2095 } 2096 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2097 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2098 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2099 if (lda1>m || lda2>m) { 2100 for (j=0; j<n; j++) { 2101 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2102 } 2103 } else { 2104 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2105 } 2106 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2107 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2108 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2109 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2110 PetscFunctionReturn(0); 2111 } 2112 2113 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2114 { 2115 PetscErrorCode ierr; 2116 2117 PetscFunctionBegin; 2118 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2119 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2120 if (!A->preallocated) { 2121 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2122 } 2123 PetscFunctionReturn(0); 2124 } 2125 2126 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2127 { 2128 PetscInt i,nz = A->rmap->n*A->cmap->n; 2129 PetscScalar *aa; 2130 PetscErrorCode ierr; 2131 2132 PetscFunctionBegin; 2133 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2134 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2135 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2136 PetscFunctionReturn(0); 2137 } 2138 2139 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2140 { 2141 PetscInt i,nz = A->rmap->n*A->cmap->n; 2142 PetscScalar *aa; 2143 PetscErrorCode ierr; 2144 2145 PetscFunctionBegin; 2146 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2147 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2148 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2149 PetscFunctionReturn(0); 2150 } 2151 2152 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2153 { 2154 PetscInt i,nz = A->rmap->n*A->cmap->n; 2155 PetscScalar *aa; 2156 PetscErrorCode ierr; 2157 2158 PetscFunctionBegin; 2159 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2160 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2161 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2162 PetscFunctionReturn(0); 2163 } 2164 2165 /* ----------------------------------------------------------------*/ 2166 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2167 { 2168 PetscErrorCode ierr; 2169 PetscInt m=A->rmap->n,n=B->cmap->n; 2170 PetscBool cisdense; 2171 2172 PetscFunctionBegin; 2173 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2174 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2175 if (!cisdense) { 2176 PetscBool flg; 2177 2178 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2179 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2180 } 2181 ierr = MatSetUp(C);CHKERRQ(ierr); 2182 PetscFunctionReturn(0); 2183 } 2184 2185 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2186 { 2187 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2188 PetscBLASInt m,n,k; 2189 const PetscScalar *av,*bv; 2190 PetscScalar *cv; 2191 PetscScalar _DOne=1.0,_DZero=0.0; 2192 PetscErrorCode ierr; 2193 2194 PetscFunctionBegin; 2195 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2196 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2197 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2198 if (!m || !n || !k) PetscFunctionReturn(0); 2199 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2200 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2201 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2202 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2203 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2204 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2205 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2206 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2207 PetscFunctionReturn(0); 2208 } 2209 2210 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2211 { 2212 PetscErrorCode ierr; 2213 PetscInt m=A->rmap->n,n=B->rmap->n; 2214 PetscBool cisdense; 2215 2216 PetscFunctionBegin; 2217 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2218 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2219 if (!cisdense) { 2220 PetscBool flg; 2221 2222 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2223 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2224 } 2225 ierr = MatSetUp(C);CHKERRQ(ierr); 2226 PetscFunctionReturn(0); 2227 } 2228 2229 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2230 { 2231 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2232 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2233 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2234 const PetscScalar *av,*bv; 2235 PetscScalar *cv; 2236 PetscBLASInt m,n,k; 2237 PetscScalar _DOne=1.0,_DZero=0.0; 2238 PetscErrorCode ierr; 2239 2240 PetscFunctionBegin; 2241 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2242 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2243 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2244 if (!m || !n || !k) PetscFunctionReturn(0); 2245 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2246 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2247 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2248 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2249 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2250 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2251 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2252 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2253 PetscFunctionReturn(0); 2254 } 2255 2256 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2257 { 2258 PetscErrorCode ierr; 2259 PetscInt m=A->cmap->n,n=B->cmap->n; 2260 PetscBool cisdense; 2261 2262 PetscFunctionBegin; 2263 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2264 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2265 if (!cisdense) { 2266 PetscBool flg; 2267 2268 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2269 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2270 } 2271 ierr = MatSetUp(C);CHKERRQ(ierr); 2272 PetscFunctionReturn(0); 2273 } 2274 2275 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2276 { 2277 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2278 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2279 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2280 const PetscScalar *av,*bv; 2281 PetscScalar *cv; 2282 PetscBLASInt m,n,k; 2283 PetscScalar _DOne=1.0,_DZero=0.0; 2284 PetscErrorCode ierr; 2285 2286 PetscFunctionBegin; 2287 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2288 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2289 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2290 if (!m || !n || !k) PetscFunctionReturn(0); 2291 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2292 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2293 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2294 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2295 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2296 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2297 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2298 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2299 PetscFunctionReturn(0); 2300 } 2301 2302 /* ----------------------------------------------- */ 2303 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2304 { 2305 PetscFunctionBegin; 2306 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2307 C->ops->productsymbolic = MatProductSymbolic_AB; 2308 PetscFunctionReturn(0); 2309 } 2310 2311 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2312 { 2313 PetscFunctionBegin; 2314 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2315 C->ops->productsymbolic = MatProductSymbolic_AtB; 2316 PetscFunctionReturn(0); 2317 } 2318 2319 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2320 { 2321 PetscFunctionBegin; 2322 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2323 C->ops->productsymbolic = MatProductSymbolic_ABt; 2324 PetscFunctionReturn(0); 2325 } 2326 2327 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2328 { 2329 PetscErrorCode ierr; 2330 Mat_Product *product = C->product; 2331 2332 PetscFunctionBegin; 2333 switch (product->type) { 2334 case MATPRODUCT_AB: 2335 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2336 break; 2337 case MATPRODUCT_AtB: 2338 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2339 break; 2340 case MATPRODUCT_ABt: 2341 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2342 break; 2343 default: 2344 break; 2345 } 2346 PetscFunctionReturn(0); 2347 } 2348 /* ----------------------------------------------- */ 2349 2350 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2351 { 2352 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2353 PetscErrorCode ierr; 2354 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2355 PetscScalar *x; 2356 const PetscScalar *aa; 2357 2358 PetscFunctionBegin; 2359 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2360 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2361 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2362 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2363 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2364 for (i=0; i<m; i++) { 2365 x[i] = aa[i]; if (idx) idx[i] = 0; 2366 for (j=1; j<n; j++) { 2367 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2368 } 2369 } 2370 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2371 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2372 PetscFunctionReturn(0); 2373 } 2374 2375 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2376 { 2377 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2378 PetscErrorCode ierr; 2379 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2380 PetscScalar *x; 2381 PetscReal atmp; 2382 const PetscScalar *aa; 2383 2384 PetscFunctionBegin; 2385 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2386 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2387 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2388 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2389 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2390 for (i=0; i<m; i++) { 2391 x[i] = PetscAbsScalar(aa[i]); 2392 for (j=1; j<n; j++) { 2393 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2394 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2395 } 2396 } 2397 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2398 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2399 PetscFunctionReturn(0); 2400 } 2401 2402 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2403 { 2404 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2405 PetscErrorCode ierr; 2406 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2407 PetscScalar *x; 2408 const PetscScalar *aa; 2409 2410 PetscFunctionBegin; 2411 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2412 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2413 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2414 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2415 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2416 for (i=0; i<m; i++) { 2417 x[i] = aa[i]; if (idx) idx[i] = 0; 2418 for (j=1; j<n; j++) { 2419 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2420 } 2421 } 2422 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2423 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2424 PetscFunctionReturn(0); 2425 } 2426 2427 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2428 { 2429 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2430 PetscErrorCode ierr; 2431 PetscScalar *x; 2432 const PetscScalar *aa; 2433 2434 PetscFunctionBegin; 2435 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2436 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2437 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2438 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2439 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2440 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2441 PetscFunctionReturn(0); 2442 } 2443 2444 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2445 { 2446 PetscErrorCode ierr; 2447 PetscInt i,j,m,n; 2448 const PetscScalar *a; 2449 2450 PetscFunctionBegin; 2451 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2452 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2453 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2454 if (type == NORM_2) { 2455 for (i=0; i<n; i++) { 2456 for (j=0; j<m; j++) { 2457 norms[i] += PetscAbsScalar(a[j]*a[j]); 2458 } 2459 a += m; 2460 } 2461 } else if (type == NORM_1) { 2462 for (i=0; i<n; i++) { 2463 for (j=0; j<m; j++) { 2464 norms[i] += PetscAbsScalar(a[j]); 2465 } 2466 a += m; 2467 } 2468 } else if (type == NORM_INFINITY) { 2469 for (i=0; i<n; i++) { 2470 for (j=0; j<m; j++) { 2471 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2472 } 2473 a += m; 2474 } 2475 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2476 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2477 if (type == NORM_2) { 2478 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2479 } 2480 PetscFunctionReturn(0); 2481 } 2482 2483 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2484 { 2485 PetscErrorCode ierr; 2486 PetscScalar *a; 2487 PetscInt lda,m,n,i,j; 2488 2489 PetscFunctionBegin; 2490 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2491 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2492 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2493 for (j=0; j<n; j++) { 2494 for (i=0; i<m; i++) { 2495 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2496 } 2497 } 2498 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2499 PetscFunctionReturn(0); 2500 } 2501 2502 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2503 { 2504 PetscFunctionBegin; 2505 *missing = PETSC_FALSE; 2506 PetscFunctionReturn(0); 2507 } 2508 2509 /* vals is not const */ 2510 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2511 { 2512 PetscErrorCode ierr; 2513 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2514 PetscScalar *v; 2515 2516 PetscFunctionBegin; 2517 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2518 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2519 *vals = v+col*a->lda; 2520 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2521 PetscFunctionReturn(0); 2522 } 2523 2524 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2525 { 2526 PetscFunctionBegin; 2527 *vals = 0; /* user cannot accidently use the array later */ 2528 PetscFunctionReturn(0); 2529 } 2530 2531 /* -------------------------------------------------------------------*/ 2532 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2533 MatGetRow_SeqDense, 2534 MatRestoreRow_SeqDense, 2535 MatMult_SeqDense, 2536 /* 4*/ MatMultAdd_SeqDense, 2537 MatMultTranspose_SeqDense, 2538 MatMultTransposeAdd_SeqDense, 2539 0, 2540 0, 2541 0, 2542 /* 10*/ 0, 2543 MatLUFactor_SeqDense, 2544 MatCholeskyFactor_SeqDense, 2545 MatSOR_SeqDense, 2546 MatTranspose_SeqDense, 2547 /* 15*/ MatGetInfo_SeqDense, 2548 MatEqual_SeqDense, 2549 MatGetDiagonal_SeqDense, 2550 MatDiagonalScale_SeqDense, 2551 MatNorm_SeqDense, 2552 /* 20*/ MatAssemblyBegin_SeqDense, 2553 MatAssemblyEnd_SeqDense, 2554 MatSetOption_SeqDense, 2555 MatZeroEntries_SeqDense, 2556 /* 24*/ MatZeroRows_SeqDense, 2557 0, 2558 0, 2559 0, 2560 0, 2561 /* 29*/ MatSetUp_SeqDense, 2562 0, 2563 0, 2564 0, 2565 0, 2566 /* 34*/ MatDuplicate_SeqDense, 2567 0, 2568 0, 2569 0, 2570 0, 2571 /* 39*/ MatAXPY_SeqDense, 2572 MatCreateSubMatrices_SeqDense, 2573 0, 2574 MatGetValues_SeqDense, 2575 MatCopy_SeqDense, 2576 /* 44*/ MatGetRowMax_SeqDense, 2577 MatScale_SeqDense, 2578 MatShift_Basic, 2579 0, 2580 MatZeroRowsColumns_SeqDense, 2581 /* 49*/ MatSetRandom_SeqDense, 2582 0, 2583 0, 2584 0, 2585 0, 2586 /* 54*/ 0, 2587 0, 2588 0, 2589 0, 2590 0, 2591 /* 59*/ 0, 2592 MatDestroy_SeqDense, 2593 MatView_SeqDense, 2594 0, 2595 0, 2596 /* 64*/ 0, 2597 0, 2598 0, 2599 0, 2600 0, 2601 /* 69*/ MatGetRowMaxAbs_SeqDense, 2602 0, 2603 0, 2604 0, 2605 0, 2606 /* 74*/ 0, 2607 0, 2608 0, 2609 0, 2610 0, 2611 /* 79*/ 0, 2612 0, 2613 0, 2614 0, 2615 /* 83*/ MatLoad_SeqDense, 2616 MatIsSymmetric_SeqDense, 2617 MatIsHermitian_SeqDense, 2618 0, 2619 0, 2620 0, 2621 /* 89*/ 0, 2622 0, 2623 MatMatMultNumeric_SeqDense_SeqDense, 2624 0, 2625 0, 2626 /* 94*/ 0, 2627 0, 2628 0, 2629 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2630 0, 2631 /* 99*/ MatProductSetFromOptions_SeqDense, 2632 0, 2633 0, 2634 MatConjugate_SeqDense, 2635 0, 2636 /*104*/ 0, 2637 MatRealPart_SeqDense, 2638 MatImaginaryPart_SeqDense, 2639 0, 2640 0, 2641 /*109*/ 0, 2642 0, 2643 MatGetRowMin_SeqDense, 2644 MatGetColumnVector_SeqDense, 2645 MatMissingDiagonal_SeqDense, 2646 /*114*/ 0, 2647 0, 2648 0, 2649 0, 2650 0, 2651 /*119*/ 0, 2652 0, 2653 0, 2654 0, 2655 0, 2656 /*124*/ 0, 2657 MatGetColumnNorms_SeqDense, 2658 0, 2659 0, 2660 0, 2661 /*129*/ 0, 2662 0, 2663 0, 2664 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2665 0, 2666 /*134*/ 0, 2667 0, 2668 0, 2669 0, 2670 0, 2671 /*139*/ 0, 2672 0, 2673 0, 2674 0, 2675 0, 2676 MatCreateMPIMatConcatenateSeqMat_SeqDense, 2677 /*145*/ 0, 2678 0, 2679 0 2680 }; 2681 2682 /*@C 2683 MatCreateSeqDense - Creates a sequential dense matrix that 2684 is stored in column major order (the usual Fortran 77 manner). Many 2685 of the matrix operations use the BLAS and LAPACK routines. 2686 2687 Collective 2688 2689 Input Parameters: 2690 + comm - MPI communicator, set to PETSC_COMM_SELF 2691 . m - number of rows 2692 . n - number of columns 2693 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2694 to control all matrix memory allocation. 2695 2696 Output Parameter: 2697 . A - the matrix 2698 2699 Notes: 2700 The data input variable is intended primarily for Fortran programmers 2701 who wish to allocate their own matrix memory space. Most users should 2702 set data=NULL. 2703 2704 Level: intermediate 2705 2706 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2707 @*/ 2708 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2709 { 2710 PetscErrorCode ierr; 2711 2712 PetscFunctionBegin; 2713 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2714 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2715 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2716 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2717 PetscFunctionReturn(0); 2718 } 2719 2720 /*@C 2721 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2722 2723 Collective 2724 2725 Input Parameters: 2726 + B - the matrix 2727 - data - the array (or NULL) 2728 2729 Notes: 2730 The data input variable is intended primarily for Fortran programmers 2731 who wish to allocate their own matrix memory space. Most users should 2732 need not call this routine. 2733 2734 Level: intermediate 2735 2736 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA() 2737 2738 @*/ 2739 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2740 { 2741 PetscErrorCode ierr; 2742 2743 PetscFunctionBegin; 2744 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2745 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2746 PetscFunctionReturn(0); 2747 } 2748 2749 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2750 { 2751 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2752 PetscErrorCode ierr; 2753 2754 PetscFunctionBegin; 2755 if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2756 B->preallocated = PETSC_TRUE; 2757 2758 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2759 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2760 2761 if (b->lda <= 0) b->lda = B->rmap->n; 2762 2763 ierr = PetscIntMultError(b->lda,B->cmap->n,NULL);CHKERRQ(ierr); 2764 if (!data) { /* petsc-allocated storage */ 2765 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2766 ierr = PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);CHKERRQ(ierr); 2767 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 2768 2769 b->user_alloc = PETSC_FALSE; 2770 } else { /* user-allocated storage */ 2771 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2772 b->v = data; 2773 b->user_alloc = PETSC_TRUE; 2774 } 2775 B->assembled = PETSC_TRUE; 2776 PetscFunctionReturn(0); 2777 } 2778 2779 #if defined(PETSC_HAVE_ELEMENTAL) 2780 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2781 { 2782 Mat mat_elemental; 2783 PetscErrorCode ierr; 2784 const PetscScalar *array; 2785 PetscScalar *v_colwise; 2786 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2787 2788 PetscFunctionBegin; 2789 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2790 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2791 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2792 k = 0; 2793 for (j=0; j<N; j++) { 2794 cols[j] = j; 2795 for (i=0; i<M; i++) { 2796 v_colwise[j*M+i] = array[k++]; 2797 } 2798 } 2799 for (i=0; i<M; i++) { 2800 rows[i] = i; 2801 } 2802 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2803 2804 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2805 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2806 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2807 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2808 2809 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2810 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2811 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2812 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2813 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2814 2815 if (reuse == MAT_INPLACE_MATRIX) { 2816 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2817 } else { 2818 *newmat = mat_elemental; 2819 } 2820 PetscFunctionReturn(0); 2821 } 2822 #endif 2823 2824 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 2825 { 2826 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2827 PetscBool data; 2828 2829 PetscFunctionBegin; 2830 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 2831 if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage"); 2832 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); 2833 b->lda = lda; 2834 PetscFunctionReturn(0); 2835 } 2836 2837 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2838 { 2839 PetscErrorCode ierr; 2840 PetscMPIInt size; 2841 2842 PetscFunctionBegin; 2843 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2844 if (size == 1) { 2845 if (scall == MAT_INITIAL_MATRIX) { 2846 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2847 } else { 2848 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2849 } 2850 } else { 2851 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2852 } 2853 PetscFunctionReturn(0); 2854 } 2855 2856 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2857 { 2858 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2859 PetscErrorCode ierr; 2860 2861 PetscFunctionBegin; 2862 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2863 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2864 if (!a->cvec) { 2865 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2866 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2867 } 2868 a->vecinuse = col + 1; 2869 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2870 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2871 *v = a->cvec; 2872 PetscFunctionReturn(0); 2873 } 2874 2875 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 2876 { 2877 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2878 PetscErrorCode ierr; 2879 2880 PetscFunctionBegin; 2881 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2882 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2883 a->vecinuse = 0; 2884 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2885 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2886 *v = NULL; 2887 PetscFunctionReturn(0); 2888 } 2889 2890 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2891 { 2892 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2893 PetscErrorCode ierr; 2894 2895 PetscFunctionBegin; 2896 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2897 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2898 if (!a->cvec) { 2899 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2900 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2901 } 2902 a->vecinuse = col + 1; 2903 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2904 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2905 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr); 2906 *v = a->cvec; 2907 PetscFunctionReturn(0); 2908 } 2909 2910 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 2911 { 2912 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2913 PetscErrorCode ierr; 2914 2915 PetscFunctionBegin; 2916 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2917 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2918 a->vecinuse = 0; 2919 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 2920 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr); 2921 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2922 *v = NULL; 2923 PetscFunctionReturn(0); 2924 } 2925 2926 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2927 { 2928 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2929 PetscErrorCode ierr; 2930 2931 PetscFunctionBegin; 2932 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2933 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2934 if (!a->cvec) { 2935 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 2936 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 2937 } 2938 a->vecinuse = col + 1; 2939 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2940 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 2941 *v = a->cvec; 2942 PetscFunctionReturn(0); 2943 } 2944 2945 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 2946 { 2947 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2948 PetscErrorCode ierr; 2949 2950 PetscFunctionBegin; 2951 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 2952 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 2953 a->vecinuse = 0; 2954 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 2955 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 2956 *v = NULL; 2957 PetscFunctionReturn(0); 2958 } 2959 2960 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 2961 { 2962 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2963 PetscErrorCode ierr; 2964 2965 PetscFunctionBegin; 2966 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 2967 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2968 if (a->cmat && cend-cbegin != a->cmat->cmap->N) { 2969 ierr = MatDestroy(&a->cmat);CHKERRQ(ierr); 2970 } 2971 if (!a->cmat) { 2972 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); 2973 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);CHKERRQ(ierr); 2974 } else { 2975 ierr = MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);CHKERRQ(ierr); 2976 } 2977 ierr = MatDenseSetLDA(a->cmat,a->lda);CHKERRQ(ierr); 2978 a->matinuse = cbegin + 1; 2979 *v = a->cmat; 2980 PetscFunctionReturn(0); 2981 } 2982 2983 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 2984 { 2985 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2986 PetscErrorCode ierr; 2987 2988 PetscFunctionBegin; 2989 if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 2990 if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 2991 if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 2992 a->matinuse = 0; 2993 ierr = MatDenseResetArray(a->cmat);CHKERRQ(ierr); 2994 *v = NULL; 2995 PetscFunctionReturn(0); 2996 } 2997 2998 /*MC 2999 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3000 3001 Options Database Keys: 3002 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 3003 3004 Level: beginner 3005 3006 .seealso: MatCreateSeqDense() 3007 3008 M*/ 3009 PetscErrorCode MatCreate_SeqDense(Mat B) 3010 { 3011 Mat_SeqDense *b; 3012 PetscErrorCode ierr; 3013 PetscMPIInt size; 3014 3015 PetscFunctionBegin; 3016 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3017 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3018 3019 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3020 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3021 B->data = (void*)b; 3022 3023 b->roworiented = PETSC_TRUE; 3024 3025 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 3026 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);CHKERRQ(ierr); 3027 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3028 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3029 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 3030 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 3031 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);CHKERRQ(ierr); 3032 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3033 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3034 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3035 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3036 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 3037 #if defined(PETSC_HAVE_ELEMENTAL) 3038 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 3039 #endif 3040 #if defined(PETSC_HAVE_SCALAPACK) 3041 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);CHKERRQ(ierr); 3042 #endif 3043 #if defined(PETSC_HAVE_CUDA) 3044 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 3045 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3046 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3048 #endif 3049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 3050 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 3051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3052 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3053 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3054 3055 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 3056 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 3057 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr); 3058 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr); 3059 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr); 3060 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr); 3061 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr); 3062 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr); 3063 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);CHKERRQ(ierr); 3064 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);CHKERRQ(ierr); 3065 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 3066 PetscFunctionReturn(0); 3067 } 3068 3069 /*@C 3070 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. 3071 3072 Not Collective 3073 3074 Input Parameters: 3075 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3076 - col - column index 3077 3078 Output Parameter: 3079 . vals - pointer to the data 3080 3081 Level: intermediate 3082 3083 .seealso: MatDenseRestoreColumn() 3084 @*/ 3085 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3086 { 3087 PetscErrorCode ierr; 3088 3089 PetscFunctionBegin; 3090 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3091 PetscValidLogicalCollectiveInt(A,col,2); 3092 PetscValidPointer(vals,3); 3093 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 3094 PetscFunctionReturn(0); 3095 } 3096 3097 /*@C 3098 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3099 3100 Not Collective 3101 3102 Input Parameter: 3103 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3104 3105 Output Parameter: 3106 . vals - pointer to the data 3107 3108 Level: intermediate 3109 3110 .seealso: MatDenseGetColumn() 3111 @*/ 3112 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3113 { 3114 PetscErrorCode ierr; 3115 3116 PetscFunctionBegin; 3117 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3118 PetscValidPointer(vals,2); 3119 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 3120 PetscFunctionReturn(0); 3121 } 3122 3123 /*@C 3124 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3125 3126 Collective 3127 3128 Input Parameters: 3129 + mat - the Mat object 3130 - col - the column index 3131 3132 Output Parameter: 3133 . v - the vector 3134 3135 Notes: 3136 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3137 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3138 3139 Level: intermediate 3140 3141 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3142 @*/ 3143 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3144 { 3145 PetscErrorCode ierr; 3146 3147 PetscFunctionBegin; 3148 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3149 PetscValidType(A,1); 3150 PetscValidLogicalCollectiveInt(A,col,2); 3151 PetscValidPointer(v,3); 3152 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3153 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); 3154 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3155 PetscFunctionReturn(0); 3156 } 3157 3158 /*@C 3159 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3160 3161 Collective 3162 3163 Input Parameters: 3164 + mat - the Mat object 3165 . col - the column index 3166 - v - the Vec object 3167 3168 Level: intermediate 3169 3170 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3171 @*/ 3172 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3173 { 3174 PetscErrorCode ierr; 3175 3176 PetscFunctionBegin; 3177 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3178 PetscValidType(A,1); 3179 PetscValidLogicalCollectiveInt(A,col,2); 3180 PetscValidPointer(v,3); 3181 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3182 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); 3183 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3184 PetscFunctionReturn(0); 3185 } 3186 3187 /*@C 3188 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3189 3190 Collective 3191 3192 Input Parameters: 3193 + mat - the Mat object 3194 - col - the column index 3195 3196 Output Parameter: 3197 . v - the vector 3198 3199 Notes: 3200 The vector is owned by PETSc and users cannot modify it. 3201 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3202 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3203 3204 Level: intermediate 3205 3206 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3207 @*/ 3208 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3209 { 3210 PetscErrorCode ierr; 3211 3212 PetscFunctionBegin; 3213 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3214 PetscValidType(A,1); 3215 PetscValidLogicalCollectiveInt(A,col,2); 3216 PetscValidPointer(v,3); 3217 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3218 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); 3219 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3220 PetscFunctionReturn(0); 3221 } 3222 3223 /*@C 3224 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3225 3226 Collective 3227 3228 Input Parameters: 3229 + mat - the Mat object 3230 . col - the column index 3231 - v - the Vec object 3232 3233 Level: intermediate 3234 3235 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite() 3236 @*/ 3237 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3238 { 3239 PetscErrorCode ierr; 3240 3241 PetscFunctionBegin; 3242 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3243 PetscValidType(A,1); 3244 PetscValidLogicalCollectiveInt(A,col,2); 3245 PetscValidPointer(v,3); 3246 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3247 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); 3248 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3249 PetscFunctionReturn(0); 3250 } 3251 3252 /*@C 3253 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3254 3255 Collective 3256 3257 Input Parameters: 3258 + mat - the Mat object 3259 - col - the column index 3260 3261 Output Parameter: 3262 . v - the vector 3263 3264 Notes: 3265 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3266 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3267 3268 Level: intermediate 3269 3270 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3271 @*/ 3272 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3273 { 3274 PetscErrorCode ierr; 3275 3276 PetscFunctionBegin; 3277 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3278 PetscValidType(A,1); 3279 PetscValidLogicalCollectiveInt(A,col,2); 3280 PetscValidPointer(v,3); 3281 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3282 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); 3283 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3284 PetscFunctionReturn(0); 3285 } 3286 3287 /*@C 3288 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3289 3290 Collective 3291 3292 Input Parameters: 3293 + mat - the Mat object 3294 . col - the column index 3295 - v - the Vec object 3296 3297 Level: intermediate 3298 3299 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead() 3300 @*/ 3301 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3302 { 3303 PetscErrorCode ierr; 3304 3305 PetscFunctionBegin; 3306 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3307 PetscValidType(A,1); 3308 PetscValidLogicalCollectiveInt(A,col,2); 3309 PetscValidPointer(v,3); 3310 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3311 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); 3312 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3313 PetscFunctionReturn(0); 3314 } 3315 3316 /*@C 3317 MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat. 3318 3319 Collective 3320 3321 Input Parameters: 3322 + mat - the Mat object 3323 . cbegin - the first index in the block 3324 - cend - the last index in the block 3325 3326 Output Parameter: 3327 . v - the matrix 3328 3329 Notes: 3330 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3331 3332 Level: intermediate 3333 3334 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix() 3335 @*/ 3336 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 3337 { 3338 PetscErrorCode ierr; 3339 3340 PetscFunctionBegin; 3341 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3342 PetscValidType(A,1); 3343 PetscValidLogicalCollectiveInt(A,cbegin,2); 3344 PetscValidLogicalCollectiveInt(A,cend,3); 3345 PetscValidPointer(v,4); 3346 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3347 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); 3348 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); 3349 ierr = PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));CHKERRQ(ierr); 3350 PetscFunctionReturn(0); 3351 } 3352 3353 /*@C 3354 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3355 3356 Collective 3357 3358 Input Parameters: 3359 + mat - the Mat object 3360 - v - the Mat object 3361 3362 Level: intermediate 3363 3364 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix() 3365 @*/ 3366 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3367 { 3368 PetscErrorCode ierr; 3369 3370 PetscFunctionBegin; 3371 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3372 PetscValidType(A,1); 3373 PetscValidPointer(v,2); 3374 ierr = PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));CHKERRQ(ierr); 3375 PetscFunctionReturn(0); 3376 } 3377