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