1 2 /* 3 Defines the basic matrix operations for sequential dense. 4 */ 5 6 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 9 #include <../src/mat/impls/aij/seq/aij.h> 10 11 PetscErrorCode MatSeqDenseSymmetrize_Private(Mat A, PetscBool hermitian) 12 { 13 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 14 PetscInt j, k, n = A->rmap->n; 15 PetscScalar *v; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 20 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 21 if (!hermitian) { 22 for (k=0;k<n;k++) { 23 for (j=k;j<n;j++) { 24 v[j*mat->lda + k] = v[k*mat->lda + j]; 25 } 26 } 27 } else { 28 for (k=0;k<n;k++) { 29 for (j=k;j<n;j++) { 30 v[j*mat->lda + k] = PetscConj(v[k*mat->lda + j]); 31 } 32 } 33 } 34 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 35 PetscFunctionReturn(0); 36 } 37 38 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 39 { 40 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 41 PetscErrorCode ierr; 42 PetscBLASInt info,n; 43 44 PetscFunctionBegin; 45 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 46 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 47 if (A->factortype == MAT_FACTOR_LU) { 48 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 49 if (!mat->fwork) { 50 mat->lfwork = n; 51 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 52 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 53 } 54 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 55 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 56 ierr = PetscFPTrapPop();CHKERRQ(ierr); 57 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 58 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 59 if (A->spd) { 60 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 61 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&n,mat->v,&mat->lda,&info)); 62 ierr = PetscFPTrapPop();CHKERRQ(ierr); 63 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 64 #if defined(PETSC_USE_COMPLEX) 65 } else if (A->hermitian) { 66 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 67 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 68 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 69 PetscStackCallBLAS("LAPACKhetri",LAPACKhetri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 70 ierr = PetscFPTrapPop();CHKERRQ(ierr); 71 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_TRUE);CHKERRQ(ierr); 72 #endif 73 } else { /* symmetric case */ 74 if (!mat->pivots) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Pivots not present"); 75 if (!mat->fwork) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Fwork not present"); 76 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 77 PetscStackCallBLAS("LAPACKsytri",LAPACKsytri_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&info)); 78 ierr = PetscFPTrapPop();CHKERRQ(ierr); 79 ierr = MatSeqDenseSymmetrize_Private(A,PETSC_FALSE);CHKERRQ(ierr); 80 } 81 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad Inversion: zero pivot in row %D",(PetscInt)info-1); 82 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 84 85 A->ops->solve = NULL; 86 A->ops->matsolve = NULL; 87 A->ops->solvetranspose = NULL; 88 A->ops->matsolvetranspose = NULL; 89 A->ops->solveadd = NULL; 90 A->ops->solvetransposeadd = NULL; 91 A->factortype = MAT_FACTOR_NONE; 92 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 97 { 98 PetscErrorCode ierr; 99 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 100 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 101 PetscScalar *slot,*bb,*v; 102 const PetscScalar *xx; 103 104 PetscFunctionBegin; 105 if (PetscDefined(USE_DEBUG)) { 106 for (i=0; i<N; i++) { 107 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 108 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 109 if (rows[i] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Col %D requested to be zeroed greater than or equal number of cols %D",rows[i],A->cmap->n); 110 } 111 } 112 if (!N) PetscFunctionReturn(0); 113 114 /* fix right hand side if needed */ 115 if (x && b) { 116 Vec xt; 117 118 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 119 ierr = VecDuplicate(x,&xt);CHKERRQ(ierr); 120 ierr = VecCopy(x,xt);CHKERRQ(ierr); 121 ierr = VecScale(xt,-1.0);CHKERRQ(ierr); 122 ierr = MatMultAdd(A,xt,b,b);CHKERRQ(ierr); 123 ierr = VecDestroy(&xt);CHKERRQ(ierr); 124 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 125 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 126 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 127 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 128 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 129 } 130 131 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 132 for (i=0; i<N; i++) { 133 slot = v + rows[i]*m; 134 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 135 } 136 for (i=0; i<N; i++) { 137 slot = v + rows[i]; 138 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 139 } 140 if (diag != 0.0) { 141 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 142 for (i=0; i<N; i++) { 143 slot = v + (m+1)*rows[i]; 144 *slot = diag; 145 } 146 } 147 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 152 { 153 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 if (c->ptapwork) { 158 ierr = (*C->ops->matmultnumeric)(A,P,c->ptapwork);CHKERRQ(ierr); 159 ierr = (*C->ops->transposematmultnumeric)(P,c->ptapwork,C);CHKERRQ(ierr); 160 } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"Must call MatPtAPSymbolic_SeqDense_SeqDense() first"); 161 PetscFunctionReturn(0); 162 } 163 164 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat C) 165 { 166 Mat_SeqDense *c; 167 PetscBool cisdense; 168 PetscErrorCode ierr; 169 170 PetscFunctionBegin; 171 ierr = MatSetSizes(C,P->cmap->n,P->cmap->n,P->cmap->N,P->cmap->N);CHKERRQ(ierr); 172 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 173 if (!cisdense) { 174 PetscBool flg; 175 176 ierr = PetscObjectTypeCompare((PetscObject)P,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 177 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 178 } 179 ierr = MatSetUp(C);CHKERRQ(ierr); 180 c = (Mat_SeqDense*)C->data; 181 ierr = MatCreate(PetscObjectComm((PetscObject)A),&c->ptapwork);CHKERRQ(ierr); 182 ierr = MatSetSizes(c->ptapwork,A->rmap->n,P->cmap->n,A->rmap->N,P->cmap->N);CHKERRQ(ierr); 183 ierr = MatSetType(c->ptapwork,((PetscObject)C)->type_name);CHKERRQ(ierr); 184 ierr = MatSetUp(c->ptapwork);CHKERRQ(ierr); 185 PetscFunctionReturn(0); 186 } 187 188 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 189 { 190 Mat B = NULL; 191 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 192 Mat_SeqDense *b; 193 PetscErrorCode ierr; 194 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 195 const MatScalar *av; 196 PetscBool isseqdense; 197 198 PetscFunctionBegin; 199 if (reuse == MAT_REUSE_MATRIX) { 200 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 201 if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type_name); 202 } 203 if (reuse != MAT_REUSE_MATRIX) { 204 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 205 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 206 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 207 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr); 208 b = (Mat_SeqDense*)(B->data); 209 } else { 210 b = (Mat_SeqDense*)((*newmat)->data); 211 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr); 212 } 213 ierr = MatSeqAIJGetArrayRead(A,&av);CHKERRQ(ierr); 214 for (i=0; i<m; i++) { 215 PetscInt j; 216 for (j=0;j<ai[1]-ai[0];j++) { 217 b->v[*aj*m+i] = *av; 218 aj++; 219 av++; 220 } 221 ai++; 222 } 223 ierr = MatSeqAIJRestoreArrayRead(A,&av);CHKERRQ(ierr); 224 225 if (reuse == MAT_INPLACE_MATRIX) { 226 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 227 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 228 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 229 } else { 230 if (B) *newmat = B; 231 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 232 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 233 } 234 PetscFunctionReturn(0); 235 } 236 237 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 238 { 239 Mat B = NULL; 240 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 241 PetscErrorCode ierr; 242 PetscInt i, j; 243 PetscInt *rows, *nnz; 244 MatScalar *aa = a->v, *vals; 245 246 PetscFunctionBegin; 247 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 248 if (reuse != MAT_REUSE_MATRIX) { 249 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 250 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 251 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 252 for (j=0; j<A->cmap->n; j++) { 253 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) ++nnz[i]; 254 aa += a->lda; 255 } 256 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 257 } else B = *newmat; 258 aa = a->v; 259 for (j=0; j<A->cmap->n; j++) { 260 PetscInt numRows = 0; 261 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || (i == j && A->cmap->n == A->rmap->n)) {rows[numRows] = i; vals[numRows++] = aa[i];} 262 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 263 aa += a->lda; 264 } 265 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 266 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 267 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 268 269 if (reuse == MAT_INPLACE_MATRIX) { 270 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 271 } else if (reuse != MAT_REUSE_MATRIX) *newmat = B; 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 276 { 277 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 278 const PetscScalar *xv; 279 PetscScalar *yv; 280 PetscBLASInt N,m,ldax = 0,lday = 0,one = 1; 281 PetscErrorCode ierr; 282 283 PetscFunctionBegin; 284 ierr = MatDenseGetArrayRead(X,&xv);CHKERRQ(ierr); 285 ierr = MatDenseGetArray(Y,&yv);CHKERRQ(ierr); 286 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 287 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 288 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 289 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 290 if (ldax>m || lday>m) { 291 PetscInt j; 292 293 for (j=0; j<X->cmap->n; j++) { 294 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&alpha,xv+j*ldax,&one,yv+j*lday,&one)); 295 } 296 } else { 297 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&alpha,xv,&one,yv,&one)); 298 } 299 ierr = MatDenseRestoreArrayRead(X,&xv);CHKERRQ(ierr); 300 ierr = MatDenseRestoreArray(Y,&yv);CHKERRQ(ierr); 301 ierr = PetscLogFlops(PetscMax(2.0*N-1,0));CHKERRQ(ierr); 302 PetscFunctionReturn(0); 303 } 304 305 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 306 { 307 PetscLogDouble N = A->rmap->n*A->cmap->n; 308 309 PetscFunctionBegin; 310 info->block_size = 1.0; 311 info->nz_allocated = N; 312 info->nz_used = N; 313 info->nz_unneeded = 0; 314 info->assemblies = A->num_ass; 315 info->mallocs = 0; 316 info->memory = ((PetscObject)A)->mem; 317 info->fill_ratio_given = 0; 318 info->fill_ratio_needed = 0; 319 info->factor_mallocs = 0; 320 PetscFunctionReturn(0); 321 } 322 323 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 324 { 325 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 326 PetscScalar *v; 327 PetscErrorCode ierr; 328 PetscBLASInt one = 1,j,nz,lda = 0; 329 330 PetscFunctionBegin; 331 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 332 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 333 if (lda>A->rmap->n) { 334 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 335 for (j=0; j<A->cmap->n; j++) { 336 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v+j*lda,&one)); 337 } 338 } else { 339 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 340 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&alpha,v,&one)); 341 } 342 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 343 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 347 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 348 { 349 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 350 PetscInt i,j,m = A->rmap->n,N = a->lda; 351 const PetscScalar *v; 352 PetscErrorCode ierr; 353 354 PetscFunctionBegin; 355 *fl = PETSC_FALSE; 356 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 357 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 358 for (i=0; i<m; i++) { 359 for (j=i; j<m; j++) { 360 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) { 361 goto restore; 362 } 363 } 364 } 365 *fl = PETSC_TRUE; 366 restore: 367 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 368 PetscFunctionReturn(0); 369 } 370 371 static PetscErrorCode MatIsSymmetric_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 372 { 373 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 374 PetscInt i,j,m = A->rmap->n,N = a->lda; 375 const PetscScalar *v; 376 PetscErrorCode ierr; 377 378 PetscFunctionBegin; 379 *fl = PETSC_FALSE; 380 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 381 ierr = MatDenseGetArrayRead(A,&v);CHKERRQ(ierr); 382 for (i=0; i<m; i++) { 383 for (j=i; j<m; j++) { 384 if (PetscAbsScalar(v[i+j*N] - v[j+i*N]) > rtol) { 385 goto restore; 386 } 387 } 388 } 389 *fl = PETSC_TRUE; 390 restore: 391 ierr = MatDenseRestoreArrayRead(A,&v);CHKERRQ(ierr); 392 PetscFunctionReturn(0); 393 } 394 395 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 396 { 397 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 398 PetscErrorCode ierr; 399 PetscInt lda = (PetscInt)mat->lda,j,m,nlda = lda; 400 401 PetscFunctionBegin; 402 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 403 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 404 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { /* propagate LDA */ 405 ierr = MatDenseSetLDA(newi,lda);CHKERRQ(ierr); 406 } 407 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 408 if (cpvalues == MAT_COPY_VALUES) { 409 const PetscScalar *av; 410 PetscScalar *v; 411 412 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 413 ierr = MatDenseGetArray(newi,&v);CHKERRQ(ierr); 414 ierr = MatDenseGetLDA(newi,&nlda);CHKERRQ(ierr); 415 m = A->rmap->n; 416 if (lda>m || nlda>m) { 417 for (j=0; j<A->cmap->n; j++) { 418 ierr = PetscArraycpy(v+j*nlda,av+j*lda,m);CHKERRQ(ierr); 419 } 420 } else { 421 ierr = PetscArraycpy(v,av,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 422 } 423 ierr = MatDenseRestoreArray(newi,&v);CHKERRQ(ierr); 424 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 425 } 426 PetscFunctionReturn(0); 427 } 428 429 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 430 { 431 PetscErrorCode ierr; 432 433 PetscFunctionBegin; 434 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 435 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 436 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 437 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 static PetscErrorCode MatSolve_SeqDense_Internal_LU(Mat A, PetscScalar *x, PetscBLASInt 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 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 452 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 453 PetscFunctionReturn(0); 454 } 455 456 static PetscErrorCode MatConjugate_SeqDense(Mat); 457 458 static PetscErrorCode MatSolve_SeqDense_Internal_Cholesky(Mat A, PetscScalar *x, PetscBLASInt 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 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 471 if (PetscDefined(USE_COMPLEX) && T) {ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);} 472 #if defined(PETSC_USE_COMPLEX) 473 } else if (A->hermitian) { 474 if (T) {ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);} 475 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 476 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 477 ierr = PetscFPTrapPop();CHKERRQ(ierr); 478 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 479 if (T) {ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);} 480 #endif 481 } else { /* symmetric case */ 482 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 483 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 484 ierr = PetscFPTrapPop();CHKERRQ(ierr); 485 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 486 } 487 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490 491 static PetscErrorCode MatSolve_SeqDense_Internal_QR(Mat A, PetscScalar *x, PetscBLASInt 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 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform"); 508 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 509 PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U", "N", "N", &mat->rank,&nrhs,mat->v,&mat->lda,x,&ldx,&info)); 510 ierr = PetscFPTrapPop();CHKERRQ(ierr); 511 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve"); 512 for (PetscInt j = 0; j < nrhs; j++) { 513 for (PetscInt i = mat->rank; i < k; i++) { 514 x[j*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 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"TRTRS - Bad triangular solve"); 533 if (PetscDefined(USE_COMPLEX)) {ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);} 534 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 535 PetscStackCallBLAS("LAPACKormqr",LAPACKormqr_("L", "N", &m,&nrhs,&mat->rank,mat->v,&mat->lda,mat->tau,x,&ldx,mat->fwork,&mat->lfwork,&info)); 536 ierr = PetscFPTrapPop();CHKERRQ(ierr); 537 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"ORMQR - Bad orthogonal transform"); 538 if (PetscDefined(USE_COMPLEX)) {ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr);} 539 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"QR factored matrix cannot be used for transpose solve"); 540 ierr = PetscLogFlops(nrhs*(4.0*m*mat->rank - PetscSqr(mat->rank)));CHKERRQ(ierr); 541 PetscFunctionReturn(0); 542 } 543 544 static PetscErrorCode MatSolve_SeqDense_SetUp(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k) 545 { 546 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 547 PetscScalar *y; 548 PetscBLASInt m=0, k=0; 549 PetscErrorCode ierr; 550 551 PetscFunctionBegin; 552 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 553 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 554 if (k < m) { 555 ierr = VecCopy(xx, mat->qrrhs);CHKERRQ(ierr); 556 ierr = VecGetArray(mat->qrrhs,&y);CHKERRQ(ierr); 557 } else { 558 ierr = VecCopy(xx, yy);CHKERRQ(ierr); 559 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 560 } 561 *_y = y; 562 *_k = k; 563 *_m = m; 564 PetscFunctionReturn(0); 565 } 566 567 static PetscErrorCode MatSolve_SeqDense_TearDown(Mat A, Vec xx, Vec yy, PetscScalar **_y, PetscBLASInt *_m, PetscBLASInt *_k) 568 { 569 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 570 PetscScalar *y = NULL; 571 PetscBLASInt m, k; 572 PetscErrorCode ierr; 573 574 PetscFunctionBegin; 575 y = *_y; 576 *_y = NULL; 577 k = *_k; 578 m = *_m; 579 if (k < m) { 580 PetscScalar *yv; 581 ierr = VecGetArray(yy,&yv);CHKERRQ(ierr); 582 ierr = PetscArraycpy(yv, y, k);CHKERRQ(ierr); 583 ierr = VecRestoreArray(yy,&yv);CHKERRQ(ierr); 584 ierr = VecRestoreArray(mat->qrrhs, &y);CHKERRQ(ierr); 585 } else { 586 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 587 } 588 PetscFunctionReturn(0); 589 } 590 591 static PetscErrorCode MatSolve_SeqDense_LU(Mat A, Vec xx, Vec yy) 592 { 593 PetscScalar *y = NULL; 594 PetscBLASInt m = 0, k = 0; 595 PetscErrorCode ierr; 596 597 PetscFunctionBegin; 598 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 599 ierr = MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_FALSE);CHKERRQ(ierr); 600 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 601 PetscFunctionReturn(0); 602 } 603 604 static PetscErrorCode MatSolveTranspose_SeqDense_LU(Mat A, Vec xx, Vec yy) 605 { 606 PetscScalar *y = NULL; 607 PetscBLASInt m = 0, k = 0; 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 612 ierr = MatSolve_SeqDense_Internal_LU(A, y, m, m, 1, k, PETSC_TRUE);CHKERRQ(ierr); 613 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 614 PetscFunctionReturn(0); 615 } 616 617 static PetscErrorCode MatSolve_SeqDense_Cholesky(Mat A, Vec xx, Vec yy) 618 { 619 PetscScalar *y = NULL; 620 PetscBLASInt m = 0, k = 0; 621 PetscErrorCode ierr; 622 623 PetscFunctionBegin; 624 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 625 ierr = MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_FALSE);CHKERRQ(ierr); 626 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 627 PetscFunctionReturn(0); 628 } 629 630 static PetscErrorCode MatSolveTranspose_SeqDense_Cholesky(Mat A, Vec xx, Vec yy) 631 { 632 PetscScalar *y = NULL; 633 PetscBLASInt m = 0, k = 0; 634 PetscErrorCode ierr; 635 636 PetscFunctionBegin; 637 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 638 ierr = MatSolve_SeqDense_Internal_Cholesky(A, y, m, m, 1, k, PETSC_TRUE);CHKERRQ(ierr); 639 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 static PetscErrorCode MatSolve_SeqDense_QR(Mat A, Vec xx, Vec yy) 644 { 645 PetscScalar *y = NULL; 646 PetscBLASInt m = 0, k = 0; 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 651 ierr = MatSolve_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);CHKERRQ(ierr); 652 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 static PetscErrorCode MatSolveTranspose_SeqDense_QR(Mat A, Vec xx, Vec yy) 657 { 658 PetscScalar *y = NULL; 659 PetscBLASInt m = 0, k = 0; 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 ierr = MatSolve_SeqDense_SetUp(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 664 ierr = MatSolveTranspose_SeqDense_Internal_QR(A, y, PetscMax(m,k), m, 1, k);CHKERRQ(ierr); 665 ierr = MatSolve_SeqDense_TearDown(A, xx, yy, &y, &m, &k);CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668 669 static PetscErrorCode MatMatSolve_SeqDense_SetUp(Mat A, Mat B, Mat X, PetscScalar **_y, PetscBLASInt *_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; 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 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 852 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 853 854 A->ops->solve = MatSolve_SeqDense_LU; 855 A->ops->matsolve = MatMatSolve_SeqDense_LU; 856 A->ops->solvetranspose = MatSolveTranspose_SeqDense_LU; 857 A->ops->matsolvetranspose = MatMatSolveTranspose_SeqDense_LU; 858 A->factortype = MAT_FACTOR_LU; 859 860 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 861 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 862 863 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 864 PetscFunctionReturn(0); 865 } 866 867 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 868 { 869 MatFactorInfo info; 870 PetscErrorCode ierr; 871 872 PetscFunctionBegin; 873 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 874 ierr = (*fact->ops->lufactor)(fact,NULL,NULL,&info);CHKERRQ(ierr); 875 PetscFunctionReturn(0); 876 } 877 878 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 879 { 880 PetscFunctionBegin; 881 fact->preallocated = PETSC_TRUE; 882 fact->assembled = PETSC_TRUE; 883 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 884 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 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(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 if (info) SETERRQ(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 if (A->offloadmask == PETSC_OFFLOAD_GPU) SETERRQ(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 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D 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 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 1281 for (i=0; i<m; i++) { 1282 if (indexm[i] < 0) {idx++; continue;} 1283 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 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 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 1291 for (i=0; i<m; i++) { 1292 if (indexm[i] < 0) {idx++; continue;} 1293 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 1294 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1295 } 1296 } 1297 } 1298 } else { 1299 if (addv == INSERT_VALUES) { 1300 for (i=0; i<m; i++) { 1301 if (indexm[i] < 0) { idx += n; continue;} 1302 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 1303 for (j=0; j<n; j++) { 1304 if (indexn[j] < 0) { idx++; continue;} 1305 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 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 if (PetscUnlikelyDebug(indexm[i] >= A->rmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 1313 for (j=0; j<n; j++) { 1314 if (indexn[j] < 0) { idx++; continue;} 1315 if (PetscUnlikelyDebug(indexn[j] >= A->cmap->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 1316 av[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1317 } 1318 } 1319 } 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 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 1346 for (j=0; j<n; j++) { 1347 if (indexn[j] < 0) {v++; continue;} 1348 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 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 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 1422 M = header[1]; N = header[2]; 1423 if (M < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%D) in file is negative",M); 1424 if (N < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%D) in file is negative",N); 1425 nz = header[3]; 1426 if (nz != MATRIX_BINARY_FORMAT_DENSE && nz < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Unknown matrix format %D in file",nz); 1427 } else { 1428 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 1429 if (M < 0 || N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Matrix binary file header was skipped, thus the user must specify the global sizes of input matrix"); 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 if (M != rows || N != cols) SETERRQ4(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 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 SETERRQ2(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 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 %D:",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," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1530 } else if (PetscRealPart(*v)) { 1531 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1532 } 1533 #else 1534 if (*v) { 1535 ierr = PetscViewerASCIIPrintf(viewer," (%D, %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 = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1555 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\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 1681 if (iascii) { 1682 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1683 } else if (isbinary) { 1684 ierr = MatView_Dense_Binary(A,viewer);CHKERRQ(ierr); 1685 } else if (isdraw) { 1686 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1687 } 1688 PetscFunctionReturn(0); 1689 } 1690 1691 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar *array) 1692 { 1693 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1694 1695 PetscFunctionBegin; 1696 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1697 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1698 if (a->unplacedarray) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreArray() first"); 1699 a->unplacedarray = a->v; 1700 a->unplaced_user_alloc = a->user_alloc; 1701 a->v = (PetscScalar*) array; 1702 a->user_alloc = PETSC_TRUE; 1703 #if defined(PETSC_HAVE_CUDA) 1704 A->offloadmask = PETSC_OFFLOAD_CPU; 1705 #endif 1706 PetscFunctionReturn(0); 1707 } 1708 1709 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1710 { 1711 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1712 1713 PetscFunctionBegin; 1714 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1715 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1716 a->v = a->unplacedarray; 1717 a->user_alloc = a->unplaced_user_alloc; 1718 a->unplacedarray = NULL; 1719 #if defined(PETSC_HAVE_CUDA) 1720 A->offloadmask = PETSC_OFFLOAD_CPU; 1721 #endif 1722 PetscFunctionReturn(0); 1723 } 1724 1725 static PetscErrorCode MatDenseReplaceArray_SeqDense(Mat A,const PetscScalar *array) 1726 { 1727 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1728 PetscErrorCode ierr; 1729 1730 PetscFunctionBegin; 1731 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1732 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1733 if (!a->user_alloc) { ierr = PetscFree(a->v);CHKERRQ(ierr); } 1734 a->v = (PetscScalar*) array; 1735 a->user_alloc = PETSC_FALSE; 1736 #if defined(PETSC_HAVE_CUDA) 1737 A->offloadmask = PETSC_OFFLOAD_CPU; 1738 #endif 1739 PetscFunctionReturn(0); 1740 } 1741 1742 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1743 { 1744 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1745 PetscErrorCode ierr; 1746 1747 PetscFunctionBegin; 1748 #if defined(PETSC_USE_LOG) 1749 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1750 #endif 1751 ierr = VecDestroy(&(l->qrrhs));CHKERRQ(ierr); 1752 ierr = PetscFree(l->tau);CHKERRQ(ierr); 1753 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1754 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1755 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1756 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1757 if (!l->unplaced_user_alloc) {ierr = PetscFree(l->unplacedarray);CHKERRQ(ierr);} 1758 if (l->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 1759 if (l->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 1760 ierr = VecDestroy(&l->cvec);CHKERRQ(ierr); 1761 ierr = MatDestroy(&l->cmat);CHKERRQ(ierr); 1762 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1763 1764 ierr = PetscObjectChangeTypeName((PetscObject)mat,NULL);CHKERRQ(ierr); 1765 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatQRFactor_C",NULL);CHKERRQ(ierr); 1766 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1767 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseSetLDA_C",NULL);CHKERRQ(ierr); 1768 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1769 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1770 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1771 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1772 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseReplaceArray_C",NULL);CHKERRQ(ierr); 1773 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayRead_C",NULL);CHKERRQ(ierr); 1774 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayRead_C",NULL);CHKERRQ(ierr); 1775 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArrayWrite_C",NULL);CHKERRQ(ierr); 1776 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArrayWrite_C",NULL);CHKERRQ(ierr); 1777 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1778 #if defined(PETSC_HAVE_ELEMENTAL) 1779 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1780 #endif 1781 #if defined(PETSC_HAVE_SCALAPACK) 1782 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_scalapack_C",NULL);CHKERRQ(ierr); 1783 #endif 1784 #if defined(PETSC_HAVE_CUDA) 1785 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqdensecuda_C",NULL);CHKERRQ(ierr); 1786 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",NULL);CHKERRQ(ierr); 1787 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdensecuda_seqdense_C",NULL);CHKERRQ(ierr); 1788 #endif 1789 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1790 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1791 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqdense_seqdense_C",NULL);CHKERRQ(ierr); 1792 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqbaij_seqdense_C",NULL);CHKERRQ(ierr); 1793 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatProductSetFromOptions_seqsbaij_seqdense_C",NULL);CHKERRQ(ierr); 1794 1795 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1796 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1797 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVec_C",NULL);CHKERRQ(ierr); 1798 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVec_C",NULL);CHKERRQ(ierr); 1799 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecRead_C",NULL);CHKERRQ(ierr); 1800 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecRead_C",NULL);CHKERRQ(ierr); 1801 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumnVecWrite_C",NULL);CHKERRQ(ierr); 1802 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumnVecWrite_C",NULL);CHKERRQ(ierr); 1803 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetSubMatrix_C",NULL);CHKERRQ(ierr); 1804 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreSubMatrix_C",NULL);CHKERRQ(ierr); 1805 PetscFunctionReturn(0); 1806 } 1807 1808 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1809 { 1810 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1811 PetscErrorCode ierr; 1812 PetscInt k,j,m = A->rmap->n, M = mat->lda, n = A->cmap->n; 1813 PetscScalar *v,tmp; 1814 1815 PetscFunctionBegin; 1816 if (reuse == MAT_INPLACE_MATRIX) { 1817 if (m == n) { /* in place transpose */ 1818 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1819 for (j=0; j<m; j++) { 1820 for (k=0; k<j; k++) { 1821 tmp = v[j + k*M]; 1822 v[j + k*M] = v[k + j*M]; 1823 v[k + j*M] = tmp; 1824 } 1825 } 1826 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1827 } else { /* reuse memory, temporary allocates new memory */ 1828 PetscScalar *v2; 1829 PetscLayout tmplayout; 1830 1831 ierr = PetscMalloc1((size_t)m*n,&v2);CHKERRQ(ierr); 1832 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 1833 for (j=0; j<n; j++) { 1834 for (k=0; k<m; k++) v2[j + (size_t)k*n] = v[k + (size_t)j*M]; 1835 } 1836 ierr = PetscArraycpy(v,v2,(size_t)m*n);CHKERRQ(ierr); 1837 ierr = PetscFree(v2);CHKERRQ(ierr); 1838 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 1839 /* cleanup size dependent quantities */ 1840 ierr = VecDestroy(&mat->cvec);CHKERRQ(ierr); 1841 ierr = MatDestroy(&mat->cmat);CHKERRQ(ierr); 1842 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 1843 ierr = PetscFree(mat->fwork);CHKERRQ(ierr); 1844 ierr = MatDestroy(&mat->ptapwork);CHKERRQ(ierr); 1845 /* swap row/col layouts */ 1846 mat->lda = n; 1847 tmplayout = A->rmap; 1848 A->rmap = A->cmap; 1849 A->cmap = tmplayout; 1850 } 1851 } else { /* out-of-place transpose */ 1852 Mat tmat; 1853 Mat_SeqDense *tmatd; 1854 PetscScalar *v2; 1855 PetscInt M2; 1856 1857 if (reuse == MAT_INITIAL_MATRIX) { 1858 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1859 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1860 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1861 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1862 } else tmat = *matout; 1863 1864 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1865 ierr = MatDenseGetArray(tmat,&v2);CHKERRQ(ierr); 1866 tmatd = (Mat_SeqDense*)tmat->data; 1867 M2 = tmatd->lda; 1868 for (j=0; j<n; j++) { 1869 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1870 } 1871 ierr = MatDenseRestoreArray(tmat,&v2);CHKERRQ(ierr); 1872 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&v);CHKERRQ(ierr); 1873 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1874 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1875 *matout = tmat; 1876 } 1877 PetscFunctionReturn(0); 1878 } 1879 1880 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1881 { 1882 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1883 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1884 PetscInt i; 1885 const PetscScalar *v1,*v2; 1886 PetscErrorCode ierr; 1887 1888 PetscFunctionBegin; 1889 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1890 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1891 ierr = MatDenseGetArrayRead(A1,&v1);CHKERRQ(ierr); 1892 ierr = MatDenseGetArrayRead(A2,&v2);CHKERRQ(ierr); 1893 for (i=0; i<A1->cmap->n; i++) { 1894 ierr = PetscArraycmp(v1,v2,A1->rmap->n,flg);CHKERRQ(ierr); 1895 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 1896 v1 += mat1->lda; 1897 v2 += mat2->lda; 1898 } 1899 ierr = MatDenseRestoreArrayRead(A1,&v1);CHKERRQ(ierr); 1900 ierr = MatDenseRestoreArrayRead(A2,&v2);CHKERRQ(ierr); 1901 *flg = PETSC_TRUE; 1902 PetscFunctionReturn(0); 1903 } 1904 1905 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1906 { 1907 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1908 PetscInt i,n,len; 1909 PetscScalar *x; 1910 const PetscScalar *vv; 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1915 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1916 len = PetscMin(A->rmap->n,A->cmap->n); 1917 ierr = MatDenseGetArrayRead(A,&vv);CHKERRQ(ierr); 1918 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1919 for (i=0; i<len; i++) { 1920 x[i] = vv[i*mat->lda + i]; 1921 } 1922 ierr = MatDenseRestoreArrayRead(A,&vv);CHKERRQ(ierr); 1923 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1924 PetscFunctionReturn(0); 1925 } 1926 1927 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1928 { 1929 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1930 const PetscScalar *l,*r; 1931 PetscScalar x,*v,*vv; 1932 PetscErrorCode ierr; 1933 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1934 1935 PetscFunctionBegin; 1936 ierr = MatDenseGetArray(A,&vv);CHKERRQ(ierr); 1937 if (ll) { 1938 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1939 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1940 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1941 for (i=0; i<m; i++) { 1942 x = l[i]; 1943 v = vv + i; 1944 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1945 } 1946 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1947 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1948 } 1949 if (rr) { 1950 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1951 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1952 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1953 for (i=0; i<n; i++) { 1954 x = r[i]; 1955 v = vv + i*mat->lda; 1956 for (j=0; j<m; j++) (*v++) *= x; 1957 } 1958 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1959 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1960 } 1961 ierr = MatDenseRestoreArray(A,&vv);CHKERRQ(ierr); 1962 PetscFunctionReturn(0); 1963 } 1964 1965 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1966 { 1967 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1968 PetscScalar *v,*vv; 1969 PetscReal sum = 0.0; 1970 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1971 PetscErrorCode ierr; 1972 1973 PetscFunctionBegin; 1974 ierr = MatDenseGetArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 1975 v = vv; 1976 if (type == NORM_FROBENIUS) { 1977 if (lda>m) { 1978 for (j=0; j<A->cmap->n; j++) { 1979 v = vv+j*lda; 1980 for (i=0; i<m; i++) { 1981 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1982 } 1983 } 1984 } else { 1985 #if defined(PETSC_USE_REAL___FP16) 1986 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1987 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&cnt,v,&one)); 1988 } 1989 #else 1990 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1991 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1992 } 1993 } 1994 *nrm = PetscSqrtReal(sum); 1995 #endif 1996 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1997 } else if (type == NORM_1) { 1998 *nrm = 0.0; 1999 for (j=0; j<A->cmap->n; j++) { 2000 v = vv + j*mat->lda; 2001 sum = 0.0; 2002 for (i=0; i<A->rmap->n; i++) { 2003 sum += PetscAbsScalar(*v); v++; 2004 } 2005 if (sum > *nrm) *nrm = sum; 2006 } 2007 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 2008 } else if (type == NORM_INFINITY) { 2009 *nrm = 0.0; 2010 for (j=0; j<A->rmap->n; j++) { 2011 v = vv + j; 2012 sum = 0.0; 2013 for (i=0; i<A->cmap->n; i++) { 2014 sum += PetscAbsScalar(*v); v += mat->lda; 2015 } 2016 if (sum > *nrm) *nrm = sum; 2017 } 2018 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 2019 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 2020 ierr = MatDenseRestoreArrayRead(A,(const PetscScalar**)&vv);CHKERRQ(ierr); 2021 PetscFunctionReturn(0); 2022 } 2023 2024 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 2025 { 2026 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 2027 PetscErrorCode ierr; 2028 2029 PetscFunctionBegin; 2030 switch (op) { 2031 case MAT_ROW_ORIENTED: 2032 aij->roworiented = flg; 2033 break; 2034 case MAT_NEW_NONZERO_LOCATIONS: 2035 case MAT_NEW_NONZERO_LOCATION_ERR: 2036 case MAT_NEW_NONZERO_ALLOCATION_ERR: 2037 case MAT_FORCE_DIAGONAL_ENTRIES: 2038 case MAT_KEEP_NONZERO_PATTERN: 2039 case MAT_IGNORE_OFF_PROC_ENTRIES: 2040 case MAT_USE_HASH_TABLE: 2041 case MAT_IGNORE_ZERO_ENTRIES: 2042 case MAT_IGNORE_LOWER_TRIANGULAR: 2043 case MAT_SORTED_FULL: 2044 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 2045 break; 2046 case MAT_SPD: 2047 case MAT_SYMMETRIC: 2048 case MAT_STRUCTURALLY_SYMMETRIC: 2049 case MAT_HERMITIAN: 2050 case MAT_SYMMETRY_ETERNAL: 2051 /* These options are handled directly by MatSetOption() */ 2052 break; 2053 default: 2054 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 2055 } 2056 PetscFunctionReturn(0); 2057 } 2058 2059 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 2060 { 2061 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2062 PetscErrorCode ierr; 2063 PetscInt lda=l->lda,m=A->rmap->n,n=A->cmap->n,j; 2064 PetscScalar *v; 2065 2066 PetscFunctionBegin; 2067 ierr = MatDenseGetArrayWrite(A,&v);CHKERRQ(ierr); 2068 if (lda>m) { 2069 for (j=0; j<n; j++) { 2070 ierr = PetscArrayzero(v+j*lda,m);CHKERRQ(ierr); 2071 } 2072 } else { 2073 ierr = PetscArrayzero(v,PetscInt64Mult(m,n));CHKERRQ(ierr); 2074 } 2075 ierr = MatDenseRestoreArrayWrite(A,&v);CHKERRQ(ierr); 2076 PetscFunctionReturn(0); 2077 } 2078 2079 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2080 { 2081 PetscErrorCode ierr; 2082 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 2083 PetscInt m = l->lda, n = A->cmap->n, i,j; 2084 PetscScalar *slot,*bb,*v; 2085 const PetscScalar *xx; 2086 2087 PetscFunctionBegin; 2088 if (PetscDefined(USE_DEBUG)) { 2089 for (i=0; i<N; i++) { 2090 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 2091 if (rows[i] >= A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested to be zeroed greater than or equal number of rows %D",rows[i],A->rmap->n); 2092 } 2093 } 2094 if (!N) PetscFunctionReturn(0); 2095 2096 /* fix right hand side if needed */ 2097 if (x && b) { 2098 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2099 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 2100 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 2101 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2102 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 2103 } 2104 2105 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2106 for (i=0; i<N; i++) { 2107 slot = v + rows[i]; 2108 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 2109 } 2110 if (diag != 0.0) { 2111 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 2112 for (i=0; i<N; i++) { 2113 slot = v + (m+1)*rows[i]; 2114 *slot = diag; 2115 } 2116 } 2117 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2118 PetscFunctionReturn(0); 2119 } 2120 2121 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 2122 { 2123 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2124 2125 PetscFunctionBegin; 2126 *lda = mat->lda; 2127 PetscFunctionReturn(0); 2128 } 2129 2130 PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar **array) 2131 { 2132 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2133 2134 PetscFunctionBegin; 2135 if (mat->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 2136 *array = mat->v; 2137 PetscFunctionReturn(0); 2138 } 2139 2140 PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar **array) 2141 { 2142 PetscFunctionBegin; 2143 *array = NULL; 2144 PetscFunctionReturn(0); 2145 } 2146 2147 /*@C 2148 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 2149 2150 Not collective 2151 2152 Input Parameter: 2153 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2154 2155 Output Parameter: 2156 . lda - the leading dimension 2157 2158 Level: intermediate 2159 2160 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseSetLDA() 2161 @*/ 2162 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 2163 { 2164 PetscErrorCode ierr; 2165 2166 PetscFunctionBegin; 2167 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2168 PetscValidPointer(lda,2); 2169 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 2170 PetscFunctionReturn(0); 2171 } 2172 2173 /*@C 2174 MatDenseSetLDA - Sets the leading dimension of the array used by the dense matrix 2175 2176 Not collective 2177 2178 Input Parameter: 2179 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2180 - lda - the leading dimension 2181 2182 Level: intermediate 2183 2184 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetLDA() 2185 @*/ 2186 PetscErrorCode MatDenseSetLDA(Mat A,PetscInt lda) 2187 { 2188 PetscErrorCode ierr; 2189 2190 PetscFunctionBegin; 2191 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2192 ierr = PetscTryMethod(A,"MatDenseSetLDA_C",(Mat,PetscInt),(A,lda));CHKERRQ(ierr); 2193 PetscFunctionReturn(0); 2194 } 2195 2196 /*@C 2197 MatDenseGetArray - gives read-write access to the array where the data for a dense matrix is stored 2198 2199 Logically Collective on Mat 2200 2201 Input Parameter: 2202 . mat - a dense matrix 2203 2204 Output Parameter: 2205 . array - pointer to the data 2206 2207 Level: intermediate 2208 2209 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 2210 @*/ 2211 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 2212 { 2213 PetscErrorCode ierr; 2214 2215 PetscFunctionBegin; 2216 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2217 PetscValidPointer(array,2); 2218 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 2219 PetscFunctionReturn(0); 2220 } 2221 2222 /*@C 2223 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 2224 2225 Logically Collective on Mat 2226 2227 Input Parameters: 2228 + mat - a dense matrix 2229 - array - pointer to the data 2230 2231 Level: intermediate 2232 2233 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 2234 @*/ 2235 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 2236 { 2237 PetscErrorCode ierr; 2238 2239 PetscFunctionBegin; 2240 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2241 PetscValidPointer(array,2); 2242 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 2243 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 2244 #if defined(PETSC_HAVE_CUDA) 2245 A->offloadmask = PETSC_OFFLOAD_CPU; 2246 #endif 2247 PetscFunctionReturn(0); 2248 } 2249 2250 /*@C 2251 MatDenseGetArrayRead - gives read-only access to the array where the data for a dense matrix is stored 2252 2253 Not Collective 2254 2255 Input Parameter: 2256 . mat - a dense matrix 2257 2258 Output Parameter: 2259 . array - pointer to the data 2260 2261 Level: intermediate 2262 2263 .seealso: MatDenseRestoreArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 2264 @*/ 2265 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 2266 { 2267 PetscErrorCode ierr; 2268 2269 PetscFunctionBegin; 2270 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2271 PetscValidPointer(array,2); 2272 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 2273 PetscFunctionReturn(0); 2274 } 2275 2276 /*@C 2277 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayRead() 2278 2279 Not Collective 2280 2281 Input Parameters: 2282 + mat - a dense matrix 2283 - array - pointer to the data 2284 2285 Level: intermediate 2286 2287 .seealso: MatDenseGetArrayRead(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayWrite(), MatDenseRestoreArrayWrite() 2288 @*/ 2289 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 2290 { 2291 PetscErrorCode ierr; 2292 2293 PetscFunctionBegin; 2294 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2295 PetscValidPointer(array,2); 2296 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 2297 PetscFunctionReturn(0); 2298 } 2299 2300 /*@C 2301 MatDenseGetArrayWrite - gives write-only access to the array where the data for a dense matrix is stored 2302 2303 Not Collective 2304 2305 Input Parameter: 2306 . mat - a dense matrix 2307 2308 Output Parameter: 2309 . array - pointer to the data 2310 2311 Level: intermediate 2312 2313 .seealso: MatDenseRestoreArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 2314 @*/ 2315 PetscErrorCode MatDenseGetArrayWrite(Mat A,PetscScalar **array) 2316 { 2317 PetscErrorCode ierr; 2318 2319 PetscFunctionBegin; 2320 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2321 PetscValidPointer(array,2); 2322 ierr = PetscUseMethod(A,"MatDenseGetArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 2323 PetscFunctionReturn(0); 2324 } 2325 2326 /*@C 2327 MatDenseRestoreArrayWrite - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArrayWrite() 2328 2329 Not Collective 2330 2331 Input Parameters: 2332 + mat - a dense matrix 2333 - array - pointer to the data 2334 2335 Level: intermediate 2336 2337 .seealso: MatDenseGetArrayWrite(), MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 2338 @*/ 2339 PetscErrorCode MatDenseRestoreArrayWrite(Mat A,PetscScalar **array) 2340 { 2341 PetscErrorCode ierr; 2342 2343 PetscFunctionBegin; 2344 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 2345 PetscValidPointer(array,2); 2346 ierr = PetscUseMethod(A,"MatDenseRestoreArrayWrite_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 2347 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 2348 #if defined(PETSC_HAVE_CUDA) 2349 A->offloadmask = PETSC_OFFLOAD_CPU; 2350 #endif 2351 PetscFunctionReturn(0); 2352 } 2353 2354 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 2355 { 2356 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 2357 PetscErrorCode ierr; 2358 PetscInt i,j,nrows,ncols,ldb; 2359 const PetscInt *irow,*icol; 2360 PetscScalar *av,*bv,*v = mat->v; 2361 Mat newmat; 2362 2363 PetscFunctionBegin; 2364 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2365 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2366 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2367 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2368 2369 /* Check submatrixcall */ 2370 if (scall == MAT_REUSE_MATRIX) { 2371 PetscInt n_cols,n_rows; 2372 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2373 if (n_rows != nrows || n_cols != ncols) { 2374 /* resize the result matrix to match number of requested rows/columns */ 2375 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2376 } 2377 newmat = *B; 2378 } else { 2379 /* Create and fill new matrix */ 2380 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 2381 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 2382 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 2383 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 2384 } 2385 2386 /* Now extract the data pointers and do the copy,column at a time */ 2387 ierr = MatDenseGetArray(newmat,&bv);CHKERRQ(ierr); 2388 ierr = MatDenseGetLDA(newmat,&ldb);CHKERRQ(ierr); 2389 for (i=0; i<ncols; i++) { 2390 av = v + mat->lda*icol[i]; 2391 for (j=0; j<nrows; j++) bv[j] = av[irow[j]]; 2392 bv += ldb; 2393 } 2394 ierr = MatDenseRestoreArray(newmat,&bv);CHKERRQ(ierr); 2395 2396 /* Assemble the matrices so that the correct flags are set */ 2397 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2398 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2399 2400 /* Free work space */ 2401 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2402 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2403 *B = newmat; 2404 PetscFunctionReturn(0); 2405 } 2406 2407 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2408 { 2409 PetscErrorCode ierr; 2410 PetscInt i; 2411 2412 PetscFunctionBegin; 2413 if (scall == MAT_INITIAL_MATRIX) { 2414 ierr = PetscCalloc1(n,B);CHKERRQ(ierr); 2415 } 2416 2417 for (i=0; i<n; i++) { 2418 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 2419 } 2420 PetscFunctionReturn(0); 2421 } 2422 2423 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 2424 { 2425 PetscFunctionBegin; 2426 PetscFunctionReturn(0); 2427 } 2428 2429 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 2430 { 2431 PetscFunctionBegin; 2432 PetscFunctionReturn(0); 2433 } 2434 2435 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2436 { 2437 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2438 PetscErrorCode ierr; 2439 const PetscScalar *va; 2440 PetscScalar *vb; 2441 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2442 2443 PetscFunctionBegin; 2444 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2445 if (A->ops->copy != B->ops->copy) { 2446 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2447 PetscFunctionReturn(0); 2448 } 2449 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2450 ierr = MatDenseGetArrayRead(A,&va);CHKERRQ(ierr); 2451 ierr = MatDenseGetArray(B,&vb);CHKERRQ(ierr); 2452 if (lda1>m || lda2>m) { 2453 for (j=0; j<n; j++) { 2454 ierr = PetscArraycpy(vb+j*lda2,va+j*lda1,m);CHKERRQ(ierr); 2455 } 2456 } else { 2457 ierr = PetscArraycpy(vb,va,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2458 } 2459 ierr = MatDenseRestoreArray(B,&vb);CHKERRQ(ierr); 2460 ierr = MatDenseRestoreArrayRead(A,&va);CHKERRQ(ierr); 2461 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2462 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2463 PetscFunctionReturn(0); 2464 } 2465 2466 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2467 { 2468 PetscErrorCode ierr; 2469 2470 PetscFunctionBegin; 2471 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 2472 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 2473 if (!A->preallocated) { 2474 ierr = MatSeqDenseSetPreallocation(A,NULL);CHKERRQ(ierr); 2475 } 2476 PetscFunctionReturn(0); 2477 } 2478 2479 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2480 { 2481 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 2482 PetscInt i,nz = A->rmap->n*A->cmap->n; 2483 PetscInt min = PetscMin(A->rmap->n,A->cmap->n); 2484 PetscScalar *aa; 2485 PetscErrorCode ierr; 2486 2487 PetscFunctionBegin; 2488 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2489 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2490 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2491 if (mat->tau) for (i = 0; i < min; i++) mat->tau[i] = PetscConj(mat->tau[i]); 2492 PetscFunctionReturn(0); 2493 } 2494 2495 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2496 { 2497 PetscInt i,nz = A->rmap->n*A->cmap->n; 2498 PetscScalar *aa; 2499 PetscErrorCode ierr; 2500 2501 PetscFunctionBegin; 2502 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2503 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2504 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2505 PetscFunctionReturn(0); 2506 } 2507 2508 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2509 { 2510 PetscInt i,nz = A->rmap->n*A->cmap->n; 2511 PetscScalar *aa; 2512 PetscErrorCode ierr; 2513 2514 PetscFunctionBegin; 2515 ierr = MatDenseGetArray(A,&aa);CHKERRQ(ierr); 2516 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2517 ierr = MatDenseRestoreArray(A,&aa);CHKERRQ(ierr); 2518 PetscFunctionReturn(0); 2519 } 2520 2521 /* ----------------------------------------------------------------*/ 2522 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2523 { 2524 PetscErrorCode ierr; 2525 PetscInt m=A->rmap->n,n=B->cmap->n; 2526 PetscBool cisdense; 2527 2528 PetscFunctionBegin; 2529 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2530 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2531 if (!cisdense) { 2532 PetscBool flg; 2533 2534 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2535 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2536 } 2537 ierr = MatSetUp(C);CHKERRQ(ierr); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2542 { 2543 Mat_SeqDense *a=(Mat_SeqDense*)A->data,*b=(Mat_SeqDense*)B->data,*c=(Mat_SeqDense*)C->data; 2544 PetscBLASInt m,n,k; 2545 const PetscScalar *av,*bv; 2546 PetscScalar *cv; 2547 PetscScalar _DOne=1.0,_DZero=0.0; 2548 PetscErrorCode ierr; 2549 2550 PetscFunctionBegin; 2551 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2552 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2553 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2554 if (!m || !n || !k) PetscFunctionReturn(0); 2555 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2556 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2557 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2558 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2559 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2560 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2561 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2562 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2563 PetscFunctionReturn(0); 2564 } 2565 2566 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2567 { 2568 PetscErrorCode ierr; 2569 PetscInt m=A->rmap->n,n=B->rmap->n; 2570 PetscBool cisdense; 2571 2572 PetscFunctionBegin; 2573 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2574 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2575 if (!cisdense) { 2576 PetscBool flg; 2577 2578 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2579 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2580 } 2581 ierr = MatSetUp(C);CHKERRQ(ierr); 2582 PetscFunctionReturn(0); 2583 } 2584 2585 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2586 { 2587 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2588 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2589 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2590 const PetscScalar *av,*bv; 2591 PetscScalar *cv; 2592 PetscBLASInt m,n,k; 2593 PetscScalar _DOne=1.0,_DZero=0.0; 2594 PetscErrorCode ierr; 2595 2596 PetscFunctionBegin; 2597 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2598 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2599 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2600 if (!m || !n || !k) PetscFunctionReturn(0); 2601 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2602 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2603 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2604 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2605 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2606 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2607 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2608 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2609 PetscFunctionReturn(0); 2610 } 2611 2612 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat C) 2613 { 2614 PetscErrorCode ierr; 2615 PetscInt m=A->cmap->n,n=B->cmap->n; 2616 PetscBool cisdense; 2617 2618 PetscFunctionBegin; 2619 ierr = MatSetSizes(C,m,n,m,n);CHKERRQ(ierr); 2620 ierr = PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"");CHKERRQ(ierr); 2621 if (!cisdense) { 2622 PetscBool flg; 2623 2624 ierr = PetscObjectTypeCompare((PetscObject)B,((PetscObject)A)->type_name,&flg);CHKERRQ(ierr); 2625 ierr = MatSetType(C,flg ? ((PetscObject)A)->type_name : MATDENSE);CHKERRQ(ierr); 2626 } 2627 ierr = MatSetUp(C);CHKERRQ(ierr); 2628 PetscFunctionReturn(0); 2629 } 2630 2631 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2632 { 2633 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2634 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2635 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2636 const PetscScalar *av,*bv; 2637 PetscScalar *cv; 2638 PetscBLASInt m,n,k; 2639 PetscScalar _DOne=1.0,_DZero=0.0; 2640 PetscErrorCode ierr; 2641 2642 PetscFunctionBegin; 2643 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2644 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2645 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2646 if (!m || !n || !k) PetscFunctionReturn(0); 2647 ierr = MatDenseGetArrayRead(A,&av);CHKERRQ(ierr); 2648 ierr = MatDenseGetArrayRead(B,&bv);CHKERRQ(ierr); 2649 ierr = MatDenseGetArrayWrite(C,&cv);CHKERRQ(ierr); 2650 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,av,&a->lda,bv,&b->lda,&_DZero,cv,&c->lda)); 2651 ierr = MatDenseRestoreArrayRead(A,&av);CHKERRQ(ierr); 2652 ierr = MatDenseRestoreArrayRead(B,&bv);CHKERRQ(ierr); 2653 ierr = MatDenseRestoreArrayWrite(C,&cv);CHKERRQ(ierr); 2654 ierr = PetscLogFlops(1.0*m*n*k + 1.0*m*n*(k-1));CHKERRQ(ierr); 2655 PetscFunctionReturn(0); 2656 } 2657 2658 /* ----------------------------------------------- */ 2659 static PetscErrorCode MatProductSetFromOptions_SeqDense_AB(Mat C) 2660 { 2661 PetscFunctionBegin; 2662 C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqDense; 2663 C->ops->productsymbolic = MatProductSymbolic_AB; 2664 PetscFunctionReturn(0); 2665 } 2666 2667 static PetscErrorCode MatProductSetFromOptions_SeqDense_AtB(Mat C) 2668 { 2669 PetscFunctionBegin; 2670 C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_SeqDense_SeqDense; 2671 C->ops->productsymbolic = MatProductSymbolic_AtB; 2672 PetscFunctionReturn(0); 2673 } 2674 2675 static PetscErrorCode MatProductSetFromOptions_SeqDense_ABt(Mat C) 2676 { 2677 PetscFunctionBegin; 2678 C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqDense_SeqDense; 2679 C->ops->productsymbolic = MatProductSymbolic_ABt; 2680 PetscFunctionReturn(0); 2681 } 2682 2683 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense(Mat C) 2684 { 2685 PetscErrorCode ierr; 2686 Mat_Product *product = C->product; 2687 2688 PetscFunctionBegin; 2689 switch (product->type) { 2690 case MATPRODUCT_AB: 2691 ierr = MatProductSetFromOptions_SeqDense_AB(C);CHKERRQ(ierr); 2692 break; 2693 case MATPRODUCT_AtB: 2694 ierr = MatProductSetFromOptions_SeqDense_AtB(C);CHKERRQ(ierr); 2695 break; 2696 case MATPRODUCT_ABt: 2697 ierr = MatProductSetFromOptions_SeqDense_ABt(C);CHKERRQ(ierr); 2698 break; 2699 default: 2700 break; 2701 } 2702 PetscFunctionReturn(0); 2703 } 2704 /* ----------------------------------------------- */ 2705 2706 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2707 { 2708 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2709 PetscErrorCode ierr; 2710 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2711 PetscScalar *x; 2712 const PetscScalar *aa; 2713 2714 PetscFunctionBegin; 2715 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2716 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2717 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2718 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2719 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2720 for (i=0; i<m; i++) { 2721 x[i] = aa[i]; if (idx) idx[i] = 0; 2722 for (j=1; j<n; j++) { 2723 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2724 } 2725 } 2726 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2727 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2728 PetscFunctionReturn(0); 2729 } 2730 2731 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2732 { 2733 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2734 PetscErrorCode ierr; 2735 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2736 PetscScalar *x; 2737 PetscReal atmp; 2738 const PetscScalar *aa; 2739 2740 PetscFunctionBegin; 2741 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2742 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2743 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2744 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2745 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2746 for (i=0; i<m; i++) { 2747 x[i] = PetscAbsScalar(aa[i]); 2748 for (j=1; j<n; j++) { 2749 atmp = PetscAbsScalar(aa[i+a->lda*j]); 2750 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2751 } 2752 } 2753 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2754 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2755 PetscFunctionReturn(0); 2756 } 2757 2758 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2759 { 2760 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2761 PetscErrorCode ierr; 2762 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2763 PetscScalar *x; 2764 const PetscScalar *aa; 2765 2766 PetscFunctionBegin; 2767 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2768 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2769 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2770 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2771 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2772 for (i=0; i<m; i++) { 2773 x[i] = aa[i]; if (idx) idx[i] = 0; 2774 for (j=1; j<n; j++) { 2775 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+a->lda*j])) {x[i] = aa[i + a->lda*j]; if (idx) idx[i] = j;} 2776 } 2777 } 2778 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2779 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2780 PetscFunctionReturn(0); 2781 } 2782 2783 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2784 { 2785 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2786 PetscErrorCode ierr; 2787 PetscScalar *x; 2788 const PetscScalar *aa; 2789 2790 PetscFunctionBegin; 2791 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2792 ierr = MatDenseGetArrayRead(A,&aa);CHKERRQ(ierr); 2793 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2794 ierr = PetscArraycpy(x,aa+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2795 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2796 ierr = MatDenseRestoreArrayRead(A,&aa);CHKERRQ(ierr); 2797 PetscFunctionReturn(0); 2798 } 2799 2800 PETSC_INTERN PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2801 { 2802 PetscErrorCode ierr; 2803 PetscInt i,j,m,n; 2804 const PetscScalar *a; 2805 2806 PetscFunctionBegin; 2807 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2808 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2809 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2810 if (type == NORM_2) { 2811 for (i=0; i<n; i++) { 2812 for (j=0; j<m; j++) { 2813 norms[i] += PetscAbsScalar(a[j]*a[j]); 2814 } 2815 a += m; 2816 } 2817 } else if (type == NORM_1) { 2818 for (i=0; i<n; i++) { 2819 for (j=0; j<m; j++) { 2820 norms[i] += PetscAbsScalar(a[j]); 2821 } 2822 a += m; 2823 } 2824 } else if (type == NORM_INFINITY) { 2825 for (i=0; i<n; i++) { 2826 for (j=0; j<m; j++) { 2827 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2828 } 2829 a += m; 2830 } 2831 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2832 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2833 if (type == NORM_2) { 2834 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2835 } 2836 PetscFunctionReturn(0); 2837 } 2838 2839 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2840 { 2841 PetscErrorCode ierr; 2842 PetscScalar *a; 2843 PetscInt lda,m,n,i,j; 2844 2845 PetscFunctionBegin; 2846 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2847 ierr = MatDenseGetLDA(x,&lda);CHKERRQ(ierr); 2848 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2849 for (j=0; j<n; j++) { 2850 for (i=0; i<m; i++) { 2851 ierr = PetscRandomGetValue(rctx,a+j*lda+i);CHKERRQ(ierr); 2852 } 2853 } 2854 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2855 PetscFunctionReturn(0); 2856 } 2857 2858 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2859 { 2860 PetscFunctionBegin; 2861 *missing = PETSC_FALSE; 2862 PetscFunctionReturn(0); 2863 } 2864 2865 /* vals is not const */ 2866 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2867 { 2868 PetscErrorCode ierr; 2869 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2870 PetscScalar *v; 2871 2872 PetscFunctionBegin; 2873 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2874 ierr = MatDenseGetArray(A,&v);CHKERRQ(ierr); 2875 *vals = v+col*a->lda; 2876 ierr = MatDenseRestoreArray(A,&v);CHKERRQ(ierr); 2877 PetscFunctionReturn(0); 2878 } 2879 2880 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2881 { 2882 PetscFunctionBegin; 2883 *vals = NULL; /* user cannot accidently use the array later */ 2884 PetscFunctionReturn(0); 2885 } 2886 2887 /* -------------------------------------------------------------------*/ 2888 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2889 MatGetRow_SeqDense, 2890 MatRestoreRow_SeqDense, 2891 MatMult_SeqDense, 2892 /* 4*/ MatMultAdd_SeqDense, 2893 MatMultTranspose_SeqDense, 2894 MatMultTransposeAdd_SeqDense, 2895 NULL, 2896 NULL, 2897 NULL, 2898 /* 10*/ NULL, 2899 MatLUFactor_SeqDense, 2900 MatCholeskyFactor_SeqDense, 2901 MatSOR_SeqDense, 2902 MatTranspose_SeqDense, 2903 /* 15*/ MatGetInfo_SeqDense, 2904 MatEqual_SeqDense, 2905 MatGetDiagonal_SeqDense, 2906 MatDiagonalScale_SeqDense, 2907 MatNorm_SeqDense, 2908 /* 20*/ MatAssemblyBegin_SeqDense, 2909 MatAssemblyEnd_SeqDense, 2910 MatSetOption_SeqDense, 2911 MatZeroEntries_SeqDense, 2912 /* 24*/ MatZeroRows_SeqDense, 2913 NULL, 2914 NULL, 2915 NULL, 2916 NULL, 2917 /* 29*/ MatSetUp_SeqDense, 2918 NULL, 2919 NULL, 2920 NULL, 2921 NULL, 2922 /* 34*/ MatDuplicate_SeqDense, 2923 NULL, 2924 NULL, 2925 NULL, 2926 NULL, 2927 /* 39*/ MatAXPY_SeqDense, 2928 MatCreateSubMatrices_SeqDense, 2929 NULL, 2930 MatGetValues_SeqDense, 2931 MatCopy_SeqDense, 2932 /* 44*/ MatGetRowMax_SeqDense, 2933 MatScale_SeqDense, 2934 MatShift_Basic, 2935 NULL, 2936 MatZeroRowsColumns_SeqDense, 2937 /* 49*/ MatSetRandom_SeqDense, 2938 NULL, 2939 NULL, 2940 NULL, 2941 NULL, 2942 /* 54*/ NULL, 2943 NULL, 2944 NULL, 2945 NULL, 2946 NULL, 2947 /* 59*/ MatCreateSubMatrix_SeqDense, 2948 MatDestroy_SeqDense, 2949 MatView_SeqDense, 2950 NULL, 2951 NULL, 2952 /* 64*/ NULL, 2953 NULL, 2954 NULL, 2955 NULL, 2956 NULL, 2957 /* 69*/ MatGetRowMaxAbs_SeqDense, 2958 NULL, 2959 NULL, 2960 NULL, 2961 NULL, 2962 /* 74*/ NULL, 2963 NULL, 2964 NULL, 2965 NULL, 2966 NULL, 2967 /* 79*/ NULL, 2968 NULL, 2969 NULL, 2970 NULL, 2971 /* 83*/ MatLoad_SeqDense, 2972 MatIsSymmetric_SeqDense, 2973 MatIsHermitian_SeqDense, 2974 NULL, 2975 NULL, 2976 NULL, 2977 /* 89*/ NULL, 2978 NULL, 2979 MatMatMultNumeric_SeqDense_SeqDense, 2980 NULL, 2981 NULL, 2982 /* 94*/ NULL, 2983 NULL, 2984 NULL, 2985 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2986 NULL, 2987 /* 99*/ MatProductSetFromOptions_SeqDense, 2988 NULL, 2989 NULL, 2990 MatConjugate_SeqDense, 2991 NULL, 2992 /*104*/ NULL, 2993 MatRealPart_SeqDense, 2994 MatImaginaryPart_SeqDense, 2995 NULL, 2996 NULL, 2997 /*109*/ NULL, 2998 NULL, 2999 MatGetRowMin_SeqDense, 3000 MatGetColumnVector_SeqDense, 3001 MatMissingDiagonal_SeqDense, 3002 /*114*/ NULL, 3003 NULL, 3004 NULL, 3005 NULL, 3006 NULL, 3007 /*119*/ NULL, 3008 NULL, 3009 NULL, 3010 NULL, 3011 NULL, 3012 /*124*/ NULL, 3013 MatGetColumnNorms_SeqDense, 3014 NULL, 3015 NULL, 3016 NULL, 3017 /*129*/ NULL, 3018 NULL, 3019 NULL, 3020 MatTransposeMatMultNumeric_SeqDense_SeqDense, 3021 NULL, 3022 /*134*/ NULL, 3023 NULL, 3024 NULL, 3025 NULL, 3026 NULL, 3027 /*139*/ NULL, 3028 NULL, 3029 NULL, 3030 NULL, 3031 NULL, 3032 MatCreateMPIMatConcatenateSeqMat_SeqDense, 3033 /*145*/ NULL, 3034 NULL, 3035 NULL 3036 }; 3037 3038 /*@C 3039 MatCreateSeqDense - Creates a sequential dense matrix that 3040 is stored in column major order (the usual Fortran 77 manner). Many 3041 of the matrix operations use the BLAS and LAPACK routines. 3042 3043 Collective 3044 3045 Input Parameters: 3046 + comm - MPI communicator, set to PETSC_COMM_SELF 3047 . m - number of rows 3048 . n - number of columns 3049 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 3050 to control all matrix memory allocation. 3051 3052 Output Parameter: 3053 . A - the matrix 3054 3055 Notes: 3056 The data input variable is intended primarily for Fortran programmers 3057 who wish to allocate their own matrix memory space. Most users should 3058 set data=NULL. 3059 3060 Level: intermediate 3061 3062 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 3063 @*/ 3064 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 3065 { 3066 PetscErrorCode ierr; 3067 3068 PetscFunctionBegin; 3069 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3070 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3071 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 3072 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 3073 PetscFunctionReturn(0); 3074 } 3075 3076 /*@C 3077 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 3078 3079 Collective 3080 3081 Input Parameters: 3082 + B - the matrix 3083 - data - the array (or NULL) 3084 3085 Notes: 3086 The data input variable is intended primarily for Fortran programmers 3087 who wish to allocate their own matrix memory space. Most users should 3088 need not call this routine. 3089 3090 Level: intermediate 3091 3092 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatDenseSetLDA() 3093 3094 @*/ 3095 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 3096 { 3097 PetscErrorCode ierr; 3098 3099 PetscFunctionBegin; 3100 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3101 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 3102 PetscFunctionReturn(0); 3103 } 3104 3105 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 3106 { 3107 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3108 PetscErrorCode ierr; 3109 3110 PetscFunctionBegin; 3111 if (b->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3112 B->preallocated = PETSC_TRUE; 3113 3114 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3115 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3116 3117 if (b->lda <= 0) b->lda = B->rmap->n; 3118 3119 ierr = PetscIntMultError(b->lda,B->cmap->n,NULL);CHKERRQ(ierr); 3120 if (!data) { /* petsc-allocated storage */ 3121 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 3122 ierr = PetscCalloc1((size_t)b->lda*B->cmap->n,&b->v);CHKERRQ(ierr); 3123 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*B->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 3124 3125 b->user_alloc = PETSC_FALSE; 3126 } else { /* user-allocated storage */ 3127 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 3128 b->v = data; 3129 b->user_alloc = PETSC_TRUE; 3130 } 3131 B->assembled = PETSC_TRUE; 3132 PetscFunctionReturn(0); 3133 } 3134 3135 #if defined(PETSC_HAVE_ELEMENTAL) 3136 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 3137 { 3138 Mat mat_elemental; 3139 PetscErrorCode ierr; 3140 const PetscScalar *array; 3141 PetscScalar *v_colwise; 3142 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 3143 3144 PetscFunctionBegin; 3145 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 3146 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 3147 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 3148 k = 0; 3149 for (j=0; j<N; j++) { 3150 cols[j] = j; 3151 for (i=0; i<M; i++) { 3152 v_colwise[j*M+i] = array[k++]; 3153 } 3154 } 3155 for (i=0; i<M; i++) { 3156 rows[i] = i; 3157 } 3158 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 3159 3160 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 3161 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3162 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 3163 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 3164 3165 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 3166 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 3167 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3168 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3169 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 3170 3171 if (reuse == MAT_INPLACE_MATRIX) { 3172 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 3173 } else { 3174 *newmat = mat_elemental; 3175 } 3176 PetscFunctionReturn(0); 3177 } 3178 #endif 3179 3180 PetscErrorCode MatDenseSetLDA_SeqDense(Mat B,PetscInt lda) 3181 { 3182 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 3183 PetscBool data; 3184 3185 PetscFunctionBegin; 3186 data = (PetscBool)((B->rmap->n > 0 && B->cmap->n > 0) ? (b->v ? PETSC_TRUE : PETSC_FALSE) : PETSC_FALSE); 3187 if (!b->user_alloc && data && b->lda!=lda) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"LDA cannot be changed after allocation of internal storage"); 3188 if (lda < B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 3189 b->lda = lda; 3190 PetscFunctionReturn(0); 3191 } 3192 3193 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 3194 { 3195 PetscErrorCode ierr; 3196 PetscMPIInt size; 3197 3198 PetscFunctionBegin; 3199 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 3200 if (size == 1) { 3201 if (scall == MAT_INITIAL_MATRIX) { 3202 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 3203 } else { 3204 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 3205 } 3206 } else { 3207 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 3208 } 3209 PetscFunctionReturn(0); 3210 } 3211 3212 PetscErrorCode MatDenseGetColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3213 { 3214 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3215 PetscErrorCode ierr; 3216 3217 PetscFunctionBegin; 3218 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3219 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3220 if (!a->cvec) { 3221 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 3222 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 3223 } 3224 a->vecinuse = col + 1; 3225 ierr = MatDenseGetArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 3226 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 3227 *v = a->cvec; 3228 PetscFunctionReturn(0); 3229 } 3230 3231 PetscErrorCode MatDenseRestoreColumnVec_SeqDense(Mat A,PetscInt col,Vec *v) 3232 { 3233 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3234 PetscErrorCode ierr; 3235 3236 PetscFunctionBegin; 3237 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3238 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3239 a->vecinuse = 0; 3240 ierr = MatDenseRestoreArray(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 3241 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 3242 *v = NULL; 3243 PetscFunctionReturn(0); 3244 } 3245 3246 PetscErrorCode MatDenseGetColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3247 { 3248 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3249 PetscErrorCode ierr; 3250 3251 PetscFunctionBegin; 3252 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3253 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3254 if (!a->cvec) { 3255 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 3256 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 3257 } 3258 a->vecinuse = col + 1; 3259 ierr = MatDenseGetArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 3260 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 3261 ierr = VecLockReadPush(a->cvec);CHKERRQ(ierr); 3262 *v = a->cvec; 3263 PetscFunctionReturn(0); 3264 } 3265 3266 PetscErrorCode MatDenseRestoreColumnVecRead_SeqDense(Mat A,PetscInt col,Vec *v) 3267 { 3268 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3269 PetscErrorCode ierr; 3270 3271 PetscFunctionBegin; 3272 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3273 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3274 a->vecinuse = 0; 3275 ierr = MatDenseRestoreArrayRead(A,&a->ptrinuse);CHKERRQ(ierr); 3276 ierr = VecLockReadPop(a->cvec);CHKERRQ(ierr); 3277 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 3278 *v = NULL; 3279 PetscFunctionReturn(0); 3280 } 3281 3282 PetscErrorCode MatDenseGetColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3283 { 3284 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3285 PetscErrorCode ierr; 3286 3287 PetscFunctionBegin; 3288 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3289 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3290 if (!a->cvec) { 3291 ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)A),A->rmap->bs,A->rmap->n,NULL,&a->cvec);CHKERRQ(ierr); 3292 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cvec);CHKERRQ(ierr); 3293 } 3294 a->vecinuse = col + 1; 3295 ierr = MatDenseGetArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 3296 ierr = VecPlaceArray(a->cvec,a->ptrinuse + (size_t)col * (size_t)a->lda);CHKERRQ(ierr); 3297 *v = a->cvec; 3298 PetscFunctionReturn(0); 3299 } 3300 3301 PetscErrorCode MatDenseRestoreColumnVecWrite_SeqDense(Mat A,PetscInt col,Vec *v) 3302 { 3303 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3304 PetscErrorCode ierr; 3305 3306 PetscFunctionBegin; 3307 if (!a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetColumnVec() first"); 3308 if (!a->cvec) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column vector"); 3309 a->vecinuse = 0; 3310 ierr = MatDenseRestoreArrayWrite(A,(PetscScalar**)&a->ptrinuse);CHKERRQ(ierr); 3311 ierr = VecResetArray(a->cvec);CHKERRQ(ierr); 3312 *v = NULL; 3313 PetscFunctionReturn(0); 3314 } 3315 3316 PetscErrorCode MatDenseGetSubMatrix_SeqDense(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 3317 { 3318 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3319 PetscErrorCode ierr; 3320 3321 PetscFunctionBegin; 3322 if (a->vecinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreColumnVec() first"); 3323 if (a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseRestoreSubMatrix() first"); 3324 if (a->cmat && cend-cbegin != a->cmat->cmap->N) { 3325 ierr = MatDestroy(&a->cmat);CHKERRQ(ierr); 3326 } 3327 if (!a->cmat) { 3328 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); 3329 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)a->cmat);CHKERRQ(ierr); 3330 } else { 3331 ierr = MatDensePlaceArray(a->cmat,a->v + (size_t)cbegin * (size_t)a->lda);CHKERRQ(ierr); 3332 } 3333 ierr = MatDenseSetLDA(a->cmat,a->lda);CHKERRQ(ierr); 3334 a->matinuse = cbegin + 1; 3335 *v = a->cmat; 3336 PetscFunctionReturn(0); 3337 } 3338 3339 PetscErrorCode MatDenseRestoreSubMatrix_SeqDense(Mat A,Mat *v) 3340 { 3341 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 3342 PetscErrorCode ierr; 3343 3344 PetscFunctionBegin; 3345 if (!a->matinuse) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Need to call MatDenseGetSubMatrix() first"); 3346 if (!a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Missing internal column matrix"); 3347 if (*v != a->cmat) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not the matrix obtained from MatDenseGetSubMatrix()"); 3348 a->matinuse = 0; 3349 ierr = MatDenseResetArray(a->cmat);CHKERRQ(ierr); 3350 *v = NULL; 3351 PetscFunctionReturn(0); 3352 } 3353 3354 /*MC 3355 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 3356 3357 Options Database Keys: 3358 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 3359 3360 Level: beginner 3361 3362 .seealso: MatCreateSeqDense() 3363 3364 M*/ 3365 PetscErrorCode MatCreate_SeqDense(Mat B) 3366 { 3367 Mat_SeqDense *b; 3368 PetscErrorCode ierr; 3369 PetscMPIInt size; 3370 3371 PetscFunctionBegin; 3372 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRMPI(ierr); 3373 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 3374 3375 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3376 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3377 B->data = (void*)b; 3378 3379 b->roworiented = PETSC_TRUE; 3380 3381 ierr = PetscObjectComposeFunction((PetscObject)B,"MatQRFactor_C",MatQRFactor_SeqDense);CHKERRQ(ierr); 3382 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 3383 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseSetLDA_C",MatDenseSetLDA_SeqDense);CHKERRQ(ierr); 3384 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3385 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 3387 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 3388 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseReplaceArray_C",MatDenseReplaceArray_SeqDense);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3390 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3391 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayWrite_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 3392 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayWrite_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 3393 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 3394 #if defined(PETSC_HAVE_ELEMENTAL) 3395 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 3396 #endif 3397 #if defined(PETSC_HAVE_SCALAPACK) 3398 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_scalapack_C",MatConvert_Dense_ScaLAPACK);CHKERRQ(ierr); 3399 #endif 3400 #if defined(PETSC_HAVE_CUDA) 3401 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqdensecuda_C",MatConvert_SeqDense_SeqDenseCUDA);CHKERRQ(ierr); 3402 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3403 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdensecuda_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3404 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdensecuda_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3405 #endif 3406 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 3407 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqdense_C",MatProductSetFromOptions_SeqAIJ_SeqDense);CHKERRQ(ierr); 3408 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqdense_C",MatProductSetFromOptions_SeqDense);CHKERRQ(ierr); 3409 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3410 ierr = PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqsbaij_seqdense_C",MatProductSetFromOptions_SeqXBAIJ_SeqDense);CHKERRQ(ierr); 3411 3412 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 3413 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 3414 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVec_C",MatDenseGetColumnVec_SeqDense);CHKERRQ(ierr); 3415 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVec_C",MatDenseRestoreColumnVec_SeqDense);CHKERRQ(ierr); 3416 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecRead_C",MatDenseGetColumnVecRead_SeqDense);CHKERRQ(ierr); 3417 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecRead_C",MatDenseRestoreColumnVecRead_SeqDense);CHKERRQ(ierr); 3418 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumnVecWrite_C",MatDenseGetColumnVecWrite_SeqDense);CHKERRQ(ierr); 3419 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumnVecWrite_C",MatDenseRestoreColumnVecWrite_SeqDense);CHKERRQ(ierr); 3420 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetSubMatrix_C",MatDenseGetSubMatrix_SeqDense);CHKERRQ(ierr); 3421 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreSubMatrix_C",MatDenseRestoreSubMatrix_SeqDense);CHKERRQ(ierr); 3422 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 3423 PetscFunctionReturn(0); 3424 } 3425 3426 /*@C 3427 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. 3428 3429 Not Collective 3430 3431 Input Parameters: 3432 + mat - a MATSEQDENSE or MATMPIDENSE matrix 3433 - col - column index 3434 3435 Output Parameter: 3436 . vals - pointer to the data 3437 3438 Level: intermediate 3439 3440 .seealso: MatDenseRestoreColumn() 3441 @*/ 3442 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 3443 { 3444 PetscErrorCode ierr; 3445 3446 PetscFunctionBegin; 3447 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3448 PetscValidLogicalCollectiveInt(A,col,2); 3449 PetscValidPointer(vals,3); 3450 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 3451 PetscFunctionReturn(0); 3452 } 3453 3454 /*@C 3455 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 3456 3457 Not Collective 3458 3459 Input Parameter: 3460 . mat - a MATSEQDENSE or MATMPIDENSE matrix 3461 3462 Output Parameter: 3463 . vals - pointer to the data 3464 3465 Level: intermediate 3466 3467 .seealso: MatDenseGetColumn() 3468 @*/ 3469 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 3470 { 3471 PetscErrorCode ierr; 3472 3473 PetscFunctionBegin; 3474 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3475 PetscValidPointer(vals,2); 3476 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 3477 PetscFunctionReturn(0); 3478 } 3479 3480 /*@C 3481 MatDenseGetColumnVec - Gives read-write access to a column of a dense matrix, represented as a Vec. 3482 3483 Collective 3484 3485 Input Parameters: 3486 + mat - the Mat object 3487 - col - the column index 3488 3489 Output Parameter: 3490 . v - the vector 3491 3492 Notes: 3493 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVec() when the vector is no longer needed. 3494 Use MatDenseGetColumnVecRead() to obtain read-only access or MatDenseGetColumnVecWrite() for write-only access. 3495 3496 Level: intermediate 3497 3498 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3499 @*/ 3500 PetscErrorCode MatDenseGetColumnVec(Mat A,PetscInt col,Vec *v) 3501 { 3502 PetscErrorCode ierr; 3503 3504 PetscFunctionBegin; 3505 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3506 PetscValidType(A,1); 3507 PetscValidLogicalCollectiveInt(A,col,2); 3508 PetscValidPointer(v,3); 3509 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3510 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3511 ierr = PetscUseMethod(A,"MatDenseGetColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3512 PetscFunctionReturn(0); 3513 } 3514 3515 /*@C 3516 MatDenseRestoreColumnVec - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVec(). 3517 3518 Collective 3519 3520 Input Parameters: 3521 + mat - the Mat object 3522 . col - the column index 3523 - v - the Vec object 3524 3525 Level: intermediate 3526 3527 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3528 @*/ 3529 PetscErrorCode MatDenseRestoreColumnVec(Mat A,PetscInt col,Vec *v) 3530 { 3531 PetscErrorCode ierr; 3532 3533 PetscFunctionBegin; 3534 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3535 PetscValidType(A,1); 3536 PetscValidLogicalCollectiveInt(A,col,2); 3537 PetscValidPointer(v,3); 3538 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3539 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3540 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVec_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3541 PetscFunctionReturn(0); 3542 } 3543 3544 /*@C 3545 MatDenseGetColumnVecRead - Gives read-only access to a column of a dense matrix, represented as a Vec. 3546 3547 Collective 3548 3549 Input Parameters: 3550 + mat - the Mat object 3551 - col - the column index 3552 3553 Output Parameter: 3554 . v - the vector 3555 3556 Notes: 3557 The vector is owned by PETSc and users cannot modify it. 3558 Users need to call MatDenseRestoreColumnVecRead() when the vector is no longer needed. 3559 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecWrite() for write-only access. 3560 3561 Level: intermediate 3562 3563 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3564 @*/ 3565 PetscErrorCode MatDenseGetColumnVecRead(Mat A,PetscInt col,Vec *v) 3566 { 3567 PetscErrorCode ierr; 3568 3569 PetscFunctionBegin; 3570 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3571 PetscValidType(A,1); 3572 PetscValidLogicalCollectiveInt(A,col,2); 3573 PetscValidPointer(v,3); 3574 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3575 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3576 ierr = PetscUseMethod(A,"MatDenseGetColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3577 PetscFunctionReturn(0); 3578 } 3579 3580 /*@C 3581 MatDenseRestoreColumnVecRead - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecRead(). 3582 3583 Collective 3584 3585 Input Parameters: 3586 + mat - the Mat object 3587 . col - the column index 3588 - v - the Vec object 3589 3590 Level: intermediate 3591 3592 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecWrite() 3593 @*/ 3594 PetscErrorCode MatDenseRestoreColumnVecRead(Mat A,PetscInt col,Vec *v) 3595 { 3596 PetscErrorCode ierr; 3597 3598 PetscFunctionBegin; 3599 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3600 PetscValidType(A,1); 3601 PetscValidLogicalCollectiveInt(A,col,2); 3602 PetscValidPointer(v,3); 3603 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3604 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3605 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecRead_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3606 PetscFunctionReturn(0); 3607 } 3608 3609 /*@C 3610 MatDenseGetColumnVecWrite - Gives write-only access to a column of a dense matrix, represented as a Vec. 3611 3612 Collective 3613 3614 Input Parameters: 3615 + mat - the Mat object 3616 - col - the column index 3617 3618 Output Parameter: 3619 . v - the vector 3620 3621 Notes: 3622 The vector is owned by PETSc. Users need to call MatDenseRestoreColumnVecWrite() when the vector is no longer needed. 3623 Use MatDenseGetColumnVec() to obtain read-write access or MatDenseGetColumnVecRead() for read-only access. 3624 3625 Level: intermediate 3626 3627 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead(), MatDenseRestoreColumnVecWrite() 3628 @*/ 3629 PetscErrorCode MatDenseGetColumnVecWrite(Mat A,PetscInt col,Vec *v) 3630 { 3631 PetscErrorCode ierr; 3632 3633 PetscFunctionBegin; 3634 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3635 PetscValidType(A,1); 3636 PetscValidLogicalCollectiveInt(A,col,2); 3637 PetscValidPointer(v,3); 3638 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3639 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3640 ierr = PetscUseMethod(A,"MatDenseGetColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3641 PetscFunctionReturn(0); 3642 } 3643 3644 /*@C 3645 MatDenseRestoreColumnVecWrite - Returns access to a column of a dense matrix obtained from MatDenseGetColumnVecWrite(). 3646 3647 Collective 3648 3649 Input Parameters: 3650 + mat - the Mat object 3651 . col - the column index 3652 - v - the Vec object 3653 3654 Level: intermediate 3655 3656 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseGetColumnVecRead(), MatDenseGetColumnVecWrite(), MatDenseRestoreColumnVec(), MatDenseRestoreColumnVecRead() 3657 @*/ 3658 PetscErrorCode MatDenseRestoreColumnVecWrite(Mat A,PetscInt col,Vec *v) 3659 { 3660 PetscErrorCode ierr; 3661 3662 PetscFunctionBegin; 3663 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3664 PetscValidType(A,1); 3665 PetscValidLogicalCollectiveInt(A,col,2); 3666 PetscValidPointer(v,3); 3667 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3668 if (col < 0 || col > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid col %D, should be in [0,%D)",col,A->cmap->N); 3669 ierr = PetscUseMethod(A,"MatDenseRestoreColumnVecWrite_C",(Mat,PetscInt,Vec*),(A,col,v));CHKERRQ(ierr); 3670 PetscFunctionReturn(0); 3671 } 3672 3673 /*@C 3674 MatDenseGetSubMatrix - Gives access to a block of columns of a dense matrix, represented as a Mat. 3675 3676 Collective 3677 3678 Input Parameters: 3679 + mat - the Mat object 3680 . cbegin - the first index in the block 3681 - cend - the last index in the block 3682 3683 Output Parameter: 3684 . v - the matrix 3685 3686 Notes: 3687 The matrix is owned by PETSc. Users need to call MatDenseRestoreSubMatrix() when the matrix is no longer needed. 3688 3689 Level: intermediate 3690 3691 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseRestoreSubMatrix() 3692 @*/ 3693 PetscErrorCode MatDenseGetSubMatrix(Mat A,PetscInt cbegin,PetscInt cend,Mat *v) 3694 { 3695 PetscErrorCode ierr; 3696 3697 PetscFunctionBegin; 3698 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3699 PetscValidType(A,1); 3700 PetscValidLogicalCollectiveInt(A,cbegin,2); 3701 PetscValidLogicalCollectiveInt(A,cend,3); 3702 PetscValidPointer(v,4); 3703 if (!A->preallocated) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Matrix not preallocated"); 3704 if (cbegin < 0 || cbegin > A->cmap->N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cbegin %D, should be in [0,%D)",cbegin,A->cmap->N); 3705 if (cend < cbegin || cend > A->cmap->N) SETERRQ3(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Invalid cend %D, should be in [%D,%D)",cend,cbegin,A->cmap->N); 3706 ierr = PetscUseMethod(A,"MatDenseGetSubMatrix_C",(Mat,PetscInt,PetscInt,Mat*),(A,cbegin,cend,v));CHKERRQ(ierr); 3707 PetscFunctionReturn(0); 3708 } 3709 3710 /*@C 3711 MatDenseRestoreSubMatrix - Returns access to a block of columns of a dense matrix obtained from MatDenseGetSubMatrix(). 3712 3713 Collective 3714 3715 Input Parameters: 3716 + mat - the Mat object 3717 - v - the Mat object 3718 3719 Level: intermediate 3720 3721 .seealso: MATDENSE, MATDENSECUDA, MatDenseGetColumnVec(), MatDenseRestoreColumnVec(), MatDenseGetSubMatrix() 3722 @*/ 3723 PetscErrorCode MatDenseRestoreSubMatrix(Mat A,Mat *v) 3724 { 3725 PetscErrorCode ierr; 3726 3727 PetscFunctionBegin; 3728 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3729 PetscValidType(A,1); 3730 PetscValidPointer(v,2); 3731 ierr = PetscUseMethod(A,"MatDenseRestoreSubMatrix_C",(Mat,Mat*),(A,v));CHKERRQ(ierr); 3732 PetscFunctionReturn(0); 3733 } 3734