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