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 static 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 16 PetscFunctionBegin; 17 if (A->rmap->n != A->cmap->n) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot symmetrize a rectangular matrix"); 18 if (!hermitian) { 19 for (k=0;k<n;k++) { 20 for (j=k;j<n;j++) { 21 mat->v[j*mat->lda + k] = mat->v[k*mat->lda + j]; 22 } 23 } 24 } else { 25 for (k=0;k<n;k++) { 26 for (j=k;j<n;j++) { 27 mat->v[j*mat->lda + k] = PetscConj(mat->v[k*mat->lda + j]); 28 } 29 } 30 } 31 PetscFunctionReturn(0); 32 } 33 34 PETSC_EXTERN PetscErrorCode MatSeqDenseInvertFactors_Private(Mat A) 35 { 36 #if defined(PETSC_MISSING_LAPACK_POTRF) 37 PetscFunctionBegin; 38 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 39 #else 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); /* TODO CHECK FLOPS */ 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); /* TODO CHECK FLOPS */ 83 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 84 #endif 85 86 A->ops->solve = NULL; 87 A->ops->matsolve = NULL; 88 A->ops->solvetranspose = NULL; 89 A->ops->matsolvetranspose = NULL; 90 A->ops->solveadd = NULL; 91 A->ops->solvetransposeadd = NULL; 92 A->factortype = MAT_FACTOR_NONE; 93 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 PetscErrorCode MatZeroRowsColumns_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 98 { 99 PetscErrorCode ierr; 100 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 101 PetscInt m = l->lda, n = A->cmap->n,r = A->rmap->n, i,j; 102 PetscScalar *slot,*bb; 103 const PetscScalar *xx; 104 105 PetscFunctionBegin; 106 #if defined(PETSC_USE_DEBUG) 107 for (i=0; i<N; i++) { 108 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 109 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); 110 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); 111 } 112 #endif 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 for (i=0; i<N; i++) { 132 slot = l->v + rows[i]*m; 133 ierr = PetscArrayzero(slot,r);CHKERRQ(ierr); 134 } 135 for (i=0; i<N; i++) { 136 slot = l->v + rows[i]; 137 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 138 } 139 if (diag != 0.0) { 140 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 141 for (i=0; i<N; i++) { 142 slot = l->v + (m+1)*rows[i]; 143 *slot = diag; 144 } 145 } 146 PetscFunctionReturn(0); 147 } 148 149 PetscErrorCode MatPtAPNumeric_SeqDense_SeqDense(Mat A,Mat P,Mat C) 150 { 151 Mat_SeqDense *c = (Mat_SeqDense*)(C->data); 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,P,c->ptapwork);CHKERRQ(ierr); 156 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(P,c->ptapwork,C);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode MatPtAPSymbolic_SeqDense_SeqDense(Mat A,Mat P,PetscReal fill,Mat *C) 161 { 162 Mat_SeqDense *c; 163 PetscErrorCode ierr; 164 165 PetscFunctionBegin; 166 ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)A),P->cmap->N,P->cmap->N,NULL,C);CHKERRQ(ierr); 167 c = (Mat_SeqDense*)((*C)->data); 168 ierr = MatCreateSeqDense(PetscObjectComm((PetscObject)A),A->rmap->N,P->cmap->N,NULL,&c->ptapwork);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 172 PETSC_INTERN PetscErrorCode MatPtAP_SeqDense_SeqDense(Mat A,Mat P,MatReuse reuse,PetscReal fill,Mat *C) 173 { 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 if (reuse == MAT_INITIAL_MATRIX) { 178 ierr = MatPtAPSymbolic_SeqDense_SeqDense(A,P,fill,C);CHKERRQ(ierr); 179 } 180 ierr = PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 181 ierr = (*(*C)->ops->ptapnumeric)(A,P,*C);CHKERRQ(ierr); 182 ierr = PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 187 { 188 Mat B = NULL; 189 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 190 Mat_SeqDense *b; 191 PetscErrorCode ierr; 192 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->N,i; 193 MatScalar *av=a->a; 194 PetscBool isseqdense; 195 196 PetscFunctionBegin; 197 if (reuse == MAT_REUSE_MATRIX) { 198 ierr = PetscObjectTypeCompare((PetscObject)*newmat,MATSEQDENSE,&isseqdense);CHKERRQ(ierr); 199 if (!isseqdense) SETERRQ1(PetscObjectComm((PetscObject)*newmat),PETSC_ERR_USER,"Cannot reuse matrix of type %s",((PetscObject)(*newmat))->type); 200 } 201 if (reuse != MAT_REUSE_MATRIX) { 202 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 203 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 204 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 205 ierr = MatSeqDenseSetPreallocation(B,NULL);CHKERRQ(ierr); 206 b = (Mat_SeqDense*)(B->data); 207 } else { 208 b = (Mat_SeqDense*)((*newmat)->data); 209 ierr = PetscArrayzero(b->v,m*n);CHKERRQ(ierr); 210 } 211 for (i=0; i<m; i++) { 212 PetscInt j; 213 for (j=0;j<ai[1]-ai[0];j++) { 214 b->v[*aj*m+i] = *av; 215 aj++; 216 av++; 217 } 218 ai++; 219 } 220 221 if (reuse == MAT_INPLACE_MATRIX) { 222 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 223 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 224 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 225 } else { 226 if (B) *newmat = B; 227 ierr = MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 228 ierr = MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 229 } 230 PetscFunctionReturn(0); 231 } 232 233 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 234 { 235 Mat B; 236 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 237 PetscErrorCode ierr; 238 PetscInt i, j; 239 PetscInt *rows, *nnz; 240 MatScalar *aa = a->v, *vals; 241 242 PetscFunctionBegin; 243 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 244 ierr = MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 245 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 246 ierr = PetscCalloc3(A->rmap->n,&rows,A->rmap->n,&nnz,A->rmap->n,&vals);CHKERRQ(ierr); 247 for (j=0; j<A->cmap->n; j++) { 248 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) ++nnz[i]; 249 aa += a->lda; 250 } 251 ierr = MatSeqAIJSetPreallocation(B,PETSC_DETERMINE,nnz);CHKERRQ(ierr); 252 aa = a->v; 253 for (j=0; j<A->cmap->n; j++) { 254 PetscInt numRows = 0; 255 for (i=0; i<A->rmap->n; i++) if (aa[i] != 0.0 || i == j) {rows[numRows] = i; vals[numRows++] = aa[i];} 256 ierr = MatSetValues(B,numRows,rows,1,&j,vals,INSERT_VALUES);CHKERRQ(ierr); 257 aa += a->lda; 258 } 259 ierr = PetscFree3(rows,nnz,vals);CHKERRQ(ierr); 260 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 261 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 262 263 if (reuse == MAT_INPLACE_MATRIX) { 264 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 265 } else { 266 *newmat = B; 267 } 268 PetscFunctionReturn(0); 269 } 270 271 static PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 272 { 273 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 274 PetscScalar oalpha = alpha; 275 PetscInt j; 276 PetscBLASInt N,m,ldax,lday,one = 1; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = PetscBLASIntCast(X->rmap->n*X->cmap->n,&N);CHKERRQ(ierr); 281 ierr = PetscBLASIntCast(X->rmap->n,&m);CHKERRQ(ierr); 282 ierr = PetscBLASIntCast(x->lda,&ldax);CHKERRQ(ierr); 283 ierr = PetscBLASIntCast(y->lda,&lday);CHKERRQ(ierr); 284 if (ldax>m || lday>m) { 285 for (j=0; j<X->cmap->n; j++) { 286 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one)); 287 } 288 } else { 289 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one)); 290 } 291 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 292 ierr = PetscLogFlops(PetscMax(2*N-1,0));CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 296 static PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 297 { 298 PetscInt N = A->rmap->n*A->cmap->n; 299 300 PetscFunctionBegin; 301 info->block_size = 1.0; 302 info->nz_allocated = (double)N; 303 info->nz_used = (double)N; 304 info->nz_unneeded = (double)0; 305 info->assemblies = (double)A->num_ass; 306 info->mallocs = 0; 307 info->memory = ((PetscObject)A)->mem; 308 info->fill_ratio_given = 0; 309 info->fill_ratio_needed = 0; 310 info->factor_mallocs = 0; 311 PetscFunctionReturn(0); 312 } 313 314 static PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 315 { 316 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 317 PetscScalar oalpha = alpha; 318 PetscErrorCode ierr; 319 PetscBLASInt one = 1,j,nz,lda; 320 321 PetscFunctionBegin; 322 ierr = PetscBLASIntCast(a->lda,&lda);CHKERRQ(ierr); 323 if (lda>A->rmap->n) { 324 ierr = PetscBLASIntCast(A->rmap->n,&nz);CHKERRQ(ierr); 325 for (j=0; j<A->cmap->n; j++) { 326 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v+j*lda,&one)); 327 } 328 } else { 329 ierr = PetscBLASIntCast(A->rmap->n*A->cmap->n,&nz);CHKERRQ(ierr); 330 PetscStackCallBLAS("BLASscal",BLASscal_(&nz,&oalpha,a->v,&one)); 331 } 332 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 333 PetscFunctionReturn(0); 334 } 335 336 static PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscBool *fl) 337 { 338 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 339 PetscInt i,j,m = A->rmap->n,N; 340 PetscScalar *v = a->v; 341 342 PetscFunctionBegin; 343 *fl = PETSC_FALSE; 344 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 345 N = a->lda; 346 347 for (i=0; i<m; i++) { 348 for (j=i+1; j<m; j++) { 349 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 350 } 351 } 352 *fl = PETSC_TRUE; 353 PetscFunctionReturn(0); 354 } 355 356 static PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat newi,Mat A,MatDuplicateOption cpvalues) 357 { 358 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 359 PetscErrorCode ierr; 360 PetscInt lda = (PetscInt)mat->lda,j,m; 361 362 PetscFunctionBegin; 363 ierr = PetscLayoutReference(A->rmap,&newi->rmap);CHKERRQ(ierr); 364 ierr = PetscLayoutReference(A->cmap,&newi->cmap);CHKERRQ(ierr); 365 ierr = MatSeqDenseSetPreallocation(newi,NULL);CHKERRQ(ierr); 366 if (cpvalues == MAT_COPY_VALUES) { 367 l = (Mat_SeqDense*)newi->data; 368 if (lda>A->rmap->n) { 369 m = A->rmap->n; 370 for (j=0; j<A->cmap->n; j++) { 371 ierr = PetscArraycpy(l->v+j*m,mat->v+j*lda,m);CHKERRQ(ierr); 372 } 373 } else { 374 ierr = PetscArraycpy(l->v,mat->v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 375 } 376 } 377 newi->assembled = PETSC_TRUE; 378 PetscFunctionReturn(0); 379 } 380 381 static PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 382 { 383 PetscErrorCode ierr; 384 385 PetscFunctionBegin; 386 ierr = MatCreate(PetscObjectComm((PetscObject)A),newmat);CHKERRQ(ierr); 387 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 388 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 389 ierr = MatDuplicateNoCreate_SeqDense(*newmat,A,cpvalues);CHKERRQ(ierr); 390 PetscFunctionReturn(0); 391 } 392 393 394 static PetscErrorCode MatLUFactor_SeqDense(Mat,IS,IS,const MatFactorInfo*); 395 396 static PetscErrorCode MatLUFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 397 { 398 MatFactorInfo info; 399 PetscErrorCode ierr; 400 401 PetscFunctionBegin; 402 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 403 ierr = MatLUFactor_SeqDense(fact,0,0,&info);CHKERRQ(ierr); 404 PetscFunctionReturn(0); 405 } 406 407 static PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 408 { 409 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 410 PetscErrorCode ierr; 411 const PetscScalar *x; 412 PetscScalar *y; 413 PetscBLASInt one = 1,info,m; 414 415 PetscFunctionBegin; 416 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 417 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 418 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 419 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 420 if (A->factortype == MAT_FACTOR_LU) { 421 #if defined(PETSC_MISSING_LAPACK_GETRS) 422 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 423 #else 424 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 425 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 426 ierr = PetscFPTrapPop();CHKERRQ(ierr); 427 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 428 #endif 429 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 430 #if defined(PETSC_MISSING_LAPACK_POTRS) 431 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 432 #else 433 if (A->spd) { 434 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 435 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 436 ierr = PetscFPTrapPop();CHKERRQ(ierr); 437 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 438 #if defined(PETSC_USE_COMPLEX) 439 } else if (A->hermitian) { 440 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 441 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 442 ierr = PetscFPTrapPop();CHKERRQ(ierr); 443 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 444 #endif 445 } else { /* symmetric case */ 446 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 447 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 448 ierr = PetscFPTrapPop();CHKERRQ(ierr); 449 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 450 } 451 #endif 452 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 453 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 454 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 455 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 static PetscErrorCode MatMatSolve_SeqDense(Mat A,Mat B,Mat X) 460 { 461 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 462 PetscErrorCode ierr; 463 const PetscScalar *b; 464 PetscScalar *x; 465 PetscInt n; 466 PetscBLASInt nrhs,info,m; 467 PetscBool flg; 468 469 PetscFunctionBegin; 470 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 471 ierr = PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 472 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix B must be MATDENSE matrix"); 473 ierr = PetscObjectTypeCompareAny((PetscObject)X,&flg,MATSEQDENSE,MATMPIDENSE,NULL);CHKERRQ(ierr); 474 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix X must be MATDENSE matrix"); 475 476 ierr = MatGetSize(B,NULL,&n);CHKERRQ(ierr); 477 ierr = PetscBLASIntCast(n,&nrhs);CHKERRQ(ierr); 478 ierr = MatDenseGetArrayRead(B,&b);CHKERRQ(ierr); 479 ierr = MatDenseGetArray(X,&x);CHKERRQ(ierr); 480 481 ierr = PetscArraycpy(x,b,m*nrhs);CHKERRQ(ierr); 482 483 if (A->factortype == MAT_FACTOR_LU) { 484 #if defined(PETSC_MISSING_LAPACK_GETRS) 485 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 486 #else 487 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 488 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 489 ierr = PetscFPTrapPop();CHKERRQ(ierr); 490 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"GETRS - Bad solve"); 491 #endif 492 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 493 #if defined(PETSC_MISSING_LAPACK_POTRS) 494 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 495 #else 496 if (A->spd) { 497 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 498 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&nrhs,mat->v,&mat->lda,x,&m,&info)); 499 ierr = PetscFPTrapPop();CHKERRQ(ierr); 500 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 501 #if defined(PETSC_USE_COMPLEX) 502 } else if (A->hermitian) { 503 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 504 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 505 ierr = PetscFPTrapPop();CHKERRQ(ierr); 506 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"HETRS Bad solve"); 507 #endif 508 } else { /* symmetric case */ 509 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 510 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&nrhs,mat->v,&mat->lda,mat->pivots,x,&m,&info)); 511 ierr = PetscFPTrapPop();CHKERRQ(ierr); 512 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 513 } 514 #endif 515 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 516 517 ierr = MatDenseRestoreArrayRead(B,&b);CHKERRQ(ierr); 518 ierr = MatDenseRestoreArray(X,&x);CHKERRQ(ierr); 519 ierr = PetscLogFlops(nrhs*(2.0*m*m - m));CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 523 static PetscErrorCode MatConjugate_SeqDense(Mat); 524 525 static PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 526 { 527 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 528 PetscErrorCode ierr; 529 const PetscScalar *x; 530 PetscScalar *y; 531 PetscBLASInt one = 1,info,m; 532 533 PetscFunctionBegin; 534 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 535 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 536 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 537 ierr = PetscArraycpy(y,x,A->rmap->n);CHKERRQ(ierr); 538 if (A->factortype == MAT_FACTOR_LU) { 539 #if defined(PETSC_MISSING_LAPACK_GETRS) 540 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 541 #else 542 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 543 PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 544 ierr = PetscFPTrapPop();CHKERRQ(ierr); 545 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS - Bad solve"); 546 #endif 547 } else if (A->factortype == MAT_FACTOR_CHOLESKY) { 548 #if defined(PETSC_MISSING_LAPACK_POTRS) 549 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 550 #else 551 if (A->spd) { 552 #if defined(PETSC_USE_COMPLEX) 553 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 554 #endif 555 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 556 PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info)); 557 ierr = PetscFPTrapPop();CHKERRQ(ierr); 558 #if defined(PETSC_USE_COMPLEX) 559 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 560 #endif 561 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"POTRS Bad solve"); 562 #if defined(PETSC_USE_COMPLEX) 563 } else if (A->hermitian) { 564 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 565 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 566 PetscStackCallBLAS("LAPACKhetrs",LAPACKhetrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 567 ierr = PetscFPTrapPop();CHKERRQ(ierr); 568 ierr = MatConjugate_SeqDense(A);CHKERRQ(ierr); 569 #endif 570 } else { /* symmetric case */ 571 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 572 PetscStackCallBLAS("LAPACKsytrs",LAPACKsytrs_("L",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info)); 573 ierr = PetscFPTrapPop();CHKERRQ(ierr); 574 if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"SYTRS Bad solve"); 575 } 576 #endif 577 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 578 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 579 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 580 ierr = PetscLogFlops(2.0*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 static PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 585 { 586 PetscErrorCode ierr; 587 const PetscScalar *x; 588 PetscScalar *y,sone = 1.0; 589 Vec tmp = 0; 590 591 PetscFunctionBegin; 592 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 593 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 594 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 595 if (yy == zz) { 596 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 597 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 598 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 599 } 600 ierr = MatSolve_SeqDense(A,xx,yy);CHKERRQ(ierr); 601 if (tmp) { 602 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 603 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 604 } else { 605 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 606 } 607 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 608 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 609 ierr = PetscLogFlops(A->cmap->n);CHKERRQ(ierr); 610 PetscFunctionReturn(0); 611 } 612 613 static PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 614 { 615 PetscErrorCode ierr; 616 const PetscScalar *x; 617 PetscScalar *y,sone = 1.0; 618 Vec tmp = NULL; 619 620 PetscFunctionBegin; 621 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 622 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 623 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 624 if (yy == zz) { 625 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 626 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)tmp);CHKERRQ(ierr); 627 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 628 } 629 ierr = MatSolveTranspose_SeqDense(A,xx,yy);CHKERRQ(ierr); 630 if (tmp) { 631 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 632 ierr = VecDestroy(&tmp);CHKERRQ(ierr); 633 } else { 634 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 635 } 636 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 637 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 638 ierr = PetscLogFlops(A->cmap->n);CHKERRQ(ierr); 639 PetscFunctionReturn(0); 640 } 641 642 /* ---------------------------------------------------------------*/ 643 /* COMMENT: I have chosen to hide row permutation in the pivots, 644 rather than put it in the Mat->row slot.*/ 645 static PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,const MatFactorInfo *minfo) 646 { 647 #if defined(PETSC_MISSING_LAPACK_GETRF) 648 PetscFunctionBegin; 649 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 650 #else 651 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 652 PetscErrorCode ierr; 653 PetscBLASInt n,m,info; 654 655 PetscFunctionBegin; 656 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 657 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 658 if (!mat->pivots) { 659 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 660 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 661 } 662 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 663 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 664 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info)); 665 ierr = PetscFPTrapPop();CHKERRQ(ierr); 666 667 if (info<0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to LU factorization"); 668 if (info>0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 669 670 A->ops->solve = MatSolve_SeqDense; 671 A->ops->matsolve = MatMatSolve_SeqDense; 672 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 673 A->ops->solveadd = MatSolveAdd_SeqDense; 674 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 675 A->factortype = MAT_FACTOR_LU; 676 677 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 678 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 679 680 ierr = PetscLogFlops((2.0*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 681 #endif 682 PetscFunctionReturn(0); 683 } 684 685 /* Cholesky as L*L^T or L*D*L^T and the symmetric/hermitian complex variants */ 686 static PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,const MatFactorInfo *factinfo) 687 { 688 #if defined(PETSC_MISSING_LAPACK_POTRF) 689 PetscFunctionBegin; 690 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 691 #else 692 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 693 PetscErrorCode ierr; 694 PetscBLASInt info,n; 695 696 PetscFunctionBegin; 697 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 698 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 699 if (A->spd) { 700 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 701 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info)); 702 ierr = PetscFPTrapPop();CHKERRQ(ierr); 703 #if defined(PETSC_USE_COMPLEX) 704 } else if (A->hermitian) { 705 if (!mat->pivots) { 706 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 707 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 708 } 709 if (!mat->fwork) { 710 PetscScalar dummy; 711 712 mat->lfwork = -1; 713 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 714 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 715 ierr = PetscFPTrapPop();CHKERRQ(ierr); 716 mat->lfwork = (PetscInt)PetscRealPart(dummy); 717 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 718 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 719 } 720 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 721 PetscStackCallBLAS("LAPACKhetrf",LAPACKhetrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 722 ierr = PetscFPTrapPop();CHKERRQ(ierr); 723 #endif 724 } else { /* symmetric case */ 725 if (!mat->pivots) { 726 ierr = PetscMalloc1(A->rmap->n,&mat->pivots);CHKERRQ(ierr); 727 ierr = PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 728 } 729 if (!mat->fwork) { 730 PetscScalar dummy; 731 732 mat->lfwork = -1; 733 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 734 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,&dummy,&mat->lfwork,&info)); 735 ierr = PetscFPTrapPop();CHKERRQ(ierr); 736 mat->lfwork = (PetscInt)PetscRealPart(dummy); 737 ierr = PetscMalloc1(mat->lfwork,&mat->fwork);CHKERRQ(ierr); 738 ierr = PetscLogObjectMemory((PetscObject)A,mat->lfwork*sizeof(PetscBLASInt));CHKERRQ(ierr); 739 } 740 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 741 PetscStackCallBLAS("LAPACKsytrf",LAPACKsytrf_("L",&n,mat->v,&mat->lda,mat->pivots,mat->fwork,&mat->lfwork,&info)); 742 ierr = PetscFPTrapPop();CHKERRQ(ierr); 743 } 744 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 745 746 A->ops->solve = MatSolve_SeqDense; 747 A->ops->matsolve = MatMatSolve_SeqDense; 748 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 749 A->ops->solveadd = MatSolveAdd_SeqDense; 750 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 751 A->factortype = MAT_FACTOR_CHOLESKY; 752 753 ierr = PetscFree(A->solvertype);CHKERRQ(ierr); 754 ierr = PetscStrallocpy(MATSOLVERPETSC,&A->solvertype);CHKERRQ(ierr); 755 756 ierr = PetscLogFlops((1.0*A->cmap->n*A->cmap->n*A->cmap->n)/3.0);CHKERRQ(ierr); 757 #endif 758 PetscFunctionReturn(0); 759 } 760 761 762 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat fact,Mat A,const MatFactorInfo *info_dummy) 763 { 764 PetscErrorCode ierr; 765 MatFactorInfo info; 766 767 PetscFunctionBegin; 768 info.fill = 1.0; 769 770 ierr = MatDuplicateNoCreate_SeqDense(fact,A,MAT_COPY_VALUES);CHKERRQ(ierr); 771 ierr = MatCholeskyFactor_SeqDense(fact,0,&info);CHKERRQ(ierr); 772 PetscFunctionReturn(0); 773 } 774 775 static PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,const MatFactorInfo *info) 776 { 777 PetscFunctionBegin; 778 fact->assembled = PETSC_TRUE; 779 fact->preallocated = PETSC_TRUE; 780 fact->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 781 fact->ops->solve = MatSolve_SeqDense; 782 fact->ops->matsolve = MatMatSolve_SeqDense; 783 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 784 fact->ops->solveadd = MatSolveAdd_SeqDense; 785 fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 786 PetscFunctionReturn(0); 787 } 788 789 static PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat fact,Mat A,IS row,IS col,const MatFactorInfo *info) 790 { 791 PetscFunctionBegin; 792 fact->preallocated = PETSC_TRUE; 793 fact->assembled = PETSC_TRUE; 794 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 795 fact->ops->solve = MatSolve_SeqDense; 796 fact->ops->matsolve = MatMatSolve_SeqDense; 797 fact->ops->solvetranspose = MatSolveTranspose_SeqDense; 798 fact->ops->solveadd = MatSolveAdd_SeqDense; 799 fact->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 800 PetscFunctionReturn(0); 801 } 802 803 PETSC_INTERN PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 804 { 805 PetscErrorCode ierr; 806 807 PetscFunctionBegin; 808 ierr = MatCreate(PetscObjectComm((PetscObject)A),fact);CHKERRQ(ierr); 809 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 810 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 811 if (ftype == MAT_FACTOR_LU) { 812 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 813 } else { 814 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 815 } 816 (*fact)->factortype = ftype; 817 818 ierr = PetscFree((*fact)->solvertype);CHKERRQ(ierr); 819 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*fact)->solvertype);CHKERRQ(ierr); 820 PetscFunctionReturn(0); 821 } 822 823 /* ------------------------------------------------------------------*/ 824 static PetscErrorCode MatSOR_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 825 { 826 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 827 PetscScalar *x,*v = mat->v,zero = 0.0,xt; 828 const PetscScalar *b; 829 PetscErrorCode ierr; 830 PetscInt m = A->rmap->n,i; 831 PetscBLASInt o = 1,bm; 832 833 PetscFunctionBegin; 834 if (shift == -1) shift = 0.0; /* negative shift indicates do not error on zero diagonal; this code never zeros on zero diagonal */ 835 ierr = PetscBLASIntCast(m,&bm);CHKERRQ(ierr); 836 if (flag & SOR_ZERO_INITIAL_GUESS) { 837 /* this is a hack fix, should have another version without the second BLASdotu */ 838 ierr = VecSet(xx,zero);CHKERRQ(ierr); 839 } 840 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 841 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 842 its = its*lits; 843 if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 844 while (its--) { 845 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 846 for (i=0; i<m; i++) { 847 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 848 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 849 } 850 } 851 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 852 for (i=m-1; i>=0; i--) { 853 PetscStackCallBLAS("BLASdotu",xt = b[i] - BLASdotu_(&bm,v+i,&bm,x,&o)); 854 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 855 } 856 } 857 } 858 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 859 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 860 PetscFunctionReturn(0); 861 } 862 863 /* -----------------------------------------------------------------*/ 864 PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 865 { 866 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 867 const PetscScalar *v = mat->v,*x; 868 PetscScalar *y; 869 PetscErrorCode ierr; 870 PetscBLASInt m, n,_One=1; 871 PetscScalar _DOne=1.0,_DZero=0.0; 872 873 PetscFunctionBegin; 874 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 875 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 876 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 877 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 878 if (!A->rmap->n || !A->cmap->n) { 879 PetscBLASInt i; 880 for (i=0; i<n; i++) y[i] = 0.0; 881 } else { 882 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One)); 883 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 884 } 885 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 886 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 PETSC_INTERN PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 891 { 892 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 893 PetscScalar *y,_DOne=1.0,_DZero=0.0; 894 PetscErrorCode ierr; 895 PetscBLASInt m, n, _One=1; 896 const PetscScalar *v = mat->v,*x; 897 898 PetscFunctionBegin; 899 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 900 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 901 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 902 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 903 if (!A->rmap->n || !A->cmap->n) { 904 PetscBLASInt i; 905 for (i=0; i<m; i++) y[i] = 0.0; 906 } else { 907 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One)); 908 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 909 } 910 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 911 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 912 PetscFunctionReturn(0); 913 } 914 915 PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 916 { 917 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 918 const PetscScalar *v = mat->v,*x; 919 PetscScalar *y,_DOne=1.0; 920 PetscErrorCode ierr; 921 PetscBLASInt m, n, _One=1; 922 923 PetscFunctionBegin; 924 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 925 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 926 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 927 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 928 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 929 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 930 PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 931 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 932 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 933 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 934 PetscFunctionReturn(0); 935 } 936 937 static PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 938 { 939 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 940 const PetscScalar *v = mat->v,*x; 941 PetscScalar *y; 942 PetscErrorCode ierr; 943 PetscBLASInt m, n, _One=1; 944 PetscScalar _DOne=1.0; 945 946 PetscFunctionBegin; 947 ierr = PetscBLASIntCast(A->rmap->n,&m);CHKERRQ(ierr); 948 ierr = PetscBLASIntCast(A->cmap->n,&n);CHKERRQ(ierr); 949 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 950 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 951 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 952 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 953 PetscStackCallBLAS("BLASgemv",BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One)); 954 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 955 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 956 ierr = PetscLogFlops(2.0*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 957 PetscFunctionReturn(0); 958 } 959 960 /* -----------------------------------------------------------------*/ 961 static PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 962 { 963 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 964 PetscScalar *v; 965 PetscErrorCode ierr; 966 PetscInt i; 967 968 PetscFunctionBegin; 969 *ncols = A->cmap->n; 970 if (cols) { 971 ierr = PetscMalloc1(A->cmap->n+1,cols);CHKERRQ(ierr); 972 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 973 } 974 if (vals) { 975 ierr = PetscMalloc1(A->cmap->n+1,vals);CHKERRQ(ierr); 976 v = mat->v + row; 977 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 978 } 979 PetscFunctionReturn(0); 980 } 981 982 static PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 983 { 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 988 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 989 PetscFunctionReturn(0); 990 } 991 /* ----------------------------------------------------------------*/ 992 static PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 993 { 994 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 995 PetscInt i,j,idx=0; 996 997 PetscFunctionBegin; 998 if (!mat->roworiented) { 999 if (addv == INSERT_VALUES) { 1000 for (j=0; j<n; j++) { 1001 if (indexn[j] < 0) {idx += m; continue;} 1002 #if defined(PETSC_USE_DEBUG) 1003 if (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); 1004 #endif 1005 for (i=0; i<m; i++) { 1006 if (indexm[i] < 0) {idx++; continue;} 1007 #if defined(PETSC_USE_DEBUG) 1008 if (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); 1009 #endif 1010 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 1011 } 1012 } 1013 } else { 1014 for (j=0; j<n; j++) { 1015 if (indexn[j] < 0) {idx += m; continue;} 1016 #if defined(PETSC_USE_DEBUG) 1017 if (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); 1018 #endif 1019 for (i=0; i<m; i++) { 1020 if (indexm[i] < 0) {idx++; continue;} 1021 #if defined(PETSC_USE_DEBUG) 1022 if (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); 1023 #endif 1024 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1025 } 1026 } 1027 } 1028 } else { 1029 if (addv == INSERT_VALUES) { 1030 for (i=0; i<m; i++) { 1031 if (indexm[i] < 0) { idx += n; continue;} 1032 #if defined(PETSC_USE_DEBUG) 1033 if (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); 1034 #endif 1035 for (j=0; j<n; j++) { 1036 if (indexn[j] < 0) { idx++; continue;} 1037 #if defined(PETSC_USE_DEBUG) 1038 if (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); 1039 #endif 1040 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 1041 } 1042 } 1043 } else { 1044 for (i=0; i<m; i++) { 1045 if (indexm[i] < 0) { idx += n; continue;} 1046 #if defined(PETSC_USE_DEBUG) 1047 if (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); 1048 #endif 1049 for (j=0; j<n; j++) { 1050 if (indexn[j] < 0) { idx++; continue;} 1051 #if defined(PETSC_USE_DEBUG) 1052 if (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); 1053 #endif 1054 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 1055 } 1056 } 1057 } 1058 } 1059 PetscFunctionReturn(0); 1060 } 1061 1062 static PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 1063 { 1064 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1065 PetscInt i,j; 1066 1067 PetscFunctionBegin; 1068 /* row-oriented output */ 1069 for (i=0; i<m; i++) { 1070 if (indexm[i] < 0) {v += n;continue;} 1071 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); 1072 for (j=0; j<n; j++) { 1073 if (indexn[j] < 0) {v++; continue;} 1074 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); 1075 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 1076 } 1077 } 1078 PetscFunctionReturn(0); 1079 } 1080 1081 /* -----------------------------------------------------------------*/ 1082 1083 static PetscErrorCode MatLoad_SeqDense_Binary(Mat newmat,PetscViewer viewer) 1084 { 1085 Mat_SeqDense *a; 1086 PetscErrorCode ierr; 1087 PetscInt *scols,i,j,nz,header[4]; 1088 int fd; 1089 PetscMPIInt size; 1090 PetscInt *rowlengths = 0,M,N,*cols,grows,gcols; 1091 PetscScalar *vals,*svals,*v,*w; 1092 MPI_Comm comm; 1093 1094 PetscFunctionBegin; 1095 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 1096 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1097 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 1098 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1099 ierr = PetscBinaryRead(fd,header,4,NULL,PETSC_INT);CHKERRQ(ierr); 1100 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 1101 M = header[1]; N = header[2]; nz = header[3]; 1102 1103 /* set global size if not set already*/ 1104 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 1105 ierr = MatSetSizes(newmat,M,N,M,N);CHKERRQ(ierr); 1106 } else { 1107 /* if sizes and type are already set, check if the vector global sizes are correct */ 1108 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 1109 if (M != grows || N != gcols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,grows,gcols); 1110 } 1111 a = (Mat_SeqDense*)newmat->data; 1112 if (!a->user_alloc) { 1113 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1114 } 1115 1116 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 1117 a = (Mat_SeqDense*)newmat->data; 1118 v = a->v; 1119 /* Allocate some temp space to read in the values and then flip them 1120 from row major to column major */ 1121 ierr = PetscMalloc1(M*N > 0 ? M*N : 1,&w);CHKERRQ(ierr); 1122 /* read in nonzero values */ 1123 ierr = PetscBinaryRead(fd,w,M*N,NULL,PETSC_SCALAR);CHKERRQ(ierr); 1124 /* now flip the values and store them in the matrix*/ 1125 for (j=0; j<N; j++) { 1126 for (i=0; i<M; i++) { 1127 *v++ =w[i*N+j]; 1128 } 1129 } 1130 ierr = PetscFree(w);CHKERRQ(ierr); 1131 } else { 1132 /* read row lengths */ 1133 ierr = PetscMalloc1(M+1,&rowlengths);CHKERRQ(ierr); 1134 ierr = PetscBinaryRead(fd,rowlengths,M,NULL,PETSC_INT);CHKERRQ(ierr); 1135 1136 a = (Mat_SeqDense*)newmat->data; 1137 v = a->v; 1138 1139 /* read column indices and nonzeros */ 1140 ierr = PetscMalloc1(nz+1,&scols);CHKERRQ(ierr); 1141 cols = scols; 1142 ierr = PetscBinaryRead(fd,cols,nz,NULL,PETSC_INT);CHKERRQ(ierr); 1143 ierr = PetscMalloc1(nz+1,&svals);CHKERRQ(ierr); 1144 vals = svals; 1145 ierr = PetscBinaryRead(fd,vals,nz,NULL,PETSC_SCALAR);CHKERRQ(ierr); 1146 1147 /* insert into matrix */ 1148 for (i=0; i<M; i++) { 1149 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 1150 svals += rowlengths[i]; scols += rowlengths[i]; 1151 } 1152 ierr = PetscFree(vals);CHKERRQ(ierr); 1153 ierr = PetscFree(cols);CHKERRQ(ierr); 1154 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 1155 } 1156 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1157 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1158 PetscFunctionReturn(0); 1159 } 1160 1161 PetscErrorCode MatLoad_SeqDense(Mat newMat, PetscViewer viewer) 1162 { 1163 PetscBool isbinary, ishdf5; 1164 PetscErrorCode ierr; 1165 1166 PetscFunctionBegin; 1167 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 1168 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 1169 /* force binary viewer to load .info file if it has not yet done so */ 1170 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 1171 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1172 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5);CHKERRQ(ierr); 1173 if (isbinary) { 1174 ierr = MatLoad_SeqDense_Binary(newMat,viewer);CHKERRQ(ierr); 1175 } else if (ishdf5) { 1176 #if defined(PETSC_HAVE_HDF5) 1177 ierr = MatLoad_Dense_HDF5(newMat,viewer);CHKERRQ(ierr); 1178 #else 1179 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 1180 #endif 1181 } else { 1182 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); 1183 } 1184 PetscFunctionReturn(0); 1185 } 1186 1187 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 1188 { 1189 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1190 PetscErrorCode ierr; 1191 PetscInt i,j; 1192 const char *name; 1193 PetscScalar *v; 1194 PetscViewerFormat format; 1195 #if defined(PETSC_USE_COMPLEX) 1196 PetscBool allreal = PETSC_TRUE; 1197 #endif 1198 1199 PetscFunctionBegin; 1200 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1201 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 1202 PetscFunctionReturn(0); /* do nothing for now */ 1203 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 1204 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1205 for (i=0; i<A->rmap->n; i++) { 1206 v = a->v + i; 1207 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 1208 for (j=0; j<A->cmap->n; j++) { 1209 #if defined(PETSC_USE_COMPLEX) 1210 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 1211 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1212 } else if (PetscRealPart(*v)) { 1213 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)PetscRealPart(*v));CHKERRQ(ierr); 1214 } 1215 #else 1216 if (*v) { 1217 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,(double)*v);CHKERRQ(ierr); 1218 } 1219 #endif 1220 v += a->lda; 1221 } 1222 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1223 } 1224 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1225 } else { 1226 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 1227 #if defined(PETSC_USE_COMPLEX) 1228 /* determine if matrix has all real values */ 1229 v = a->v; 1230 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 1231 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break;} 1232 } 1233 #endif 1234 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1235 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 1236 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1237 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 1239 } 1240 1241 for (i=0; i<A->rmap->n; i++) { 1242 v = a->v + i; 1243 for (j=0; j<A->cmap->n; j++) { 1244 #if defined(PETSC_USE_COMPLEX) 1245 if (allreal) { 1246 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));CHKERRQ(ierr); 1247 } else { 1248 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));CHKERRQ(ierr); 1249 } 1250 #else 1251 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);CHKERRQ(ierr); 1252 #endif 1253 v += a->lda; 1254 } 1255 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 1256 } 1257 if (format == PETSC_VIEWER_ASCII_MATLAB) { 1258 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 1259 } 1260 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 1261 } 1262 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1263 PetscFunctionReturn(0); 1264 } 1265 1266 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 1267 { 1268 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1269 PetscErrorCode ierr; 1270 int fd; 1271 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 1272 PetscScalar *v,*anonz,*vals; 1273 PetscViewerFormat format; 1274 1275 PetscFunctionBegin; 1276 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1277 1278 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1279 if (format == PETSC_VIEWER_NATIVE) { 1280 /* store the matrix as a dense matrix */ 1281 ierr = PetscMalloc1(4,&col_lens);CHKERRQ(ierr); 1282 1283 col_lens[0] = MAT_FILE_CLASSID; 1284 col_lens[1] = m; 1285 col_lens[2] = n; 1286 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 1287 1288 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1289 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1290 1291 /* write out matrix, by rows */ 1292 ierr = PetscMalloc1(m*n+1,&vals);CHKERRQ(ierr); 1293 v = a->v; 1294 for (j=0; j<n; j++) { 1295 for (i=0; i<m; i++) { 1296 vals[j + i*n] = *v++; 1297 } 1298 } 1299 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1300 ierr = PetscFree(vals);CHKERRQ(ierr); 1301 } else { 1302 ierr = PetscMalloc1(4+nz,&col_lens);CHKERRQ(ierr); 1303 1304 col_lens[0] = MAT_FILE_CLASSID; 1305 col_lens[1] = m; 1306 col_lens[2] = n; 1307 col_lens[3] = nz; 1308 1309 /* store lengths of each row and write (including header) to file */ 1310 for (i=0; i<m; i++) col_lens[4+i] = n; 1311 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 1312 1313 /* Possibly should write in smaller increments, not whole matrix at once? */ 1314 /* store column indices (zero start index) */ 1315 ict = 0; 1316 for (i=0; i<m; i++) { 1317 for (j=0; j<n; j++) col_lens[ict++] = j; 1318 } 1319 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 1320 ierr = PetscFree(col_lens);CHKERRQ(ierr); 1321 1322 /* store nonzero values */ 1323 ierr = PetscMalloc1(nz+1,&anonz);CHKERRQ(ierr); 1324 ict = 0; 1325 for (i=0; i<m; i++) { 1326 v = a->v + i; 1327 for (j=0; j<n; j++) { 1328 anonz[ict++] = *v; v += a->lda; 1329 } 1330 } 1331 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 1332 ierr = PetscFree(anonz);CHKERRQ(ierr); 1333 } 1334 PetscFunctionReturn(0); 1335 } 1336 1337 #include <petscdraw.h> 1338 static PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 1339 { 1340 Mat A = (Mat) Aa; 1341 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1342 PetscErrorCode ierr; 1343 PetscInt m = A->rmap->n,n = A->cmap->n,i,j; 1344 int color = PETSC_DRAW_WHITE; 1345 PetscScalar *v = a->v; 1346 PetscViewer viewer; 1347 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 1348 PetscViewerFormat format; 1349 1350 PetscFunctionBegin; 1351 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 1352 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 1353 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 1354 1355 /* Loop over matrix elements drawing boxes */ 1356 1357 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1358 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1359 /* Blue for negative and Red for positive */ 1360 for (j = 0; j < n; j++) { 1361 x_l = j; x_r = x_l + 1.0; 1362 for (i = 0; i < m; i++) { 1363 y_l = m - i - 1.0; 1364 y_r = y_l + 1.0; 1365 if (PetscRealPart(v[j*m+i]) > 0.) { 1366 color = PETSC_DRAW_RED; 1367 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1368 color = PETSC_DRAW_BLUE; 1369 } else { 1370 continue; 1371 } 1372 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1373 } 1374 } 1375 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1376 } else { 1377 /* use contour shading to indicate magnitude of values */ 1378 /* first determine max of all nonzero values */ 1379 PetscReal minv = 0.0, maxv = 0.0; 1380 PetscDraw popup; 1381 1382 for (i=0; i < m*n; i++) { 1383 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1384 } 1385 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1386 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1387 ierr = PetscDrawScalePopup(popup,minv,maxv);CHKERRQ(ierr); 1388 1389 ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr); 1390 for (j=0; j<n; j++) { 1391 x_l = j; 1392 x_r = x_l + 1.0; 1393 for (i=0; i<m; i++) { 1394 y_l = m - i - 1.0; 1395 y_r = y_l + 1.0; 1396 color = PetscDrawRealToColor(PetscAbsScalar(v[j*m+i]),minv,maxv); 1397 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1398 } 1399 } 1400 ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr); 1401 } 1402 PetscFunctionReturn(0); 1403 } 1404 1405 static PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1406 { 1407 PetscDraw draw; 1408 PetscBool isnull; 1409 PetscReal xr,yr,xl,yl,h,w; 1410 PetscErrorCode ierr; 1411 1412 PetscFunctionBegin; 1413 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1414 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1415 if (isnull) PetscFunctionReturn(0); 1416 1417 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1418 xr += w; yr += h; xl = -w; yl = -h; 1419 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1420 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1421 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1422 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL);CHKERRQ(ierr); 1423 ierr = PetscDrawSave(draw);CHKERRQ(ierr); 1424 PetscFunctionReturn(0); 1425 } 1426 1427 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1428 { 1429 PetscErrorCode ierr; 1430 PetscBool iascii,isbinary,isdraw; 1431 1432 PetscFunctionBegin; 1433 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1434 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1435 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1436 1437 if (iascii) { 1438 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1439 } else if (isbinary) { 1440 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1441 } else if (isdraw) { 1442 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1443 } 1444 PetscFunctionReturn(0); 1445 } 1446 1447 static PetscErrorCode MatDensePlaceArray_SeqDense(Mat A,const PetscScalar array[]) 1448 { 1449 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1450 1451 PetscFunctionBegin; 1452 a->unplacedarray = a->v; 1453 a->unplaced_user_alloc = a->user_alloc; 1454 a->v = (PetscScalar*) array; 1455 PetscFunctionReturn(0); 1456 } 1457 1458 static PetscErrorCode MatDenseResetArray_SeqDense(Mat A) 1459 { 1460 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1461 1462 PetscFunctionBegin; 1463 a->v = a->unplacedarray; 1464 a->user_alloc = a->unplaced_user_alloc; 1465 a->unplacedarray = NULL; 1466 PetscFunctionReturn(0); 1467 } 1468 1469 static PetscErrorCode MatDestroy_SeqDense(Mat mat) 1470 { 1471 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1472 PetscErrorCode ierr; 1473 1474 PetscFunctionBegin; 1475 #if defined(PETSC_USE_LOG) 1476 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1477 #endif 1478 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1479 ierr = PetscFree(l->fwork);CHKERRQ(ierr); 1480 ierr = MatDestroy(&l->ptapwork);CHKERRQ(ierr); 1481 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1482 ierr = PetscFree(mat->data);CHKERRQ(ierr); 1483 1484 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1485 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetLDA_C",NULL);CHKERRQ(ierr); 1486 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetArray_C",NULL);CHKERRQ(ierr); 1487 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDensePlaceArray_C",NULL);CHKERRQ(ierr); 1488 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseResetArray_C",NULL);CHKERRQ(ierr); 1489 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreArray_C",NULL);CHKERRQ(ierr); 1490 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_seqaij_C",NULL);CHKERRQ(ierr); 1491 #if defined(PETSC_HAVE_ELEMENTAL) 1492 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatConvert_seqdense_elemental_C",NULL);CHKERRQ(ierr); 1493 #endif 1494 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatSeqDenseSetPreallocation_C",NULL);CHKERRQ(ierr); 1495 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1496 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1497 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1498 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatPtAP_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1499 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMult_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1500 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1501 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatTransposeMatMultNumeric_seqaij_seqdense_C",NULL);CHKERRQ(ierr); 1502 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseGetColumn_C",NULL);CHKERRQ(ierr); 1503 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatDenseRestoreColumn_C",NULL);CHKERRQ(ierr); 1504 PetscFunctionReturn(0); 1505 } 1506 1507 static PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1508 { 1509 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1510 PetscErrorCode ierr; 1511 PetscInt k,j,m,n,M; 1512 PetscScalar *v,tmp; 1513 1514 PetscFunctionBegin; 1515 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1516 if (reuse == MAT_INPLACE_MATRIX && m == n) { /* in place transpose */ 1517 for (j=0; j<m; j++) { 1518 for (k=0; k<j; k++) { 1519 tmp = v[j + k*M]; 1520 v[j + k*M] = v[k + j*M]; 1521 v[k + j*M] = tmp; 1522 } 1523 } 1524 } else { /* out-of-place transpose */ 1525 Mat tmat; 1526 Mat_SeqDense *tmatd; 1527 PetscScalar *v2; 1528 PetscInt M2; 1529 1530 if (reuse != MAT_REUSE_MATRIX) { 1531 ierr = MatCreate(PetscObjectComm((PetscObject)A),&tmat);CHKERRQ(ierr); 1532 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1533 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1534 ierr = MatSeqDenseSetPreallocation(tmat,NULL);CHKERRQ(ierr); 1535 } else { 1536 tmat = *matout; 1537 } 1538 tmatd = (Mat_SeqDense*)tmat->data; 1539 v = mat->v; v2 = tmatd->v; M2 = tmatd->lda; 1540 for (j=0; j<n; j++) { 1541 for (k=0; k<m; k++) v2[j + k*M2] = v[k + j*M]; 1542 } 1543 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1544 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1545 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) *matout = tmat; 1546 else { 1547 ierr = MatHeaderMerge(A,&tmat);CHKERRQ(ierr); 1548 } 1549 } 1550 PetscFunctionReturn(0); 1551 } 1552 1553 static PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1554 { 1555 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1556 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1557 PetscInt i,j; 1558 PetscScalar *v1,*v2; 1559 1560 PetscFunctionBegin; 1561 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1562 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1563 for (i=0; i<A1->rmap->n; i++) { 1564 v1 = mat1->v+i; v2 = mat2->v+i; 1565 for (j=0; j<A1->cmap->n; j++) { 1566 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1567 v1 += mat1->lda; v2 += mat2->lda; 1568 } 1569 } 1570 *flg = PETSC_TRUE; 1571 PetscFunctionReturn(0); 1572 } 1573 1574 static PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1575 { 1576 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1577 PetscErrorCode ierr; 1578 PetscInt i,n,len; 1579 PetscScalar *x,zero = 0.0; 1580 1581 PetscFunctionBegin; 1582 ierr = VecSet(v,zero);CHKERRQ(ierr); 1583 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1584 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1585 len = PetscMin(A->rmap->n,A->cmap->n); 1586 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1587 for (i=0; i<len; i++) { 1588 x[i] = mat->v[i*mat->lda + i]; 1589 } 1590 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1591 PetscFunctionReturn(0); 1592 } 1593 1594 static PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1595 { 1596 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1597 const PetscScalar *l,*r; 1598 PetscScalar x,*v; 1599 PetscErrorCode ierr; 1600 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1601 1602 PetscFunctionBegin; 1603 if (ll) { 1604 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1605 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 1606 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1607 for (i=0; i<m; i++) { 1608 x = l[i]; 1609 v = mat->v + i; 1610 for (j=0; j<n; j++) { (*v) *= x; v+= mat->lda;} 1611 } 1612 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 1613 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1614 } 1615 if (rr) { 1616 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1617 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 1618 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1619 for (i=0; i<n; i++) { 1620 x = r[i]; 1621 v = mat->v + i*mat->lda; 1622 for (j=0; j<m; j++) (*v++) *= x; 1623 } 1624 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 1625 ierr = PetscLogFlops(1.0*n*m);CHKERRQ(ierr); 1626 } 1627 PetscFunctionReturn(0); 1628 } 1629 1630 static PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1631 { 1632 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1633 PetscScalar *v = mat->v; 1634 PetscReal sum = 0.0; 1635 PetscInt lda =mat->lda,m=A->rmap->n,i,j; 1636 PetscErrorCode ierr; 1637 1638 PetscFunctionBegin; 1639 if (type == NORM_FROBENIUS) { 1640 if (lda>m) { 1641 for (j=0; j<A->cmap->n; j++) { 1642 v = mat->v+j*lda; 1643 for (i=0; i<m; i++) { 1644 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1645 } 1646 } 1647 } else { 1648 #if defined(PETSC_USE_REAL___FP16) 1649 PetscBLASInt one = 1,cnt = A->cmap->n*A->rmap->n; 1650 *nrm = BLASnrm2_(&cnt,v,&one); 1651 } 1652 #else 1653 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1654 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1655 } 1656 } 1657 *nrm = PetscSqrtReal(sum); 1658 #endif 1659 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1660 } else if (type == NORM_1) { 1661 *nrm = 0.0; 1662 for (j=0; j<A->cmap->n; j++) { 1663 v = mat->v + j*mat->lda; 1664 sum = 0.0; 1665 for (i=0; i<A->rmap->n; i++) { 1666 sum += PetscAbsScalar(*v); v++; 1667 } 1668 if (sum > *nrm) *nrm = sum; 1669 } 1670 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1671 } else if (type == NORM_INFINITY) { 1672 *nrm = 0.0; 1673 for (j=0; j<A->rmap->n; j++) { 1674 v = mat->v + j; 1675 sum = 0.0; 1676 for (i=0; i<A->cmap->n; i++) { 1677 sum += PetscAbsScalar(*v); v += mat->lda; 1678 } 1679 if (sum > *nrm) *nrm = sum; 1680 } 1681 ierr = PetscLogFlops(1.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1682 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 static PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1687 { 1688 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1689 PetscErrorCode ierr; 1690 1691 PetscFunctionBegin; 1692 switch (op) { 1693 case MAT_ROW_ORIENTED: 1694 aij->roworiented = flg; 1695 break; 1696 case MAT_NEW_NONZERO_LOCATIONS: 1697 case MAT_NEW_NONZERO_LOCATION_ERR: 1698 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1699 case MAT_NEW_DIAGONALS: 1700 case MAT_KEEP_NONZERO_PATTERN: 1701 case MAT_IGNORE_OFF_PROC_ENTRIES: 1702 case MAT_USE_HASH_TABLE: 1703 case MAT_IGNORE_ZERO_ENTRIES: 1704 case MAT_IGNORE_LOWER_TRIANGULAR: 1705 case MAT_SORTED_FULL: 1706 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1707 break; 1708 case MAT_SPD: 1709 case MAT_SYMMETRIC: 1710 case MAT_STRUCTURALLY_SYMMETRIC: 1711 case MAT_HERMITIAN: 1712 case MAT_SYMMETRY_ETERNAL: 1713 /* These options are handled directly by MatSetOption() */ 1714 break; 1715 default: 1716 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1717 } 1718 PetscFunctionReturn(0); 1719 } 1720 1721 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1722 { 1723 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1724 PetscErrorCode ierr; 1725 PetscInt lda=l->lda,m=A->rmap->n,j; 1726 1727 PetscFunctionBegin; 1728 if (lda>m) { 1729 for (j=0; j<A->cmap->n; j++) { 1730 ierr = PetscArrayzero(l->v+j*lda,m);CHKERRQ(ierr); 1731 } 1732 } else { 1733 ierr = PetscArrayzero(l->v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1734 } 1735 PetscFunctionReturn(0); 1736 } 1737 1738 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1739 { 1740 PetscErrorCode ierr; 1741 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1742 PetscInt m = l->lda, n = A->cmap->n, i,j; 1743 PetscScalar *slot,*bb; 1744 const PetscScalar *xx; 1745 1746 PetscFunctionBegin; 1747 #if defined(PETSC_USE_DEBUG) 1748 for (i=0; i<N; i++) { 1749 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1750 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); 1751 } 1752 #endif 1753 1754 /* fix right hand side if needed */ 1755 if (x && b) { 1756 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1757 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1758 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1759 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1760 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1761 } 1762 1763 for (i=0; i<N; i++) { 1764 slot = l->v + rows[i]; 1765 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1766 } 1767 if (diag != 0.0) { 1768 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1769 for (i=0; i<N; i++) { 1770 slot = l->v + (m+1)*rows[i]; 1771 *slot = diag; 1772 } 1773 } 1774 PetscFunctionReturn(0); 1775 } 1776 1777 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1778 { 1779 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1780 1781 PetscFunctionBegin; 1782 *lda = mat->lda; 1783 PetscFunctionReturn(0); 1784 } 1785 1786 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1787 { 1788 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1789 1790 PetscFunctionBegin; 1791 if (mat->lda != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1792 *array = mat->v; 1793 PetscFunctionReturn(0); 1794 } 1795 1796 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1797 { 1798 PetscFunctionBegin; 1799 PetscFunctionReturn(0); 1800 } 1801 1802 /*@C 1803 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1804 1805 Logically Collective on Mat 1806 1807 Input Parameter: 1808 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1809 1810 Output Parameter: 1811 . lda - the leading dimension 1812 1813 Level: intermediate 1814 1815 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1816 @*/ 1817 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1818 { 1819 PetscErrorCode ierr; 1820 1821 PetscFunctionBegin; 1822 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1823 PetscFunctionReturn(0); 1824 } 1825 1826 /*@C 1827 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1828 1829 Logically Collective on Mat 1830 1831 Input Parameter: 1832 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1833 1834 Output Parameter: 1835 . array - pointer to the data 1836 1837 Level: intermediate 1838 1839 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1840 @*/ 1841 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1842 { 1843 PetscErrorCode ierr; 1844 1845 PetscFunctionBegin; 1846 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 /*@C 1851 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1852 1853 Logically Collective on Mat 1854 1855 Input Parameters: 1856 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1857 . array - pointer to the data 1858 1859 Level: intermediate 1860 1861 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1862 @*/ 1863 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1864 { 1865 PetscErrorCode ierr; 1866 1867 PetscFunctionBegin; 1868 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1869 if (array) *array = NULL; 1870 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1871 PetscFunctionReturn(0); 1872 } 1873 1874 /*@C 1875 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1876 1877 Not Collective 1878 1879 Input Parameter: 1880 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1881 1882 Output Parameter: 1883 . array - pointer to the data 1884 1885 Level: intermediate 1886 1887 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1888 @*/ 1889 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1890 { 1891 PetscErrorCode ierr; 1892 1893 PetscFunctionBegin; 1894 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1895 PetscFunctionReturn(0); 1896 } 1897 1898 /*@C 1899 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1900 1901 Not Collective 1902 1903 Input Parameters: 1904 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1905 . array - pointer to the data 1906 1907 Level: intermediate 1908 1909 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1910 @*/ 1911 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1912 { 1913 PetscErrorCode ierr; 1914 1915 PetscFunctionBegin; 1916 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1917 if (array) *array = NULL; 1918 PetscFunctionReturn(0); 1919 } 1920 1921 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1922 { 1923 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1924 PetscErrorCode ierr; 1925 PetscInt i,j,nrows,ncols; 1926 const PetscInt *irow,*icol; 1927 PetscScalar *av,*bv,*v = mat->v; 1928 Mat newmat; 1929 1930 PetscFunctionBegin; 1931 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1932 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1933 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1934 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1935 1936 /* Check submatrixcall */ 1937 if (scall == MAT_REUSE_MATRIX) { 1938 PetscInt n_cols,n_rows; 1939 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1940 if (n_rows != nrows || n_cols != ncols) { 1941 /* resize the result matrix to match number of requested rows/columns */ 1942 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1943 } 1944 newmat = *B; 1945 } else { 1946 /* Create and fill new matrix */ 1947 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1948 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1949 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1950 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1951 } 1952 1953 /* Now extract the data pointers and do the copy,column at a time */ 1954 bv = ((Mat_SeqDense*)newmat->data)->v; 1955 1956 for (i=0; i<ncols; i++) { 1957 av = v + mat->lda*icol[i]; 1958 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1959 } 1960 1961 /* Assemble the matrices so that the correct flags are set */ 1962 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1963 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1964 1965 /* Free work space */ 1966 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1967 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1968 *B = newmat; 1969 PetscFunctionReturn(0); 1970 } 1971 1972 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1973 { 1974 PetscErrorCode ierr; 1975 PetscInt i; 1976 1977 PetscFunctionBegin; 1978 if (scall == MAT_INITIAL_MATRIX) { 1979 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1980 } 1981 1982 for (i=0; i<n; i++) { 1983 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1984 } 1985 PetscFunctionReturn(0); 1986 } 1987 1988 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1989 { 1990 PetscFunctionBegin; 1991 PetscFunctionReturn(0); 1992 } 1993 1994 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1995 { 1996 PetscFunctionBegin; 1997 PetscFunctionReturn(0); 1998 } 1999 2000 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2001 { 2002 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2003 PetscErrorCode ierr; 2004 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2005 2006 PetscFunctionBegin; 2007 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2008 if (A->ops->copy != B->ops->copy) { 2009 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2010 PetscFunctionReturn(0); 2011 } 2012 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2013 if (lda1>m || lda2>m) { 2014 for (j=0; j<n; j++) { 2015 ierr = PetscArraycpy(b->v+j*lda2,a->v+j*lda1,m);CHKERRQ(ierr); 2016 } 2017 } else { 2018 ierr = PetscArraycpy(b->v,a->v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2019 } 2020 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2021 PetscFunctionReturn(0); 2022 } 2023 2024 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2025 { 2026 PetscErrorCode ierr; 2027 2028 PetscFunctionBegin; 2029 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2030 PetscFunctionReturn(0); 2031 } 2032 2033 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2034 { 2035 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2036 PetscInt i,nz = A->rmap->n*A->cmap->n; 2037 PetscScalar *aa = a->v; 2038 2039 PetscFunctionBegin; 2040 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2041 PetscFunctionReturn(0); 2042 } 2043 2044 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2045 { 2046 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2047 PetscInt i,nz = A->rmap->n*A->cmap->n; 2048 PetscScalar *aa = a->v; 2049 2050 PetscFunctionBegin; 2051 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2052 PetscFunctionReturn(0); 2053 } 2054 2055 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2056 { 2057 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2058 PetscInt i,nz = A->rmap->n*A->cmap->n; 2059 PetscScalar *aa = a->v; 2060 2061 PetscFunctionBegin; 2062 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 /* ----------------------------------------------------------------*/ 2067 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2068 { 2069 PetscErrorCode ierr; 2070 2071 PetscFunctionBegin; 2072 if (scall == MAT_INITIAL_MATRIX) { 2073 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2074 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2075 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2076 } 2077 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2078 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2079 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2080 PetscFunctionReturn(0); 2081 } 2082 2083 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2084 { 2085 PetscErrorCode ierr; 2086 PetscInt m=A->rmap->n,n=B->cmap->n; 2087 Mat Cmat; 2088 2089 PetscFunctionBegin; 2090 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 2091 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2092 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2093 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2094 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2095 2096 *C = Cmat; 2097 PetscFunctionReturn(0); 2098 } 2099 2100 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2101 { 2102 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2103 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2104 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2105 PetscBLASInt m,n,k; 2106 PetscScalar _DOne=1.0,_DZero=0.0; 2107 PetscErrorCode ierr; 2108 PetscBool flg; 2109 2110 PetscFunctionBegin; 2111 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 2112 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 2113 2114 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2115 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2116 if (flg) { 2117 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2118 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 2123 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 2124 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2125 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2126 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2127 if (!m || !n || !k) PetscFunctionReturn(0); 2128 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2129 PetscFunctionReturn(0); 2130 } 2131 2132 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2133 { 2134 PetscErrorCode ierr; 2135 2136 PetscFunctionBegin; 2137 if (scall == MAT_INITIAL_MATRIX) { 2138 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2139 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2140 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2141 } 2142 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2143 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2144 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2145 PetscFunctionReturn(0); 2146 } 2147 2148 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2149 { 2150 PetscErrorCode ierr; 2151 PetscInt m=A->rmap->n,n=B->rmap->n; 2152 Mat Cmat; 2153 2154 PetscFunctionBegin; 2155 if (A->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->cmap->n %d\n",A->cmap->n,B->cmap->n); 2156 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2157 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2158 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2159 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2160 2161 Cmat->assembled = PETSC_TRUE; 2162 2163 *C = Cmat; 2164 PetscFunctionReturn(0); 2165 } 2166 2167 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2168 { 2169 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2170 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2171 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2172 PetscBLASInt m,n,k; 2173 PetscScalar _DOne=1.0,_DZero=0.0; 2174 PetscErrorCode ierr; 2175 2176 PetscFunctionBegin; 2177 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2178 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2179 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2180 if (!m || !n || !k) PetscFunctionReturn(0); 2181 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2182 PetscFunctionReturn(0); 2183 } 2184 2185 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2186 { 2187 PetscErrorCode ierr; 2188 2189 PetscFunctionBegin; 2190 if (scall == MAT_INITIAL_MATRIX) { 2191 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2192 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2193 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2194 } 2195 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2196 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2197 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2198 PetscFunctionReturn(0); 2199 } 2200 2201 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2202 { 2203 PetscErrorCode ierr; 2204 PetscInt m=A->cmap->n,n=B->cmap->n; 2205 Mat Cmat; 2206 2207 PetscFunctionBegin; 2208 if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n); 2209 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2210 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2211 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2212 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2213 2214 Cmat->assembled = PETSC_TRUE; 2215 2216 *C = Cmat; 2217 PetscFunctionReturn(0); 2218 } 2219 2220 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2221 { 2222 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2223 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2224 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2225 PetscBLASInt m,n,k; 2226 PetscScalar _DOne=1.0,_DZero=0.0; 2227 PetscErrorCode ierr; 2228 2229 PetscFunctionBegin; 2230 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2231 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2232 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2233 if (!m || !n || !k) PetscFunctionReturn(0); 2234 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2235 PetscFunctionReturn(0); 2236 } 2237 2238 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2239 { 2240 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2241 PetscErrorCode ierr; 2242 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2243 PetscScalar *x; 2244 MatScalar *aa = a->v; 2245 2246 PetscFunctionBegin; 2247 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2248 2249 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2250 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2251 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2252 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2253 for (i=0; i<m; i++) { 2254 x[i] = aa[i]; if (idx) idx[i] = 0; 2255 for (j=1; j<n; j++) { 2256 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2257 } 2258 } 2259 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2260 PetscFunctionReturn(0); 2261 } 2262 2263 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2264 { 2265 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2266 PetscErrorCode ierr; 2267 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2268 PetscScalar *x; 2269 PetscReal atmp; 2270 MatScalar *aa = a->v; 2271 2272 PetscFunctionBegin; 2273 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2274 2275 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2276 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2277 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2278 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2279 for (i=0; i<m; i++) { 2280 x[i] = PetscAbsScalar(aa[i]); 2281 for (j=1; j<n; j++) { 2282 atmp = PetscAbsScalar(aa[i+m*j]); 2283 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2284 } 2285 } 2286 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2287 PetscFunctionReturn(0); 2288 } 2289 2290 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2291 { 2292 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2293 PetscErrorCode ierr; 2294 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2295 PetscScalar *x; 2296 MatScalar *aa = a->v; 2297 2298 PetscFunctionBegin; 2299 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2300 2301 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2302 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2303 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2304 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2305 for (i=0; i<m; i++) { 2306 x[i] = aa[i]; if (idx) idx[i] = 0; 2307 for (j=1; j<n; j++) { 2308 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2309 } 2310 } 2311 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 2315 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2316 { 2317 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2318 PetscErrorCode ierr; 2319 PetscScalar *x; 2320 2321 PetscFunctionBegin; 2322 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2323 2324 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2325 ierr = PetscArraycpy(x,a->v+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2326 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2327 PetscFunctionReturn(0); 2328 } 2329 2330 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2331 { 2332 PetscErrorCode ierr; 2333 PetscInt i,j,m,n; 2334 const PetscScalar *a; 2335 2336 PetscFunctionBegin; 2337 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2338 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2339 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2340 if (type == NORM_2) { 2341 for (i=0; i<n; i++) { 2342 for (j=0; j<m; j++) { 2343 norms[i] += PetscAbsScalar(a[j]*a[j]); 2344 } 2345 a += m; 2346 } 2347 } else if (type == NORM_1) { 2348 for (i=0; i<n; i++) { 2349 for (j=0; j<m; j++) { 2350 norms[i] += PetscAbsScalar(a[j]); 2351 } 2352 a += m; 2353 } 2354 } else if (type == NORM_INFINITY) { 2355 for (i=0; i<n; i++) { 2356 for (j=0; j<m; j++) { 2357 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2358 } 2359 a += m; 2360 } 2361 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2362 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2363 if (type == NORM_2) { 2364 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2365 } 2366 PetscFunctionReturn(0); 2367 } 2368 2369 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2370 { 2371 PetscErrorCode ierr; 2372 PetscScalar *a; 2373 PetscInt m,n,i; 2374 2375 PetscFunctionBegin; 2376 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2377 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2378 for (i=0; i<m*n; i++) { 2379 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2380 } 2381 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2382 PetscFunctionReturn(0); 2383 } 2384 2385 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2386 { 2387 PetscFunctionBegin; 2388 *missing = PETSC_FALSE; 2389 PetscFunctionReturn(0); 2390 } 2391 2392 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2393 { 2394 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2395 2396 PetscFunctionBegin; 2397 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2398 *vals = a->v+col*a->lda; 2399 PetscFunctionReturn(0); 2400 } 2401 2402 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2403 { 2404 PetscFunctionBegin; 2405 *vals = 0; /* user cannot accidently use the array later */ 2406 PetscFunctionReturn(0); 2407 } 2408 2409 /* -------------------------------------------------------------------*/ 2410 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2411 MatGetRow_SeqDense, 2412 MatRestoreRow_SeqDense, 2413 MatMult_SeqDense, 2414 /* 4*/ MatMultAdd_SeqDense, 2415 MatMultTranspose_SeqDense, 2416 MatMultTransposeAdd_SeqDense, 2417 0, 2418 0, 2419 0, 2420 /* 10*/ 0, 2421 MatLUFactor_SeqDense, 2422 MatCholeskyFactor_SeqDense, 2423 MatSOR_SeqDense, 2424 MatTranspose_SeqDense, 2425 /* 15*/ MatGetInfo_SeqDense, 2426 MatEqual_SeqDense, 2427 MatGetDiagonal_SeqDense, 2428 MatDiagonalScale_SeqDense, 2429 MatNorm_SeqDense, 2430 /* 20*/ MatAssemblyBegin_SeqDense, 2431 MatAssemblyEnd_SeqDense, 2432 MatSetOption_SeqDense, 2433 MatZeroEntries_SeqDense, 2434 /* 24*/ MatZeroRows_SeqDense, 2435 0, 2436 0, 2437 0, 2438 0, 2439 /* 29*/ MatSetUp_SeqDense, 2440 0, 2441 0, 2442 0, 2443 0, 2444 /* 34*/ MatDuplicate_SeqDense, 2445 0, 2446 0, 2447 0, 2448 0, 2449 /* 39*/ MatAXPY_SeqDense, 2450 MatCreateSubMatrices_SeqDense, 2451 0, 2452 MatGetValues_SeqDense, 2453 MatCopy_SeqDense, 2454 /* 44*/ MatGetRowMax_SeqDense, 2455 MatScale_SeqDense, 2456 MatShift_Basic, 2457 0, 2458 MatZeroRowsColumns_SeqDense, 2459 /* 49*/ MatSetRandom_SeqDense, 2460 0, 2461 0, 2462 0, 2463 0, 2464 /* 54*/ 0, 2465 0, 2466 0, 2467 0, 2468 0, 2469 /* 59*/ 0, 2470 MatDestroy_SeqDense, 2471 MatView_SeqDense, 2472 0, 2473 0, 2474 /* 64*/ 0, 2475 0, 2476 0, 2477 0, 2478 0, 2479 /* 69*/ MatGetRowMaxAbs_SeqDense, 2480 0, 2481 0, 2482 0, 2483 0, 2484 /* 74*/ 0, 2485 0, 2486 0, 2487 0, 2488 0, 2489 /* 79*/ 0, 2490 0, 2491 0, 2492 0, 2493 /* 83*/ MatLoad_SeqDense, 2494 0, 2495 MatIsHermitian_SeqDense, 2496 0, 2497 0, 2498 0, 2499 /* 89*/ MatMatMult_SeqDense_SeqDense, 2500 MatMatMultSymbolic_SeqDense_SeqDense, 2501 MatMatMultNumeric_SeqDense_SeqDense, 2502 MatPtAP_SeqDense_SeqDense, 2503 MatPtAPSymbolic_SeqDense_SeqDense, 2504 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2505 MatMatTransposeMult_SeqDense_SeqDense, 2506 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2507 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2508 0, 2509 /* 99*/ 0, 2510 0, 2511 0, 2512 MatConjugate_SeqDense, 2513 0, 2514 /*104*/ 0, 2515 MatRealPart_SeqDense, 2516 MatImaginaryPart_SeqDense, 2517 0, 2518 0, 2519 /*109*/ 0, 2520 0, 2521 MatGetRowMin_SeqDense, 2522 MatGetColumnVector_SeqDense, 2523 MatMissingDiagonal_SeqDense, 2524 /*114*/ 0, 2525 0, 2526 0, 2527 0, 2528 0, 2529 /*119*/ 0, 2530 0, 2531 0, 2532 0, 2533 0, 2534 /*124*/ 0, 2535 MatGetColumnNorms_SeqDense, 2536 0, 2537 0, 2538 0, 2539 /*129*/ 0, 2540 MatTransposeMatMult_SeqDense_SeqDense, 2541 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2542 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2543 0, 2544 /*134*/ 0, 2545 0, 2546 0, 2547 0, 2548 0, 2549 /*139*/ 0, 2550 0, 2551 0, 2552 0, 2553 0, 2554 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2555 }; 2556 2557 /*@C 2558 MatCreateSeqDense - Creates a sequential dense matrix that 2559 is stored in column major order (the usual Fortran 77 manner). Many 2560 of the matrix operations use the BLAS and LAPACK routines. 2561 2562 Collective 2563 2564 Input Parameters: 2565 + comm - MPI communicator, set to PETSC_COMM_SELF 2566 . m - number of rows 2567 . n - number of columns 2568 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2569 to control all matrix memory allocation. 2570 2571 Output Parameter: 2572 . A - the matrix 2573 2574 Notes: 2575 The data input variable is intended primarily for Fortran programmers 2576 who wish to allocate their own matrix memory space. Most users should 2577 set data=NULL. 2578 2579 Level: intermediate 2580 2581 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2582 @*/ 2583 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2584 { 2585 PetscErrorCode ierr; 2586 2587 PetscFunctionBegin; 2588 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2589 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2590 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2591 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2592 PetscFunctionReturn(0); 2593 } 2594 2595 /*@C 2596 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2597 2598 Collective 2599 2600 Input Parameters: 2601 + B - the matrix 2602 - data - the array (or NULL) 2603 2604 Notes: 2605 The data input variable is intended primarily for Fortran programmers 2606 who wish to allocate their own matrix memory space. Most users should 2607 need not call this routine. 2608 2609 Level: intermediate 2610 2611 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2612 2613 @*/ 2614 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2615 { 2616 PetscErrorCode ierr; 2617 2618 PetscFunctionBegin; 2619 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 2623 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2624 { 2625 Mat_SeqDense *b; 2626 PetscErrorCode ierr; 2627 2628 PetscFunctionBegin; 2629 B->preallocated = PETSC_TRUE; 2630 2631 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2632 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2633 2634 b = (Mat_SeqDense*)B->data; 2635 b->Mmax = B->rmap->n; 2636 b->Nmax = B->cmap->n; 2637 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2638 2639 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2640 if (!data) { /* petsc-allocated storage */ 2641 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2642 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2643 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2644 2645 b->user_alloc = PETSC_FALSE; 2646 } else { /* user-allocated storage */ 2647 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2648 b->v = data; 2649 b->user_alloc = PETSC_TRUE; 2650 } 2651 B->assembled = PETSC_TRUE; 2652 PetscFunctionReturn(0); 2653 } 2654 2655 #if defined(PETSC_HAVE_ELEMENTAL) 2656 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2657 { 2658 Mat mat_elemental; 2659 PetscErrorCode ierr; 2660 const PetscScalar *array; 2661 PetscScalar *v_colwise; 2662 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2663 2664 PetscFunctionBegin; 2665 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2666 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2667 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2668 k = 0; 2669 for (j=0; j<N; j++) { 2670 cols[j] = j; 2671 for (i=0; i<M; i++) { 2672 v_colwise[j*M+i] = array[k++]; 2673 } 2674 } 2675 for (i=0; i<M; i++) { 2676 rows[i] = i; 2677 } 2678 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2679 2680 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2681 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2682 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2683 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2684 2685 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2686 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2687 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2688 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2689 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2690 2691 if (reuse == MAT_INPLACE_MATRIX) { 2692 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2693 } else { 2694 *newmat = mat_elemental; 2695 } 2696 PetscFunctionReturn(0); 2697 } 2698 #endif 2699 2700 /*@C 2701 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2702 2703 Input parameter: 2704 + A - the matrix 2705 - lda - the leading dimension 2706 2707 Notes: 2708 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2709 it asserts that the preallocation has a leading dimension (the LDA parameter 2710 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2711 2712 Level: intermediate 2713 2714 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2715 2716 @*/ 2717 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2718 { 2719 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2720 2721 PetscFunctionBegin; 2722 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); 2723 b->lda = lda; 2724 b->changelda = PETSC_FALSE; 2725 b->Mmax = PetscMax(b->Mmax,lda); 2726 PetscFunctionReturn(0); 2727 } 2728 2729 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2730 { 2731 PetscErrorCode ierr; 2732 PetscMPIInt size; 2733 2734 PetscFunctionBegin; 2735 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2736 if (size == 1) { 2737 if (scall == MAT_INITIAL_MATRIX) { 2738 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2739 } else { 2740 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2741 } 2742 } else { 2743 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2744 } 2745 PetscFunctionReturn(0); 2746 } 2747 2748 /*MC 2749 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2750 2751 Options Database Keys: 2752 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2753 2754 Level: beginner 2755 2756 .seealso: MatCreateSeqDense() 2757 2758 M*/ 2759 2760 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2761 { 2762 Mat_SeqDense *b; 2763 PetscErrorCode ierr; 2764 PetscMPIInt size; 2765 2766 PetscFunctionBegin; 2767 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2768 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2769 2770 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2771 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2772 B->data = (void*)b; 2773 2774 b->roworiented = PETSC_TRUE; 2775 2776 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2777 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2778 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2779 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2780 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2781 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2782 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2783 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2784 #if defined(PETSC_HAVE_ELEMENTAL) 2785 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2786 #endif 2787 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2788 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2789 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2790 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2791 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2792 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2793 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2794 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2795 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2796 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2797 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2798 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2799 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2800 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2801 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2802 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2803 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2804 2805 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2806 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2807 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2808 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2809 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2810 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2811 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2812 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2813 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2814 2815 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2816 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2817 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2818 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2819 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2820 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2821 PetscFunctionReturn(0); 2822 } 2823 2824 /*@C 2825 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. 2826 2827 Not Collective 2828 2829 Input Parameter: 2830 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2831 - col - column index 2832 2833 Output Parameter: 2834 . vals - pointer to the data 2835 2836 Level: intermediate 2837 2838 .seealso: MatDenseRestoreColumn() 2839 @*/ 2840 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2841 { 2842 PetscErrorCode ierr; 2843 2844 PetscFunctionBegin; 2845 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2846 PetscFunctionReturn(0); 2847 } 2848 2849 /*@C 2850 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2851 2852 Not Collective 2853 2854 Input Parameter: 2855 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2856 2857 Output Parameter: 2858 . vals - pointer to the data 2859 2860 Level: intermediate 2861 2862 .seealso: MatDenseGetColumn() 2863 @*/ 2864 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2865 { 2866 PetscErrorCode ierr; 2867 2868 PetscFunctionBegin; 2869 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2870 PetscFunctionReturn(0); 2871 } 2872