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