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