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