1 2 /* 3 Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1 4 5 When build with PETSC_USE_64BIT_INDICES this will use UF_Long as the 6 integer type in UMFPACK, otherwise it will use int. This means 7 all integers in this file as simply declared as PetscInt. Also it means 8 that UMFPACK UL_Long version MUST be built with 64 bit integers when used. 9 10 */ 11 12 #include <../src/mat/impls/sbaij/seq/sbaij.h> 13 #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h> 14 15 /* 16 This is a terrible hack, but it allows the error handler to retain a context. 17 Note that this hack really cannot be made both reentrant and concurrent. 18 */ 19 static Mat static_F; 20 21 #undef __FUNCT__ 22 #define __FUNCT__ "CholmodErrorHandler" 23 static void CholmodErrorHandler(int status,const char *file,int line,const char *message) 24 { 25 26 PetscFunctionBegin; 27 if (status > CHOLMOD_OK) { 28 PetscInfo4(static_F,"CHOLMOD warning %d at %s:%d: %s",status,file,line,message); 29 } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */ 30 PetscInfo3(static_F,"CHOLMOD OK at %s:%d: %s",file,line,message); 31 } else { 32 PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message); 33 } 34 PetscFunctionReturnVoid(); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "CholmodStart" 39 PetscErrorCode CholmodStart(Mat F) 40 { 41 PetscErrorCode ierr; 42 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->spptr; 43 cholmod_common *c; 44 PetscBool flg; 45 46 PetscFunctionBegin; 47 if (chol->common) PetscFunctionReturn(0); 48 ierr = PetscMalloc1(1,&chol->common);CHKERRQ(ierr); 49 ierr = !cholmod_X_start(chol->common);CHKERRQ(ierr); 50 51 c = chol->common; 52 c->error_handler = CholmodErrorHandler; 53 54 #define CHOLMOD_OPTION_DOUBLE(name,help) do { \ 55 PetscReal tmp = (PetscReal)c->name; \ 56 ierr = PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 57 c->name = (double)tmp; \ 58 } while (0) 59 60 #define CHOLMOD_OPTION_INT(name,help) do { \ 61 PetscInt tmp = (PetscInt)c->name; \ 62 ierr = PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 63 c->name = (int)tmp; \ 64 } while (0) 65 66 #define CHOLMOD_OPTION_SIZE_T(name,help) do { \ 67 PetscInt tmp = (PetscInt)c->name; \ 68 ierr = PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 69 if (tmp < 0) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"value must be positive"); \ 70 c->name = (size_t)tmp; \ 71 } while (0) 72 73 #define CHOLMOD_OPTION_TRUTH(name,help) do { \ 74 PetscBool tmp = (PetscBool) !!c->name; \ 75 ierr = PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 76 c->name = (int)tmp; \ 77 } while (0) 78 79 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat");CHKERRQ(ierr); 80 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 81 chol->pack = (PetscBool)c->final_pack; 82 83 ierr = PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL);CHKERRQ(ierr); 84 c->final_pack = (int)chol->pack; 85 86 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 87 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 88 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 89 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 90 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 91 { 92 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 93 ierr = PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL);CHKERRQ(ierr); 94 } 95 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 96 CHOLMOD_OPTION_TRUTH(final_asis,"Leave factors \"as is\""); 97 CHOLMOD_OPTION_TRUTH(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 98 if (!c->final_asis) { 99 CHOLMOD_OPTION_TRUTH(final_super,"Leave supernodal factors instead of converting to simplicial"); 100 CHOLMOD_OPTION_TRUTH(final_ll,"Turn LDL' factorization into LL'"); 101 CHOLMOD_OPTION_TRUTH(final_monotonic,"Ensure columns are monotonic when done"); 102 CHOLMOD_OPTION_TRUTH(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 103 } 104 { 105 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 106 PetscInt n = 3; 107 ierr = PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 108 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 109 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 110 } 111 { 112 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 113 ierr = PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 114 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 115 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 116 } 117 CHOLMOD_OPTION_TRUTH(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 118 CHOLMOD_OPTION_TRUTH(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 119 CHOLMOD_OPTION_INT(print,"Verbosity level"); 120 ierr = PetscOptionsEnd();CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "MatWrapCholmod_seqsbaij" 126 static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc) 127 { 128 Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)A->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = PetscMemzero(C,sizeof(*C));CHKERRQ(ierr); 133 /* CHOLMOD uses column alignment, SBAIJ stores the upper factor, so we pass it on as a lower factor, swapping the meaning of row and column */ 134 C->nrow = (size_t)A->cmap->n; 135 C->ncol = (size_t)A->rmap->n; 136 C->nzmax = (size_t)sbaij->maxnz; 137 C->p = sbaij->i; 138 C->i = sbaij->j; 139 C->x = sbaij->a; 140 C->stype = -1; 141 C->itype = CHOLMOD_INT_TYPE; 142 C->xtype = CHOLMOD_SCALAR_TYPE; 143 C->dtype = CHOLMOD_DOUBLE; 144 C->sorted = 1; 145 C->packed = 1; 146 *aijalloc = PETSC_FALSE; 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "VecWrapCholmod" 152 static PetscErrorCode VecWrapCholmod(Vec X,cholmod_dense *Y) 153 { 154 PetscErrorCode ierr; 155 PetscScalar *x; 156 PetscInt n; 157 158 PetscFunctionBegin; 159 ierr = PetscMemzero(Y,sizeof(*Y));CHKERRQ(ierr); 160 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 161 ierr = VecGetSize(X,&n);CHKERRQ(ierr); 162 163 Y->x = (double*)x; 164 Y->nrow = n; 165 Y->ncol = 1; 166 Y->nzmax = n; 167 Y->d = n; 168 Y->x = (double*)x; 169 Y->xtype = CHOLMOD_SCALAR_TYPE; 170 Y->dtype = CHOLMOD_DOUBLE; 171 PetscFunctionReturn(0); 172 } 173 174 #undef __FUNCT__ 175 #define __FUNCT__ "MatDestroy_CHOLMOD" 176 PetscErrorCode MatDestroy_CHOLMOD(Mat F) 177 { 178 PetscErrorCode ierr; 179 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->spptr; 180 181 PetscFunctionBegin; 182 if (chol) { 183 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 184 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 185 ierr = PetscFree(chol->common);CHKERRQ(ierr); 186 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 187 ierr = (*chol->Destroy)(F);CHKERRQ(ierr); 188 } 189 ierr = PetscFree(F->spptr);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 193 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 194 195 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "MatFactorInfo_CHOLMOD" 199 static PetscErrorCode MatFactorInfo_CHOLMOD(Mat F,PetscViewer viewer) 200 { 201 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 202 const cholmod_common *c = chol->common; 203 PetscErrorCode ierr; 204 PetscInt i; 205 206 PetscFunctionBegin; 207 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 208 ierr = PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");CHKERRQ(ierr); 209 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");CHKERRQ(ierr); 211 ierr = PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound);CHKERRQ(ierr); 212 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0);CHKERRQ(ierr); 213 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1);CHKERRQ(ierr); 214 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2);CHKERRQ(ierr); 215 ierr = PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank);CHKERRQ(ierr); 216 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);CHKERRQ(ierr); 217 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal);CHKERRQ(ierr); 218 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis);CHKERRQ(ierr); 219 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super);CHKERRQ(ierr); 220 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll);CHKERRQ(ierr); 221 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack);CHKERRQ(ierr); 222 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic);CHKERRQ(ierr); 223 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol);CHKERRQ(ierr); 224 ierr = PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);CHKERRQ(ierr); 225 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);CHKERRQ(ierr); 226 ierr = PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper);CHKERRQ(ierr); 227 ierr = PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print);CHKERRQ(ierr); 228 for (i=0; i<c->nmethods; i++) { 229 ierr = PetscViewerASCIIPrintf(viewer,"Ordering method %D%s:\n",i,i==c->selected ? " [SELECTED]" : "");CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 231 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2);CHKERRQ(ierr); 232 } 233 ierr = PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder);CHKERRQ(ierr); 234 ierr = PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);CHKERRQ(ierr); 235 /* Statistics */ 236 ierr = PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl);CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz);CHKERRQ(ierr); 238 ierr = PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz);CHKERRQ(ierr); 239 ierr = PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl);CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count);CHKERRQ(ierr);CHKERRQ(ierr); 241 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage);CHKERRQ(ierr);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse);CHKERRQ(ierr);CHKERRQ(ierr); 243 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col);CHKERRQ(ierr);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);CHKERRQ(ierr);CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);CHKERRQ(ierr);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);CHKERRQ(ierr);CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);CHKERRQ(ierr);CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNCT__ 253 #define __FUNCT__ "MatView_CHOLMOD" 254 PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 255 { 256 PetscErrorCode ierr; 257 PetscBool iascii; 258 PetscViewerFormat format; 259 260 PetscFunctionBegin; 261 ierr = MatView_SeqSBAIJ(F,viewer);CHKERRQ(ierr); 262 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 263 if (iascii) { 264 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 265 if (format == PETSC_VIEWER_ASCII_INFO) { 266 ierr = MatFactorInfo_CHOLMOD(F,viewer);CHKERRQ(ierr); 267 } 268 } 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "MatSolve_CHOLMOD" 274 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 275 { 276 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 277 cholmod_dense cholB,*cholX; 278 PetscScalar *x; 279 PetscErrorCode ierr; 280 281 PetscFunctionBegin; 282 ierr = VecWrapCholmod(B,&cholB);CHKERRQ(ierr); 283 static_F = F; 284 cholX = cholmod_X_solve(CHOLMOD_A,chol->factor,&cholB,chol->common); 285 if (!cholX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CHOLMOD failed"); 286 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 287 ierr = PetscMemcpy(x,cholX->x,cholX->nrow*sizeof(*x));CHKERRQ(ierr); 288 ierr = !cholmod_X_free_dense(&cholX,chol->common);CHKERRQ(ierr); 289 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 #undef __FUNCT__ 294 #define __FUNCT__ "MatCholeskyFactorNumeric_CHOLMOD" 295 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 296 { 297 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 298 cholmod_sparse cholA; 299 PetscBool aijalloc; 300 PetscErrorCode ierr; 301 302 PetscFunctionBegin; 303 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc);CHKERRQ(ierr); 304 static_F = F; 305 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 306 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 307 if (chol->common->status == CHOLMOD_NOT_POSDEF) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_MAT_CH_ZRPVT,"CHOLMOD detected that the matrix is not positive definite, failure at column %u",(unsigned)chol->factor->minor); 308 309 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 310 311 F->ops->solve = MatSolve_CHOLMOD; 312 F->ops->solvetranspose = MatSolve_CHOLMOD; 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "MatCholeskyFactorSymbolic_CHOLMOD" 318 PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 319 { 320 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 321 PetscErrorCode ierr; 322 cholmod_sparse cholA; 323 PetscBool aijalloc; 324 PetscInt *fset = 0; 325 size_t fsize = 0; 326 327 PetscFunctionBegin; 328 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);CHKERRQ(ierr); 329 static_F = F; 330 if (chol->factor) { 331 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 332 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 333 } else if (perm) { 334 const PetscInt *ip; 335 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 336 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 337 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 338 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 339 } else { 340 chol->factor = cholmod_X_analyze(&cholA,chol->common); 341 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 342 } 343 344 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 345 346 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 347 PetscFunctionReturn(0); 348 } 349 350 #undef __FUNCT__ 351 #define __FUNCT__ "MatFactorGetSolverPackage_seqsbaij_cholmod" 352 PetscErrorCode MatFactorGetSolverPackage_seqsbaij_cholmod(Mat A,const MatSolverPackage *type) 353 { 354 PetscFunctionBegin; 355 *type = MATSOLVERCHOLMOD; 356 PetscFunctionReturn(0); 357 } 358 359 /*MC 360 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 361 via the external package CHOLMOD. 362 363 ./configure --download-suitesparse to install PETSc to use CHOLMOD 364 365 Consult CHOLMOD documentation for more information about the Common parameters 366 which correspond to the options database keys below. 367 368 Options Database Keys: 369 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 370 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 371 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 372 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 373 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 374 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 375 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 376 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 377 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 378 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 379 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 380 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 381 - -mat_cholmod_print <3> - Verbosity level (None) 382 383 Level: beginner 384 385 .seealso: PCCHOLESKY, PCFactorSetMatSolverPackage(), MatSolverPackage 386 M*/ 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "MatGetFactor_seqsbaij_cholmod" 390 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 391 { 392 Mat B; 393 Mat_CHOLMOD *chol; 394 PetscErrorCode ierr; 395 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 396 397 PetscFunctionBegin; 398 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s", 399 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 400 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 401 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 402 /* Create the factorization matrix F */ 403 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 404 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 405 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 406 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 407 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 408 409 chol->Wrap = MatWrapCholmod_seqsbaij; 410 chol->Destroy = MatDestroy_SeqSBAIJ; 411 B->spptr = chol; 412 413 B->ops->view = MatView_CHOLMOD; 414 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 415 B->ops->destroy = MatDestroy_CHOLMOD; 416 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqsbaij_cholmod);CHKERRQ(ierr); 417 B->factortype = MAT_FACTOR_CHOLESKY; 418 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 419 B->preallocated = PETSC_TRUE; 420 421 ierr = CholmodStart(B);CHKERRQ(ierr); 422 *F = B; 423 PetscFunctionReturn(0); 424 } 425