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