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