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