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