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