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