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 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1706 break; 1707 case MAT_SPD: 1708 case MAT_SYMMETRIC: 1709 case MAT_STRUCTURALLY_SYMMETRIC: 1710 case MAT_HERMITIAN: 1711 case MAT_SYMMETRY_ETERNAL: 1712 /* These options are handled directly by MatSetOption() */ 1713 break; 1714 default: 1715 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1716 } 1717 PetscFunctionReturn(0); 1718 } 1719 1720 static PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1721 { 1722 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1723 PetscErrorCode ierr; 1724 PetscInt lda=l->lda,m=A->rmap->n,j; 1725 1726 PetscFunctionBegin; 1727 if (lda>m) { 1728 for (j=0; j<A->cmap->n; j++) { 1729 ierr = PetscArrayzero(l->v+j*lda,m);CHKERRQ(ierr); 1730 } 1731 } else { 1732 ierr = PetscArrayzero(l->v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 1733 } 1734 PetscFunctionReturn(0); 1735 } 1736 1737 static PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1738 { 1739 PetscErrorCode ierr; 1740 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1741 PetscInt m = l->lda, n = A->cmap->n, i,j; 1742 PetscScalar *slot,*bb; 1743 const PetscScalar *xx; 1744 1745 PetscFunctionBegin; 1746 #if defined(PETSC_USE_DEBUG) 1747 for (i=0; i<N; i++) { 1748 if (rows[i] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row requested to be zeroed"); 1749 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); 1750 } 1751 #endif 1752 1753 /* fix right hand side if needed */ 1754 if (x && b) { 1755 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1756 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1757 for (i=0; i<N; i++) bb[rows[i]] = diag*xx[rows[i]]; 1758 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1759 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1760 } 1761 1762 for (i=0; i<N; i++) { 1763 slot = l->v + rows[i]; 1764 for (j=0; j<n; j++) { *slot = 0.0; slot += m;} 1765 } 1766 if (diag != 0.0) { 1767 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only coded for square matrices"); 1768 for (i=0; i<N; i++) { 1769 slot = l->v + (m+1)*rows[i]; 1770 *slot = diag; 1771 } 1772 } 1773 PetscFunctionReturn(0); 1774 } 1775 1776 static PetscErrorCode MatDenseGetLDA_SeqDense(Mat A,PetscInt *lda) 1777 { 1778 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1779 1780 PetscFunctionBegin; 1781 *lda = mat->lda; 1782 PetscFunctionReturn(0); 1783 } 1784 1785 static PetscErrorCode MatDenseGetArray_SeqDense(Mat A,PetscScalar *array[]) 1786 { 1787 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1788 1789 PetscFunctionBegin; 1790 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"); 1791 *array = mat->v; 1792 PetscFunctionReturn(0); 1793 } 1794 1795 static PetscErrorCode MatDenseRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1796 { 1797 PetscFunctionBegin; 1798 PetscFunctionReturn(0); 1799 } 1800 1801 /*@C 1802 MatDenseGetLDA - gets the leading dimension of the array returned from MatDenseGetArray() 1803 1804 Logically Collective on Mat 1805 1806 Input Parameter: 1807 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1808 1809 Output Parameter: 1810 . lda - the leading dimension 1811 1812 Level: intermediate 1813 1814 .seealso: MatDenseGetArray(), MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead(), MatSeqDenseSetLDA() 1815 @*/ 1816 PetscErrorCode MatDenseGetLDA(Mat A,PetscInt *lda) 1817 { 1818 PetscErrorCode ierr; 1819 1820 PetscFunctionBegin; 1821 ierr = PetscUseMethod(A,"MatDenseGetLDA_C",(Mat,PetscInt*),(A,lda));CHKERRQ(ierr); 1822 PetscFunctionReturn(0); 1823 } 1824 1825 /*@C 1826 MatDenseGetArray - gives access to the array where the data for a SeqDense matrix is stored 1827 1828 Logically Collective on Mat 1829 1830 Input Parameter: 1831 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1832 1833 Output Parameter: 1834 . array - pointer to the data 1835 1836 Level: intermediate 1837 1838 .seealso: MatDenseRestoreArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1839 @*/ 1840 PetscErrorCode MatDenseGetArray(Mat A,PetscScalar **array) 1841 { 1842 PetscErrorCode ierr; 1843 1844 PetscFunctionBegin; 1845 ierr = PetscUseMethod(A,"MatDenseGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1846 PetscFunctionReturn(0); 1847 } 1848 1849 /*@C 1850 MatDenseRestoreArray - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1851 1852 Logically Collective on Mat 1853 1854 Input Parameters: 1855 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1856 . array - pointer to the data 1857 1858 Level: intermediate 1859 1860 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArrayRead() 1861 @*/ 1862 PetscErrorCode MatDenseRestoreArray(Mat A,PetscScalar **array) 1863 { 1864 PetscErrorCode ierr; 1865 1866 PetscFunctionBegin; 1867 ierr = PetscUseMethod(A,"MatDenseRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 1868 if (array) *array = NULL; 1869 ierr = PetscObjectStateIncrease((PetscObject)A);CHKERRQ(ierr); 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /*@C 1874 MatDenseGetArrayRead - gives access to the array where the data for a SeqDense matrix is stored 1875 1876 Not Collective 1877 1878 Input Parameter: 1879 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1880 1881 Output Parameter: 1882 . array - pointer to the data 1883 1884 Level: intermediate 1885 1886 .seealso: MatDenseRestoreArray(), MatDenseGetArray(), MatDenseRestoreArrayRead() 1887 @*/ 1888 PetscErrorCode MatDenseGetArrayRead(Mat A,const PetscScalar **array) 1889 { 1890 PetscErrorCode ierr; 1891 1892 PetscFunctionBegin; 1893 ierr = PetscUseMethod(A,"MatDenseGetArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1894 PetscFunctionReturn(0); 1895 } 1896 1897 /*@C 1898 MatDenseRestoreArrayRead - returns access to the array where the data for a dense matrix is stored obtained by MatDenseGetArray() 1899 1900 Not Collective 1901 1902 Input Parameters: 1903 . mat - a MATSEQDENSE or MATMPIDENSE matrix 1904 . array - pointer to the data 1905 1906 Level: intermediate 1907 1908 .seealso: MatDenseGetArray(), MatDenseGetArrayRead(), MatDenseRestoreArray() 1909 @*/ 1910 PetscErrorCode MatDenseRestoreArrayRead(Mat A,const PetscScalar **array) 1911 { 1912 PetscErrorCode ierr; 1913 1914 PetscFunctionBegin; 1915 ierr = PetscUseMethod(A,"MatDenseRestoreArrayRead_C",(Mat,const PetscScalar**),(A,array));CHKERRQ(ierr); 1916 if (array) *array = NULL; 1917 PetscFunctionReturn(0); 1918 } 1919 1920 static PetscErrorCode MatCreateSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1921 { 1922 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1923 PetscErrorCode ierr; 1924 PetscInt i,j,nrows,ncols; 1925 const PetscInt *irow,*icol; 1926 PetscScalar *av,*bv,*v = mat->v; 1927 Mat newmat; 1928 1929 PetscFunctionBegin; 1930 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1931 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1932 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1933 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1934 1935 /* Check submatrixcall */ 1936 if (scall == MAT_REUSE_MATRIX) { 1937 PetscInt n_cols,n_rows; 1938 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1939 if (n_rows != nrows || n_cols != ncols) { 1940 /* resize the result matrix to match number of requested rows/columns */ 1941 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1942 } 1943 newmat = *B; 1944 } else { 1945 /* Create and fill new matrix */ 1946 ierr = MatCreate(PetscObjectComm((PetscObject)A),&newmat);CHKERRQ(ierr); 1947 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1948 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1949 ierr = MatSeqDenseSetPreallocation(newmat,NULL);CHKERRQ(ierr); 1950 } 1951 1952 /* Now extract the data pointers and do the copy,column at a time */ 1953 bv = ((Mat_SeqDense*)newmat->data)->v; 1954 1955 for (i=0; i<ncols; i++) { 1956 av = v + mat->lda*icol[i]; 1957 for (j=0; j<nrows; j++) *bv++ = av[irow[j]]; 1958 } 1959 1960 /* Assemble the matrices so that the correct flags are set */ 1961 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1962 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1963 1964 /* Free work space */ 1965 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1966 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1967 *B = newmat; 1968 PetscFunctionReturn(0); 1969 } 1970 1971 static PetscErrorCode MatCreateSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1972 { 1973 PetscErrorCode ierr; 1974 PetscInt i; 1975 1976 PetscFunctionBegin; 1977 if (scall == MAT_INITIAL_MATRIX) { 1978 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 1979 } 1980 1981 for (i=0; i<n; i++) { 1982 ierr = MatCreateSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1983 } 1984 PetscFunctionReturn(0); 1985 } 1986 1987 static PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1988 { 1989 PetscFunctionBegin; 1990 PetscFunctionReturn(0); 1991 } 1992 1993 static PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1994 { 1995 PetscFunctionBegin; 1996 PetscFunctionReturn(0); 1997 } 1998 1999 static PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 2000 { 2001 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense*)B->data; 2002 PetscErrorCode ierr; 2003 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 2004 2005 PetscFunctionBegin; 2006 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 2007 if (A->ops->copy != B->ops->copy) { 2008 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2009 PetscFunctionReturn(0); 2010 } 2011 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 2012 if (lda1>m || lda2>m) { 2013 for (j=0; j<n; j++) { 2014 ierr = PetscArraycpy(b->v+j*lda2,a->v+j*lda1,m);CHKERRQ(ierr); 2015 } 2016 } else { 2017 ierr = PetscArraycpy(b->v,a->v,A->rmap->n*A->cmap->n);CHKERRQ(ierr); 2018 } 2019 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2020 PetscFunctionReturn(0); 2021 } 2022 2023 static PetscErrorCode MatSetUp_SeqDense(Mat A) 2024 { 2025 PetscErrorCode ierr; 2026 2027 PetscFunctionBegin; 2028 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 2032 static PetscErrorCode MatConjugate_SeqDense(Mat A) 2033 { 2034 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2035 PetscInt i,nz = A->rmap->n*A->cmap->n; 2036 PetscScalar *aa = a->v; 2037 2038 PetscFunctionBegin; 2039 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 2040 PetscFunctionReturn(0); 2041 } 2042 2043 static PetscErrorCode MatRealPart_SeqDense(Mat A) 2044 { 2045 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2046 PetscInt i,nz = A->rmap->n*A->cmap->n; 2047 PetscScalar *aa = a->v; 2048 2049 PetscFunctionBegin; 2050 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 2055 { 2056 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2057 PetscInt i,nz = A->rmap->n*A->cmap->n; 2058 PetscScalar *aa = a->v; 2059 2060 PetscFunctionBegin; 2061 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 2062 PetscFunctionReturn(0); 2063 } 2064 2065 /* ----------------------------------------------------------------*/ 2066 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2067 { 2068 PetscErrorCode ierr; 2069 2070 PetscFunctionBegin; 2071 if (scall == MAT_INITIAL_MATRIX) { 2072 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2073 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2074 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2075 } 2076 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2077 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2078 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2079 PetscFunctionReturn(0); 2080 } 2081 2082 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2083 { 2084 PetscErrorCode ierr; 2085 PetscInt m=A->rmap->n,n=B->cmap->n; 2086 Mat Cmat; 2087 2088 PetscFunctionBegin; 2089 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); 2090 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2091 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2092 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2093 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2094 2095 *C = Cmat; 2096 PetscFunctionReturn(0); 2097 } 2098 2099 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2100 { 2101 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2102 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2103 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2104 PetscBLASInt m,n,k; 2105 PetscScalar _DOne=1.0,_DZero=0.0; 2106 PetscErrorCode ierr; 2107 PetscBool flg; 2108 2109 PetscFunctionBegin; 2110 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQDENSE,&flg);CHKERRQ(ierr); 2111 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Second matrix must be dense"); 2112 2113 /* Handle case where where user provided the final C matrix rather than calling MatMatMult() with MAT_INITIAL_MATRIX*/ 2114 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 2115 if (flg) { 2116 C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense; 2117 ierr = (*C->ops->matmultnumeric)(A,B,C);CHKERRQ(ierr); 2118 PetscFunctionReturn(0); 2119 } 2120 2121 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 2122 if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"First matrix must be dense"); 2123 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2124 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2125 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2126 if (!m || !n || !k) PetscFunctionReturn(0); 2127 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2128 PetscFunctionReturn(0); 2129 } 2130 2131 PetscErrorCode MatMatTransposeMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2132 { 2133 PetscErrorCode ierr; 2134 2135 PetscFunctionBegin; 2136 if (scall == MAT_INITIAL_MATRIX) { 2137 ierr = PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2138 ierr = MatMatTransposeMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2139 ierr = PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2140 } 2141 ierr = PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2142 ierr = MatMatTransposeMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2143 ierr = PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);CHKERRQ(ierr); 2144 PetscFunctionReturn(0); 2145 } 2146 2147 PetscErrorCode MatMatTransposeMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2148 { 2149 PetscErrorCode ierr; 2150 PetscInt m=A->rmap->n,n=B->rmap->n; 2151 Mat Cmat; 2152 2153 PetscFunctionBegin; 2154 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); 2155 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2156 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2157 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2158 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2159 2160 Cmat->assembled = PETSC_TRUE; 2161 2162 *C = Cmat; 2163 PetscFunctionReturn(0); 2164 } 2165 2166 PetscErrorCode MatMatTransposeMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2167 { 2168 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2169 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2170 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2171 PetscBLASInt m,n,k; 2172 PetscScalar _DOne=1.0,_DZero=0.0; 2173 PetscErrorCode ierr; 2174 2175 PetscFunctionBegin; 2176 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2177 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2178 ierr = PetscBLASIntCast(A->cmap->n,&k);CHKERRQ(ierr); 2179 if (!m || !n || !k) PetscFunctionReturn(0); 2180 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","T",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 PetscErrorCode MatTransposeMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 2185 { 2186 PetscErrorCode ierr; 2187 2188 PetscFunctionBegin; 2189 if (scall == MAT_INITIAL_MATRIX) { 2190 ierr = PetscLogEventBegin(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2191 ierr = MatTransposeMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 2192 ierr = PetscLogEventEnd(MAT_TransposeMatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 2193 } 2194 ierr = PetscLogEventBegin(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2195 ierr = MatTransposeMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 2196 ierr = PetscLogEventEnd(MAT_TransposeMatMultNumeric,A,B,0,0);CHKERRQ(ierr); 2197 PetscFunctionReturn(0); 2198 } 2199 2200 PetscErrorCode MatTransposeMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 2201 { 2202 PetscErrorCode ierr; 2203 PetscInt m=A->cmap->n,n=B->cmap->n; 2204 Mat Cmat; 2205 2206 PetscFunctionBegin; 2207 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); 2208 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 2209 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 2210 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 2211 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 2212 2213 Cmat->assembled = PETSC_TRUE; 2214 2215 *C = Cmat; 2216 PetscFunctionReturn(0); 2217 } 2218 2219 PetscErrorCode MatTransposeMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 2220 { 2221 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2222 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2223 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 2224 PetscBLASInt m,n,k; 2225 PetscScalar _DOne=1.0,_DZero=0.0; 2226 PetscErrorCode ierr; 2227 2228 PetscFunctionBegin; 2229 ierr = PetscBLASIntCast(C->rmap->n,&m);CHKERRQ(ierr); 2230 ierr = PetscBLASIntCast(C->cmap->n,&n);CHKERRQ(ierr); 2231 ierr = PetscBLASIntCast(A->rmap->n,&k);CHKERRQ(ierr); 2232 if (!m || !n || !k) PetscFunctionReturn(0); 2233 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda)); 2234 PetscFunctionReturn(0); 2235 } 2236 2237 static PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 2238 { 2239 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2240 PetscErrorCode ierr; 2241 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2242 PetscScalar *x; 2243 MatScalar *aa = a->v; 2244 2245 PetscFunctionBegin; 2246 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2247 2248 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2249 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2250 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2251 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2252 for (i=0; i<m; i++) { 2253 x[i] = aa[i]; if (idx) idx[i] = 0; 2254 for (j=1; j<n; j++) { 2255 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2256 } 2257 } 2258 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2259 PetscFunctionReturn(0); 2260 } 2261 2262 static PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 2263 { 2264 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2265 PetscErrorCode ierr; 2266 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2267 PetscScalar *x; 2268 PetscReal atmp; 2269 MatScalar *aa = a->v; 2270 2271 PetscFunctionBegin; 2272 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2273 2274 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2275 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2276 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2277 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2278 for (i=0; i<m; i++) { 2279 x[i] = PetscAbsScalar(aa[i]); 2280 for (j=1; j<n; j++) { 2281 atmp = PetscAbsScalar(aa[i+m*j]); 2282 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 2283 } 2284 } 2285 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2286 PetscFunctionReturn(0); 2287 } 2288 2289 static PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 2290 { 2291 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2292 PetscErrorCode ierr; 2293 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 2294 PetscScalar *x; 2295 MatScalar *aa = a->v; 2296 2297 PetscFunctionBegin; 2298 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2299 2300 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2301 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2302 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 2303 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2304 for (i=0; i<m; i++) { 2305 x[i] = aa[i]; if (idx) idx[i] = 0; 2306 for (j=1; j<n; j++) { 2307 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 2308 } 2309 } 2310 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 2314 static PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 2315 { 2316 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2317 PetscErrorCode ierr; 2318 PetscScalar *x; 2319 2320 PetscFunctionBegin; 2321 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2322 2323 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2324 ierr = PetscArraycpy(x,a->v+col*a->lda,A->rmap->n);CHKERRQ(ierr); 2325 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2326 PetscFunctionReturn(0); 2327 } 2328 2329 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 2330 { 2331 PetscErrorCode ierr; 2332 PetscInt i,j,m,n; 2333 const PetscScalar *a; 2334 2335 PetscFunctionBegin; 2336 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 2337 ierr = PetscArrayzero(norms,n);CHKERRQ(ierr); 2338 ierr = MatDenseGetArrayRead(A,&a);CHKERRQ(ierr); 2339 if (type == NORM_2) { 2340 for (i=0; i<n; i++) { 2341 for (j=0; j<m; j++) { 2342 norms[i] += PetscAbsScalar(a[j]*a[j]); 2343 } 2344 a += m; 2345 } 2346 } else if (type == NORM_1) { 2347 for (i=0; i<n; i++) { 2348 for (j=0; j<m; j++) { 2349 norms[i] += PetscAbsScalar(a[j]); 2350 } 2351 a += m; 2352 } 2353 } else if (type == NORM_INFINITY) { 2354 for (i=0; i<n; i++) { 2355 for (j=0; j<m; j++) { 2356 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 2357 } 2358 a += m; 2359 } 2360 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 2361 ierr = MatDenseRestoreArrayRead(A,&a);CHKERRQ(ierr); 2362 if (type == NORM_2) { 2363 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 2364 } 2365 PetscFunctionReturn(0); 2366 } 2367 2368 static PetscErrorCode MatSetRandom_SeqDense(Mat x,PetscRandom rctx) 2369 { 2370 PetscErrorCode ierr; 2371 PetscScalar *a; 2372 PetscInt m,n,i; 2373 2374 PetscFunctionBegin; 2375 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 2376 ierr = MatDenseGetArray(x,&a);CHKERRQ(ierr); 2377 for (i=0; i<m*n; i++) { 2378 ierr = PetscRandomGetValue(rctx,a+i);CHKERRQ(ierr); 2379 } 2380 ierr = MatDenseRestoreArray(x,&a);CHKERRQ(ierr); 2381 PetscFunctionReturn(0); 2382 } 2383 2384 static PetscErrorCode MatMissingDiagonal_SeqDense(Mat A,PetscBool *missing,PetscInt *d) 2385 { 2386 PetscFunctionBegin; 2387 *missing = PETSC_FALSE; 2388 PetscFunctionReturn(0); 2389 } 2390 2391 static PetscErrorCode MatDenseGetColumn_SeqDense(Mat A,PetscInt col,PetscScalar **vals) 2392 { 2393 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 2394 2395 PetscFunctionBegin; 2396 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2397 *vals = a->v+col*a->lda; 2398 PetscFunctionReturn(0); 2399 } 2400 2401 static PetscErrorCode MatDenseRestoreColumn_SeqDense(Mat A,PetscScalar **vals) 2402 { 2403 PetscFunctionBegin; 2404 *vals = 0; /* user cannot accidently use the array later */ 2405 PetscFunctionReturn(0); 2406 } 2407 2408 /* -------------------------------------------------------------------*/ 2409 static struct _MatOps MatOps_Values = { MatSetValues_SeqDense, 2410 MatGetRow_SeqDense, 2411 MatRestoreRow_SeqDense, 2412 MatMult_SeqDense, 2413 /* 4*/ MatMultAdd_SeqDense, 2414 MatMultTranspose_SeqDense, 2415 MatMultTransposeAdd_SeqDense, 2416 0, 2417 0, 2418 0, 2419 /* 10*/ 0, 2420 MatLUFactor_SeqDense, 2421 MatCholeskyFactor_SeqDense, 2422 MatSOR_SeqDense, 2423 MatTranspose_SeqDense, 2424 /* 15*/ MatGetInfo_SeqDense, 2425 MatEqual_SeqDense, 2426 MatGetDiagonal_SeqDense, 2427 MatDiagonalScale_SeqDense, 2428 MatNorm_SeqDense, 2429 /* 20*/ MatAssemblyBegin_SeqDense, 2430 MatAssemblyEnd_SeqDense, 2431 MatSetOption_SeqDense, 2432 MatZeroEntries_SeqDense, 2433 /* 24*/ MatZeroRows_SeqDense, 2434 0, 2435 0, 2436 0, 2437 0, 2438 /* 29*/ MatSetUp_SeqDense, 2439 0, 2440 0, 2441 0, 2442 0, 2443 /* 34*/ MatDuplicate_SeqDense, 2444 0, 2445 0, 2446 0, 2447 0, 2448 /* 39*/ MatAXPY_SeqDense, 2449 MatCreateSubMatrices_SeqDense, 2450 0, 2451 MatGetValues_SeqDense, 2452 MatCopy_SeqDense, 2453 /* 44*/ MatGetRowMax_SeqDense, 2454 MatScale_SeqDense, 2455 MatShift_Basic, 2456 0, 2457 MatZeroRowsColumns_SeqDense, 2458 /* 49*/ MatSetRandom_SeqDense, 2459 0, 2460 0, 2461 0, 2462 0, 2463 /* 54*/ 0, 2464 0, 2465 0, 2466 0, 2467 0, 2468 /* 59*/ 0, 2469 MatDestroy_SeqDense, 2470 MatView_SeqDense, 2471 0, 2472 0, 2473 /* 64*/ 0, 2474 0, 2475 0, 2476 0, 2477 0, 2478 /* 69*/ MatGetRowMaxAbs_SeqDense, 2479 0, 2480 0, 2481 0, 2482 0, 2483 /* 74*/ 0, 2484 0, 2485 0, 2486 0, 2487 0, 2488 /* 79*/ 0, 2489 0, 2490 0, 2491 0, 2492 /* 83*/ MatLoad_SeqDense, 2493 0, 2494 MatIsHermitian_SeqDense, 2495 0, 2496 0, 2497 0, 2498 /* 89*/ MatMatMult_SeqDense_SeqDense, 2499 MatMatMultSymbolic_SeqDense_SeqDense, 2500 MatMatMultNumeric_SeqDense_SeqDense, 2501 MatPtAP_SeqDense_SeqDense, 2502 MatPtAPSymbolic_SeqDense_SeqDense, 2503 /* 94*/ MatPtAPNumeric_SeqDense_SeqDense, 2504 MatMatTransposeMult_SeqDense_SeqDense, 2505 MatMatTransposeMultSymbolic_SeqDense_SeqDense, 2506 MatMatTransposeMultNumeric_SeqDense_SeqDense, 2507 0, 2508 /* 99*/ 0, 2509 0, 2510 0, 2511 MatConjugate_SeqDense, 2512 0, 2513 /*104*/ 0, 2514 MatRealPart_SeqDense, 2515 MatImaginaryPart_SeqDense, 2516 0, 2517 0, 2518 /*109*/ 0, 2519 0, 2520 MatGetRowMin_SeqDense, 2521 MatGetColumnVector_SeqDense, 2522 MatMissingDiagonal_SeqDense, 2523 /*114*/ 0, 2524 0, 2525 0, 2526 0, 2527 0, 2528 /*119*/ 0, 2529 0, 2530 0, 2531 0, 2532 0, 2533 /*124*/ 0, 2534 MatGetColumnNorms_SeqDense, 2535 0, 2536 0, 2537 0, 2538 /*129*/ 0, 2539 MatTransposeMatMult_SeqDense_SeqDense, 2540 MatTransposeMatMultSymbolic_SeqDense_SeqDense, 2541 MatTransposeMatMultNumeric_SeqDense_SeqDense, 2542 0, 2543 /*134*/ 0, 2544 0, 2545 0, 2546 0, 2547 0, 2548 /*139*/ 0, 2549 0, 2550 0, 2551 0, 2552 0, 2553 /*144*/ MatCreateMPIMatConcatenateSeqMat_SeqDense 2554 }; 2555 2556 /*@C 2557 MatCreateSeqDense - Creates a sequential dense matrix that 2558 is stored in column major order (the usual Fortran 77 manner). Many 2559 of the matrix operations use the BLAS and LAPACK routines. 2560 2561 Collective 2562 2563 Input Parameters: 2564 + comm - MPI communicator, set to PETSC_COMM_SELF 2565 . m - number of rows 2566 . n - number of columns 2567 - data - optional location of matrix data in column major order. Set data=NULL for PETSc 2568 to control all matrix memory allocation. 2569 2570 Output Parameter: 2571 . A - the matrix 2572 2573 Notes: 2574 The data input variable is intended primarily for Fortran programmers 2575 who wish to allocate their own matrix memory space. Most users should 2576 set data=NULL. 2577 2578 Level: intermediate 2579 2580 .seealso: MatCreate(), MatCreateDense(), MatSetValues() 2581 @*/ 2582 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 2583 { 2584 PetscErrorCode ierr; 2585 2586 PetscFunctionBegin; 2587 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2588 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2589 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 2590 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 2591 PetscFunctionReturn(0); 2592 } 2593 2594 /*@C 2595 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 2596 2597 Collective 2598 2599 Input Parameters: 2600 + B - the matrix 2601 - data - the array (or NULL) 2602 2603 Notes: 2604 The data input variable is intended primarily for Fortran programmers 2605 who wish to allocate their own matrix memory space. Most users should 2606 need not call this routine. 2607 2608 Level: intermediate 2609 2610 .seealso: MatCreate(), MatCreateDense(), MatSetValues(), MatSeqDenseSetLDA() 2611 2612 @*/ 2613 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 2614 { 2615 PetscErrorCode ierr; 2616 2617 PetscFunctionBegin; 2618 ierr = PetscTryMethod(B,"MatSeqDenseSetPreallocation_C",(Mat,PetscScalar[]),(B,data));CHKERRQ(ierr); 2619 PetscFunctionReturn(0); 2620 } 2621 2622 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 2623 { 2624 Mat_SeqDense *b; 2625 PetscErrorCode ierr; 2626 2627 PetscFunctionBegin; 2628 B->preallocated = PETSC_TRUE; 2629 2630 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2631 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2632 2633 b = (Mat_SeqDense*)B->data; 2634 b->Mmax = B->rmap->n; 2635 b->Nmax = B->cmap->n; 2636 if (b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 2637 2638 ierr = PetscIntMultError(b->lda,b->Nmax,NULL);CHKERRQ(ierr); 2639 if (!data) { /* petsc-allocated storage */ 2640 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2641 ierr = PetscCalloc1((size_t)b->lda*b->Nmax,&b->v);CHKERRQ(ierr); 2642 ierr = PetscLogObjectMemory((PetscObject)B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2643 2644 b->user_alloc = PETSC_FALSE; 2645 } else { /* user-allocated storage */ 2646 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2647 b->v = data; 2648 b->user_alloc = PETSC_TRUE; 2649 } 2650 B->assembled = PETSC_TRUE; 2651 PetscFunctionReturn(0); 2652 } 2653 2654 #if defined(PETSC_HAVE_ELEMENTAL) 2655 PETSC_INTERN PetscErrorCode MatConvert_SeqDense_Elemental(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 2656 { 2657 Mat mat_elemental; 2658 PetscErrorCode ierr; 2659 const PetscScalar *array; 2660 PetscScalar *v_colwise; 2661 PetscInt M=A->rmap->N,N=A->cmap->N,i,j,k,*rows,*cols; 2662 2663 PetscFunctionBegin; 2664 ierr = PetscMalloc3(M*N,&v_colwise,M,&rows,N,&cols);CHKERRQ(ierr); 2665 ierr = MatDenseGetArrayRead(A,&array);CHKERRQ(ierr); 2666 /* convert column-wise array into row-wise v_colwise, see MatSetValues_Elemental() */ 2667 k = 0; 2668 for (j=0; j<N; j++) { 2669 cols[j] = j; 2670 for (i=0; i<M; i++) { 2671 v_colwise[j*M+i] = array[k++]; 2672 } 2673 } 2674 for (i=0; i<M; i++) { 2675 rows[i] = i; 2676 } 2677 ierr = MatDenseRestoreArrayRead(A,&array);CHKERRQ(ierr); 2678 2679 ierr = MatCreate(PetscObjectComm((PetscObject)A), &mat_elemental);CHKERRQ(ierr); 2680 ierr = MatSetSizes(mat_elemental,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2681 ierr = MatSetType(mat_elemental,MATELEMENTAL);CHKERRQ(ierr); 2682 ierr = MatSetUp(mat_elemental);CHKERRQ(ierr); 2683 2684 /* PETSc-Elemental interaface uses axpy for setting off-processor entries, only ADD_VALUES is allowed */ 2685 ierr = MatSetValues(mat_elemental,M,rows,N,cols,v_colwise,ADD_VALUES);CHKERRQ(ierr); 2686 ierr = MatAssemblyBegin(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2687 ierr = MatAssemblyEnd(mat_elemental, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2688 ierr = PetscFree3(v_colwise,rows,cols);CHKERRQ(ierr); 2689 2690 if (reuse == MAT_INPLACE_MATRIX) { 2691 ierr = MatHeaderReplace(A,&mat_elemental);CHKERRQ(ierr); 2692 } else { 2693 *newmat = mat_elemental; 2694 } 2695 PetscFunctionReturn(0); 2696 } 2697 #endif 2698 2699 /*@C 2700 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2701 2702 Input parameter: 2703 + A - the matrix 2704 - lda - the leading dimension 2705 2706 Notes: 2707 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2708 it asserts that the preallocation has a leading dimension (the LDA parameter 2709 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2710 2711 Level: intermediate 2712 2713 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2714 2715 @*/ 2716 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 2717 { 2718 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2719 2720 PetscFunctionBegin; 2721 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); 2722 b->lda = lda; 2723 b->changelda = PETSC_FALSE; 2724 b->Mmax = PetscMax(b->Mmax,lda); 2725 PetscFunctionReturn(0); 2726 } 2727 2728 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqDense(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 2729 { 2730 PetscErrorCode ierr; 2731 PetscMPIInt size; 2732 2733 PetscFunctionBegin; 2734 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2735 if (size == 1) { 2736 if (scall == MAT_INITIAL_MATRIX) { 2737 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 2738 } else { 2739 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2740 } 2741 } else { 2742 ierr = MatCreateMPIMatConcatenateSeqMat_MPIDense(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 2743 } 2744 PetscFunctionReturn(0); 2745 } 2746 2747 /*MC 2748 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2749 2750 Options Database Keys: 2751 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2752 2753 Level: beginner 2754 2755 .seealso: MatCreateSeqDense() 2756 2757 M*/ 2758 2759 PETSC_EXTERN PetscErrorCode MatCreate_SeqDense(Mat B) 2760 { 2761 Mat_SeqDense *b; 2762 PetscErrorCode ierr; 2763 PetscMPIInt size; 2764 2765 PetscFunctionBegin; 2766 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 2767 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2768 2769 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 2770 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2771 B->data = (void*)b; 2772 2773 b->roworiented = PETSC_TRUE; 2774 2775 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetLDA_C",MatDenseGetLDA_SeqDense);CHKERRQ(ierr); 2776 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArray_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2777 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArray_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2778 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDensePlaceArray_C",MatDensePlaceArray_SeqDense);CHKERRQ(ierr); 2779 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseResetArray_C",MatDenseResetArray_SeqDense);CHKERRQ(ierr); 2780 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetArrayRead_C",MatDenseGetArray_SeqDense);CHKERRQ(ierr); 2781 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreArrayRead_C",MatDenseRestoreArray_SeqDense);CHKERRQ(ierr); 2782 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_seqaij_C",MatConvert_SeqDense_SeqAIJ);CHKERRQ(ierr); 2783 #if defined(PETSC_HAVE_ELEMENTAL) 2784 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqdense_elemental_C",MatConvert_SeqDense_Elemental);CHKERRQ(ierr); 2785 #endif 2786 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2787 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaij_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2788 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2789 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2790 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2791 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijperm_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2792 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijperm_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2793 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijperm_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2794 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijperm_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2795 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijsell_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2796 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijsell_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2797 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijsell_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2798 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijsell_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2799 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqaijmkl_seqdense_C",MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2800 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqaijmkl_seqdense_C",MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2801 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqaijmkl_seqdense_C",MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2802 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaijmkl_seqdense_C",MatPtAP_SeqDense_SeqDense);CHKERRQ(ierr); 2803 2804 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaij_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2805 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaij_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2806 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaij_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2807 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijperm_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2808 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijperm_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2809 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijperm_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2810 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijsell_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2811 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijsell_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2812 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijsell_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2813 2814 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMult_seqaijmkl_seqdense_C",MatTransposeMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2815 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultSymbolic_seqaijmkl_seqdense_C",MatTransposeMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2816 ierr = PetscObjectComposeFunction((PetscObject)B,"MatTransposeMatMultNumeric_seqaijmkl_seqdense_C",MatTransposeMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2817 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseGetColumn_C",MatDenseGetColumn_SeqDense);CHKERRQ(ierr); 2818 ierr = PetscObjectComposeFunction((PetscObject)B,"MatDenseRestoreColumn_C",MatDenseRestoreColumn_SeqDense);CHKERRQ(ierr); 2819 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2820 PetscFunctionReturn(0); 2821 } 2822 2823 /*@C 2824 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. 2825 2826 Not Collective 2827 2828 Input Parameter: 2829 + mat - a MATSEQDENSE or MATMPIDENSE matrix 2830 - col - column index 2831 2832 Output Parameter: 2833 . vals - pointer to the data 2834 2835 Level: intermediate 2836 2837 .seealso: MatDenseRestoreColumn() 2838 @*/ 2839 PetscErrorCode MatDenseGetColumn(Mat A,PetscInt col,PetscScalar **vals) 2840 { 2841 PetscErrorCode ierr; 2842 2843 PetscFunctionBegin; 2844 ierr = PetscUseMethod(A,"MatDenseGetColumn_C",(Mat,PetscInt,PetscScalar**),(A,col,vals));CHKERRQ(ierr); 2845 PetscFunctionReturn(0); 2846 } 2847 2848 /*@C 2849 MatDenseRestoreColumn - returns access to a column of a dense matrix which is returned by MatDenseGetColumn(). 2850 2851 Not Collective 2852 2853 Input Parameter: 2854 . mat - a MATSEQDENSE or MATMPIDENSE matrix 2855 2856 Output Parameter: 2857 . vals - pointer to the data 2858 2859 Level: intermediate 2860 2861 .seealso: MatDenseGetColumn() 2862 @*/ 2863 PetscErrorCode MatDenseRestoreColumn(Mat A,PetscScalar **vals) 2864 { 2865 PetscErrorCode ierr; 2866 2867 PetscFunctionBegin; 2868 ierr = PetscUseMethod(A,"MatDenseRestoreColumn_C",(Mat,PetscScalar**),(A,vals));CHKERRQ(ierr); 2869 PetscFunctionReturn(0); 2870 } 2871