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 Suitesparse_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 one cannot use 64BIT_INDICES on 32bit machines [as Suitesparse_long is 32bit only] 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 static void CholmodErrorHandler(int status,const char *file,int line,const char *message) 22 { 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 if (status > CHOLMOD_OK) { 27 ierr = PetscInfo4(static_F,"CHOLMOD warning %d at %s:%d: %s\n",status,file,line,message);CHKERRV(ierr); 28 } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */ 29 ierr = PetscInfo3(static_F,"CHOLMOD OK at %s:%d: %s\n",file,line,message);CHKERRV(ierr); 30 } else { 31 ierr = PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message);CHKERRV(ierr); 32 } 33 PetscFunctionReturnVoid(); 34 } 35 36 PetscErrorCode CholmodStart(Mat F) 37 { 38 PetscErrorCode ierr; 39 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 40 cholmod_common *c; 41 PetscBool flg; 42 43 PetscFunctionBegin; 44 if (chol->common) PetscFunctionReturn(0); 45 ierr = PetscMalloc1(1,&chol->common);CHKERRQ(ierr); 46 ierr = !cholmod_X_start(chol->common);CHKERRQ(ierr); 47 48 c = chol->common; 49 c->error_handler = CholmodErrorHandler; 50 51 #define CHOLMOD_OPTION_DOUBLE(name,help) do { \ 52 PetscReal tmp = (PetscReal)c->name; \ 53 ierr = PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 54 c->name = (double)tmp; \ 55 } while (0) 56 57 #define CHOLMOD_OPTION_INT(name,help) do { \ 58 PetscInt tmp = (PetscInt)c->name; \ 59 ierr = PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 60 c->name = (int)tmp; \ 61 } while (0) 62 63 #define CHOLMOD_OPTION_SIZE_T(name,help) do { \ 64 PetscInt tmp = (PetscInt)c->name; \ 65 ierr = PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 66 if (tmp < 0) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"value must be positive"); \ 67 c->name = (size_t)tmp; \ 68 } while (0) 69 70 #define CHOLMOD_OPTION_BOOL(name,help) do { \ 71 PetscBool tmp = (PetscBool) !!c->name; \ 72 ierr = PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL);CHKERRQ(ierr); \ 73 c->name = (int)tmp; \ 74 } while (0) 75 76 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat");CHKERRQ(ierr); 77 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 78 chol->pack = (PetscBool)c->final_pack; 79 80 #if defined(PETSC_USE_SUITESPARSE_GPU) 81 c->useGPU = 1; 82 CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0"); 83 #endif 84 85 ierr = PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL);CHKERRQ(ierr); 86 c->final_pack = (int)chol->pack; 87 88 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 89 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 90 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 91 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 92 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 93 { 94 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 95 ierr = PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL);CHKERRQ(ierr); 96 } 97 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 98 CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\""); 99 CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 100 if (!c->final_asis) { 101 CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial"); 102 CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'"); 103 CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done"); 104 CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 105 } 106 { 107 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 108 PetscInt n = 3; 109 ierr = PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 110 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 111 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 112 } 113 { 114 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 115 ierr = PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 116 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 117 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 118 } 119 CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 120 CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 121 CHOLMOD_OPTION_INT(print,"Verbosity level"); 122 ierr = PetscOptionsEnd();CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 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 static PetscErrorCode VecWrapCholmodRead(Vec X,cholmod_dense *Y) 151 { 152 PetscErrorCode ierr; 153 const PetscScalar *x; 154 PetscInt n; 155 156 PetscFunctionBegin; 157 ierr = PetscMemzero(Y,sizeof(*Y));CHKERRQ(ierr); 158 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 159 ierr = VecGetSize(X,&n);CHKERRQ(ierr); 160 161 Y->x = (double*)x; 162 Y->nrow = n; 163 Y->ncol = 1; 164 Y->nzmax = n; 165 Y->d = n; 166 Y->x = (double*)x; 167 Y->xtype = CHOLMOD_SCALAR_TYPE; 168 Y->dtype = CHOLMOD_DOUBLE; 169 PetscFunctionReturn(0); 170 } 171 172 static PetscErrorCode VecUnWrapCholmodRead(Vec X,cholmod_dense *Y) 173 { 174 PetscErrorCode ierr; 175 176 PetscFunctionBegin; 177 ierr = VecRestoreArrayRead(X,NULL);CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 182 { 183 PetscErrorCode ierr; 184 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 185 186 PetscFunctionBegin; 187 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 188 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 189 ierr = PetscFree(chol->common);CHKERRQ(ierr); 190 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 191 ierr = PetscFree(F->data);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 196 197 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 198 199 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 200 { 201 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 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 #if defined(PETSC_USE_SUITESPARSE_GPU) 249 ierr = PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);CHKERRQ(ierr);CHKERRQ(ierr); 250 #endif 251 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 256 { 257 PetscErrorCode ierr; 258 PetscBool iascii; 259 PetscViewerFormat format; 260 261 PetscFunctionBegin; 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 = MatView_Info_CHOLMOD(F,viewer);CHKERRQ(ierr); 267 } 268 } 269 PetscFunctionReturn(0); 270 } 271 272 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 273 { 274 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 275 cholmod_dense cholB,*cholX; 276 PetscScalar *x; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = VecWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 281 static_F = F; 282 cholX = cholmod_X_solve(CHOLMOD_A,chol->factor,&cholB,chol->common); 283 if (!cholX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CHOLMOD failed"); 284 ierr = VecUnWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 285 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 286 ierr = PetscMemcpy(x,cholX->x,cholX->nrow*sizeof(*x));CHKERRQ(ierr); 287 ierr = !cholmod_X_free_dense(&cholX,chol->common);CHKERRQ(ierr); 288 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 289 PetscFunctionReturn(0); 290 } 291 292 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 293 { 294 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 295 cholmod_sparse cholA; 296 PetscBool aijalloc; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc);CHKERRQ(ierr); 301 static_F = F; 302 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 303 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 304 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); 305 306 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 307 308 F->ops->solve = MatSolve_CHOLMOD; 309 F->ops->solvetranspose = MatSolve_CHOLMOD; 310 PetscFunctionReturn(0); 311 } 312 313 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 314 { 315 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 316 PetscErrorCode ierr; 317 cholmod_sparse cholA; 318 PetscBool aijalloc; 319 PetscInt *fset = 0; 320 size_t fsize = 0; 321 322 PetscFunctionBegin; 323 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);CHKERRQ(ierr); 324 static_F = F; 325 if (chol->factor) { 326 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 327 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 328 } else if (perm) { 329 const PetscInt *ip; 330 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 331 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 332 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 333 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 334 } else { 335 chol->factor = cholmod_X_analyze(&cholA,chol->common); 336 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 337 } 338 339 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 340 341 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 342 PetscFunctionReturn(0); 343 } 344 345 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 346 { 347 PetscFunctionBegin; 348 *type = MATSOLVERCHOLMOD; 349 PetscFunctionReturn(0); 350 } 351 352 /*MC 353 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 354 via the external package CHOLMOD. 355 356 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 357 358 Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver 359 360 Consult CHOLMOD documentation for more information about the Common parameters 361 which correspond to the options database keys below. 362 363 Options Database Keys: 364 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 365 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 366 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 367 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 368 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 369 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 370 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 371 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 372 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 373 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 374 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 375 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 376 - -mat_cholmod_print <3> - Verbosity level (None) 377 378 Level: beginner 379 380 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 381 382 .seealso: PCCHOLESKY, PCFactorSetMatSolverType(), MatSolverType 383 M*/ 384 385 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 386 { 387 Mat B; 388 Mat_CHOLMOD *chol; 389 PetscErrorCode ierr; 390 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 391 392 PetscFunctionBegin; 393 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s", 394 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 395 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 396 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 397 /* Create the factorization matrix F */ 398 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 399 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 400 ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr); 401 ierr = MatSetUp(B);CHKERRQ(ierr); 402 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 403 404 chol->Wrap = MatWrapCholmod_seqsbaij; 405 B->data = chol; 406 407 B->ops->getinfo = MatGetInfo_External; 408 B->ops->view = MatView_CHOLMOD; 409 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 410 B->ops->destroy = MatDestroy_CHOLMOD; 411 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);CHKERRQ(ierr); 412 B->factortype = MAT_FACTOR_CHOLESKY; 413 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 414 B->preallocated = PETSC_TRUE; 415 416 ierr = CholmodStart(B);CHKERRQ(ierr); 417 418 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 419 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 420 421 *F = B; 422 PetscFunctionReturn(0); 423 } 424