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