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