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 #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__ "VecWrapCholmodRead" 158 static PetscErrorCode VecWrapCholmodRead(Vec X,cholmod_dense *Y) 159 { 160 PetscErrorCode ierr; 161 const PetscScalar *x; 162 PetscInt n; 163 164 PetscFunctionBegin; 165 ierr = PetscMemzero(Y,sizeof(*Y));CHKERRQ(ierr); 166 ierr = VecGetArrayRead(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__ "VecUnWrapCholmodRead" 182 static PetscErrorCode VecUnWrapCholmodRead(Vec X,cholmod_dense *Y) 183 { 184 PetscErrorCode ierr; 185 186 PetscFunctionBegin; 187 ierr = VecRestoreArrayRead(X,PETSC_NULL);CHKERRQ(ierr); 188 PetscFunctionReturn(0); 189 } 190 191 #undef __FUNCT__ 192 #define __FUNCT__ "MatDestroy_CHOLMOD" 193 PetscErrorCode MatDestroy_CHOLMOD(Mat F) 194 { 195 PetscErrorCode ierr; 196 Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->spptr; 197 198 PetscFunctionBegin; 199 if (chol) { 200 ierr = !cholmod_X_free_factor(&chol->factor,chol->common);CHKERRQ(ierr); 201 ierr = !cholmod_X_finish(chol->common);CHKERRQ(ierr); 202 ierr = PetscFree(chol->common);CHKERRQ(ierr); 203 ierr = PetscFree(chol->matrix);CHKERRQ(ierr); 204 ierr = (*chol->Destroy)(F);CHKERRQ(ierr); 205 } 206 ierr = PetscFree(F->spptr);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec); 211 212 /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/ 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatFactorInfo_CHOLMOD" 216 static PetscErrorCode MatFactorInfo_CHOLMOD(Mat F,PetscViewer viewer) 217 { 218 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 219 const cholmod_common *c = chol->common; 220 PetscErrorCode ierr; 221 PetscInt i; 222 223 PetscFunctionBegin; 224 if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(0); 225 ierr = PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");CHKERRQ(ierr); 226 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 227 ierr = PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");CHKERRQ(ierr); 228 ierr = PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound);CHKERRQ(ierr); 229 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0);CHKERRQ(ierr); 230 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1);CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2);CHKERRQ(ierr); 232 ierr = PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank);CHKERRQ(ierr); 233 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);CHKERRQ(ierr); 234 ierr = PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal);CHKERRQ(ierr); 235 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis);CHKERRQ(ierr); 236 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super);CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll);CHKERRQ(ierr); 238 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack);CHKERRQ(ierr); 239 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic);CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol);CHKERRQ(ierr); 241 ierr = PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);CHKERRQ(ierr); 243 ierr = PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print);CHKERRQ(ierr); 245 for (i=0; i<c->nmethods; i++) { 246 ierr = PetscViewerASCIIPrintf(viewer,"Ordering method %D%s:\n",i,i==c->selected ? " [SELECTED]" : "");CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", 248 c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2);CHKERRQ(ierr); 249 } 250 ierr = PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder);CHKERRQ(ierr); 251 ierr = PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);CHKERRQ(ierr); 252 /* Statistics */ 253 ierr = PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl);CHKERRQ(ierr); 254 ierr = PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz);CHKERRQ(ierr); 255 ierr = PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl);CHKERRQ(ierr); 257 ierr = PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count);CHKERRQ(ierr);CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage);CHKERRQ(ierr);CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse);CHKERRQ(ierr);CHKERRQ(ierr); 260 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col);CHKERRQ(ierr);CHKERRQ(ierr); 261 ierr = PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);CHKERRQ(ierr);CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);CHKERRQ(ierr);CHKERRQ(ierr); 263 ierr = PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);CHKERRQ(ierr);CHKERRQ(ierr); 264 ierr = PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);CHKERRQ(ierr);CHKERRQ(ierr); 265 #if defined(PETSC_USE_SUITESPARSE_GPU) 266 ierr = PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);CHKERRQ(ierr);CHKERRQ(ierr); 267 #endif 268 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "MatView_CHOLMOD" 274 PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer) 275 { 276 PetscErrorCode ierr; 277 PetscBool iascii; 278 PetscViewerFormat format; 279 280 PetscFunctionBegin; 281 ierr = MatView_SeqSBAIJ(F,viewer);CHKERRQ(ierr); 282 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 283 if (iascii) { 284 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 285 if (format == PETSC_VIEWER_ASCII_INFO) { 286 ierr = MatFactorInfo_CHOLMOD(F,viewer);CHKERRQ(ierr); 287 } 288 } 289 PetscFunctionReturn(0); 290 } 291 292 #undef __FUNCT__ 293 #define __FUNCT__ "MatSolve_CHOLMOD" 294 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 295 { 296 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 297 cholmod_dense cholB,*cholX; 298 PetscScalar *x; 299 PetscErrorCode ierr; 300 301 PetscFunctionBegin; 302 ierr = VecWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 303 static_F = F; 304 cholX = cholmod_X_solve(CHOLMOD_A,chol->factor,&cholB,chol->common); 305 if (!cholX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CHOLMOD failed"); 306 ierr = VecUnWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 307 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 308 ierr = PetscMemcpy(x,cholX->x,cholX->nrow*sizeof(*x));CHKERRQ(ierr); 309 ierr = !cholmod_X_free_dense(&cholX,chol->common);CHKERRQ(ierr); 310 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatCholeskyFactorNumeric_CHOLMOD" 316 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 317 { 318 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 319 cholmod_sparse cholA; 320 PetscBool aijalloc; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc);CHKERRQ(ierr); 325 static_F = F; 326 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 327 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 328 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); 329 330 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 331 332 F->ops->solve = MatSolve_CHOLMOD; 333 F->ops->solvetranspose = MatSolve_CHOLMOD; 334 PetscFunctionReturn(0); 335 } 336 337 #undef __FUNCT__ 338 #define __FUNCT__ "MatCholeskyFactorSymbolic_CHOLMOD" 339 PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 340 { 341 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->spptr; 342 PetscErrorCode ierr; 343 cholmod_sparse cholA; 344 PetscBool aijalloc; 345 PetscInt *fset = 0; 346 size_t fsize = 0; 347 348 PetscFunctionBegin; 349 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);CHKERRQ(ierr); 350 static_F = F; 351 if (chol->factor) { 352 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 353 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 354 } else if (perm) { 355 const PetscInt *ip; 356 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 357 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 358 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 359 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 360 } else { 361 chol->factor = cholmod_X_analyze(&cholA,chol->common); 362 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 363 } 364 365 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 366 367 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNCT__ 372 #define __FUNCT__ "MatFactorGetSolverPackage_seqsbaij_cholmod" 373 PetscErrorCode MatFactorGetSolverPackage_seqsbaij_cholmod(Mat A,const MatSolverPackage *type) 374 { 375 PetscFunctionBegin; 376 *type = MATSOLVERCHOLMOD; 377 PetscFunctionReturn(0); 378 } 379 380 /*MC 381 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 382 via the external package CHOLMOD. 383 384 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 385 386 Use -pc_type lu -pc_factor_mat_solver_package cholmod to use this direct solver 387 388 Consult CHOLMOD documentation for more information about the Common parameters 389 which correspond to the options database keys below. 390 391 Options Database Keys: 392 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 393 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 394 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 395 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 396 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 397 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 398 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 399 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 400 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 401 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 402 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 403 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 404 - -mat_cholmod_print <3> - Verbosity level (None) 405 406 Level: beginner 407 408 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 409 410 .seealso: PCCHOLESKY, PCFactorSetMatSolverPackage(), MatSolverPackage 411 M*/ 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "MatGetFactor_seqsbaij_cholmod" 415 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 416 { 417 Mat B; 418 Mat_CHOLMOD *chol; 419 PetscErrorCode ierr; 420 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 421 422 PetscFunctionBegin; 423 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s", 424 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 425 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 426 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 427 /* Create the factorization matrix F */ 428 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 429 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 430 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 431 ierr = MatSeqSBAIJSetPreallocation(B,1,0,NULL);CHKERRQ(ierr); 432 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 433 434 chol->Wrap = MatWrapCholmod_seqsbaij; 435 chol->Destroy = MatDestroy_SeqSBAIJ; 436 B->spptr = chol; 437 438 B->ops->view = MatView_CHOLMOD; 439 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 440 B->ops->destroy = MatDestroy_CHOLMOD; 441 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqsbaij_cholmod);CHKERRQ(ierr); 442 B->factortype = MAT_FACTOR_CHOLESKY; 443 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 444 B->preallocated = PETSC_TRUE; 445 446 ierr = CholmodStart(B);CHKERRQ(ierr); 447 *F = B; 448 PetscFunctionReturn(0); 449 } 450