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