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 = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 279 if (iascii) { 280 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 281 if (format == PETSC_VIEWER_ASCII_INFO) { 282 ierr = MatFactorInfo_CHOLMOD(F,viewer);CHKERRQ(ierr); 283 } 284 } 285 PetscFunctionReturn(0); 286 } 287 288 #undef __FUNCT__ 289 #define __FUNCT__ "MatSolve_CHOLMOD" 290 static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X) 291 { 292 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 293 cholmod_dense cholB,*cholX; 294 PetscScalar *x; 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 ierr = VecWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 299 static_F = F; 300 cholX = cholmod_X_solve(CHOLMOD_A,chol->factor,&cholB,chol->common); 301 if (!cholX) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"CHOLMOD failed"); 302 ierr = VecUnWrapCholmodRead(B,&cholB);CHKERRQ(ierr); 303 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 304 ierr = PetscMemcpy(x,cholX->x,cholX->nrow*sizeof(*x));CHKERRQ(ierr); 305 ierr = !cholmod_X_free_dense(&cholX,chol->common);CHKERRQ(ierr); 306 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 307 PetscFunctionReturn(0); 308 } 309 310 #undef __FUNCT__ 311 #define __FUNCT__ "MatCholeskyFactorNumeric_CHOLMOD" 312 static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info) 313 { 314 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 315 cholmod_sparse cholA; 316 PetscBool aijalloc; 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 ierr = (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc);CHKERRQ(ierr); 321 static_F = F; 322 ierr = !cholmod_X_factorize(&cholA,chol->factor,chol->common); 323 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD factorization failed with status %d",chol->common->status); 324 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); 325 326 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 327 328 F->ops->solve = MatSolve_CHOLMOD; 329 F->ops->solvetranspose = MatSolve_CHOLMOD; 330 PetscFunctionReturn(0); 331 } 332 333 #undef __FUNCT__ 334 #define __FUNCT__ "MatCholeskyFactorSymbolic_CHOLMOD" 335 PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info) 336 { 337 Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data; 338 PetscErrorCode ierr; 339 cholmod_sparse cholA; 340 PetscBool aijalloc; 341 PetscInt *fset = 0; 342 size_t fsize = 0; 343 344 PetscFunctionBegin; 345 ierr = (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);CHKERRQ(ierr); 346 static_F = F; 347 if (chol->factor) { 348 ierr = !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common); 349 if (ierr) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 350 } else if (perm) { 351 const PetscInt *ip; 352 ierr = ISGetIndices(perm,&ip);CHKERRQ(ierr); 353 chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common); 354 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 355 ierr = ISRestoreIndices(perm,&ip);CHKERRQ(ierr); 356 } else { 357 chol->factor = cholmod_X_analyze(&cholA,chol->common); 358 if (!chol->factor) SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_LIB,"CHOLMOD analysis failed with status %d",chol->common->status); 359 } 360 361 if (aijalloc) {ierr = PetscFree3(cholA.p,cholA.i,cholA.x);CHKERRQ(ierr);} 362 363 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD; 364 PetscFunctionReturn(0); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "MatFactorGetSolverPackage_seqsbaij_cholmod" 369 static PetscErrorCode MatFactorGetSolverPackage_seqsbaij_cholmod(Mat A,const MatSolverPackage *type) 370 { 371 PetscFunctionBegin; 372 *type = MATSOLVERCHOLMOD; 373 PetscFunctionReturn(0); 374 } 375 376 /*MC 377 MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices 378 via the external package CHOLMOD. 379 380 Use ./configure --download-suitesparse to install PETSc to use CHOLMOD 381 382 Use -pc_type lu -pc_factor_mat_solver_package cholmod to use this direct solver 383 384 Consult CHOLMOD documentation for more information about the Common parameters 385 which correspond to the options database keys below. 386 387 Options Database Keys: 388 + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None) 389 . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None) 390 . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None) 391 . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None) 392 . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None) 393 . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL 394 . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None) 395 . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None) 396 . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None) 397 . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None) 398 . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None) 399 . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None) 400 - -mat_cholmod_print <3> - Verbosity level (None) 401 402 Level: beginner 403 404 Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html 405 406 .seealso: PCCHOLESKY, PCFactorSetMatSolverPackage(), MatSolverPackage 407 M*/ 408 409 #undef __FUNCT__ 410 #define __FUNCT__ "MatGetFactor_seqsbaij_cholmod" 411 PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F) 412 { 413 Mat B; 414 Mat_CHOLMOD *chol; 415 PetscErrorCode ierr; 416 PetscInt m=A->rmap->n,n=A->cmap->n,bs; 417 418 PetscFunctionBegin; 419 if (ftype != MAT_FACTOR_CHOLESKY) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"CHOLMOD cannot do %s factorization with SBAIJ, only %s", 420 MatFactorTypes[ftype],MatFactorTypes[MAT_FACTOR_CHOLESKY]); 421 ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr); 422 if (bs != 1) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"CHOLMOD only supports block size=1, given %D",bs); 423 /* Create the factorization matrix F */ 424 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 425 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 426 ierr = PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);CHKERRQ(ierr); 427 ierr = MatSetUp(B);CHKERRQ(ierr); 428 ierr = PetscNewLog(B,&chol);CHKERRQ(ierr); 429 430 chol->Wrap = MatWrapCholmod_seqsbaij; 431 B->data = chol; 432 433 B->ops->getinfo = MatGetInfo_External; 434 B->ops->view = MatView_CHOLMOD; 435 B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD; 436 B->ops->destroy = MatDestroy_CHOLMOD; 437 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqsbaij_cholmod);CHKERRQ(ierr); 438 B->factortype = MAT_FACTOR_CHOLESKY; 439 B->assembled = PETSC_TRUE; /* required by -ksp_view */ 440 B->preallocated = PETSC_TRUE; 441 442 ierr = CholmodStart(B);CHKERRQ(ierr); 443 444 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 445 ierr = PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);CHKERRQ(ierr); 446 447 *F = B; 448 PetscFunctionReturn(0); 449 } 450