1 2 /* 3 Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1 4 5 When built 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 PetscReal tmp = (PetscInt)c->name; \ 65 ierr = PetscOptionsReal("-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_OPTION_INT(nmethods,"Number of different ordering methods to try"); 78 79 #if defined(PETSC_USE_SUITESPARSE_GPU) 80 c->useGPU = 1; 81 CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0"); 82 CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes,"Maximum memory to allocate on the GPU"); 83 CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction,"Fraction of available GPU memory to allocate"); 84 #endif 85 86 /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */ 87 chol->pack = (PetscBool)c->final_pack; 88 ierr = PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL);CHKERRQ(ierr); 89 c->final_pack = (int)chol->pack; 90 91 CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D"); 92 CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified"); 93 CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified"); 94 CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified"); 95 CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]"); 96 { 97 static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0}; 98 ierr = PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL);CHKERRQ(ierr); 99 } 100 if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization"); 101 CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\""); 102 CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)"); 103 if (!c->final_asis) { 104 CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial"); 105 CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'"); 106 CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done"); 107 CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation"); 108 } 109 { 110 PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]}; 111 PetscInt n = 3; 112 ierr = PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 113 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_zrelax"); 114 if (flg) while (n--) c->zrelax[n] = (double)tmp[n]; 115 } 116 { 117 PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]}; 118 ierr = PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);CHKERRQ(ierr); 119 if (flg && n != 3) SETERRQ(PetscObjectComm((PetscObject)F),PETSC_ERR_ARG_OUTOFRANGE,"must provide exactly 3 parameters to -mat_cholmod_nrelax"); 120 if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n]; 121 } 122 CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]"); 123 CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection"); 124 CHOLMOD_OPTION_INT(print,"Verbosity level"); 125 ierr = PetscOptionsEnd();CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 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 static PetscErrorCode VecWrapCholmodRead(Vec X,cholmod_dense *Y) 154 { 155 PetscErrorCode ierr; 156 const PetscScalar *x; 157 PetscInt n; 158 159 PetscFunctionBegin; 160 ierr = PetscMemzero(Y,sizeof(*Y));CHKERRQ(ierr); 161 ierr = VecGetArrayRead(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 static PetscErrorCode VecUnWrapCholmodRead(Vec X,cholmod_dense *Y) 176 { 177 PetscErrorCode ierr; 178 179 PetscFunctionBegin; 180 ierr = VecRestoreArrayRead(X,NULL);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F) 185 { 186 PetscErrorCode ierr; 187 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data; 188 189 PetscFunctionBegin; 190 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 191 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 192 ierr = PetscFree(chol->common);CHKERRQ(ierr); 193 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 194 ierr = PetscFree(F->data);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 199 200 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 201 202 static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer) 203 { 204 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 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 #if defined(PETSC_USE_SUITESPARSE_GPU) 252 ierr = PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);CHKERRQ(ierr);CHKERRQ(ierr); 253 #endif 254 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257 258 PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 259 { 260 PetscErrorCode ierr; 261 PetscBool iascii; 262 PetscViewerFormat format; 263 264 PetscFunctionBegin; 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 = MatView_Info_CHOLMOD(F,viewer);CHKERRQ(ierr); 270 } 271 } 272 PetscFunctionReturn(0); 273 } 274 275 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 276 { 277 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 278 cholmod_dense cholB,*cholX; 279 PetscScalar *x; 280 PetscErrorCode ierr; 281 282 PetscFunctionBegin; 283 ierr = VecWrapCholmodRead(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 = VecUnWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 288 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 289 ierr = PetscArraycpy(x,(PetscScalar*)cholX->x,cholX->nrow);CHKERRQ(ierr); 290 ierr = !cholmod_X_free_dense(&cholX,chol->common);CHKERRQ(ierr); 291 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 292 PetscFunctionReturn(0); 293 } 294 295 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 296 { 297 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 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 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 317 { 318 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 319 PetscErrorCode ierr; 320 cholmod_sparse cholA; 321 PetscBool aijalloc; 322 PetscInt *fset = 0; 323 size_t fsize = 0; 324 325 PetscFunctionBegin; 326 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);CHKERRQ(ierr); 327 static_F = F; 328 if (chol->factor) { 329 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 330 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 331 } else if (perm) { 332 const PetscInt *ip; 333 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 334 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 335 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 336 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 337 } else { 338 chol->factor = cholmod_X_analyze(&cholA,chol->common); 339 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 340 } 341 342 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 343 344 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 345 PetscFunctionReturn(0); 346 } 347 348 static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type) 349 { 350 PetscFunctionBegin; 351 *type = MATSOLVERCHOLMOD; 352 PetscFunctionReturn(0); 353 } 354 355 /*MC 356 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 357 via the external package CHOLMOD. 358 359 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 360 361 Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver 362 363 Consult CHOLMOD documentation for more information about the Common parameters 364 which correspond to the options database keys below. 365 366 Options Database Keys: 367 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 368 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 369 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 370 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 371 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 372 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 373 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 374 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 375 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 376 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 377 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 378 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 379 - -mat_cholmod_print <3> - Verbosity level (None) 380 381 Level: beginner 382 383 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 384 385 .seealso: PCCHOLESKY, PCFactorSetMatSolverType(), MatSolverType 386 M*/ 387 388 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 389 { 390 Mat B; 391 Mat_CHOLMOD *chol; 392 PetscErrorCode ierr; 393 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 394 395 PetscFunctionBegin; 396 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s", 397 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 398 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 399 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 400 /* Create the factorization matrix F */ 401 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 402 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 403 ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr); 404 ierr = MatSetUp(B);CHKERRQ(ierr); 405 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 406 407 chol->Wrap = MatWrapCholmod_seqsbaij; 408 B->data = chol; 409 410 B->ops->getinfo = MatGetInfo_External; 411 B->ops->view = MatView_CHOLMOD; 412 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 413 B->ops->destroy = MatDestroy_CHOLMOD; 414 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);CHKERRQ(ierr); 415 B->factortype = MAT_FACTOR_CHOLESKY; 416 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 417 B->preallocated = PETSC_TRUE; 418 419 ierr = CholmodStart(B);CHKERRQ(ierr); 420 421 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 422 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 423 424 *F = B; 425 PetscFunctionReturn(0); 426 } 427