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