Actual source code: sbaijcholmod.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  2: /*
  3:    Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1

  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]

 10: */

 12: #include <../src/mat/impls/sbaij/seq/sbaij.h>
 13: #include <../src/mat/impls/sbaij/seq/cholmod/cholmodimpl.h>

 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;

 23: static void CholmodErrorHandler(int status,const char *file,int line,const char *message)
 24: {

 28:   if (status > CHOLMOD_OK) {
 29:     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:     PetscInfo3(static_F,"CHOLMOD OK at %s:%d: %s\n",file,line,message);CHKERRV(ierr);
 32:   } else {
 33:     PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message);CHKERRV(ierr);
 34:   }
 35:   PetscFunctionReturnVoid();
 36: }

 40: PetscErrorCode  CholmodStart(Mat F)
 41: {
 43:   Mat_CHOLMOD    *chol=(Mat_CHOLMOD*)F->spptr;
 44:   cholmod_common *c;
 45:   PetscBool      flg;

 48:   if (chol->common) return(0);
 49:   PetscMalloc1(1,&chol->common);
 50:   !cholmod_X_start(chol->common);

 52:   c                = chol->common;
 53:   c->error_handler = CholmodErrorHandler;

 55: #define CHOLMOD_OPTION_DOUBLE(name,help) do {                            \
 56:     PetscReal tmp = (PetscReal)c->name;                                  \
 57:     PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
 58:     c->name = (double)tmp;                                               \
 59: } while (0)

 61: #define CHOLMOD_OPTION_INT(name,help) do {                               \
 62:     PetscInt tmp = (PetscInt)c->name;                                    \
 63:     PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
 64:     c->name = (int)tmp;                                                  \
 65: } while (0)

 67: #define CHOLMOD_OPTION_SIZE_T(name,help) do {                            \
 68:     PetscInt tmp = (PetscInt)c->name;                                    \
 69:     PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
 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)

 74: #define CHOLMOD_OPTION_BOOL(name,help) do {                             \
 75:     PetscBool tmp = (PetscBool) !!c->name;                              \
 76:     PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
 77:     c->name = (int)tmp;                                                  \
 78: } while (0)

 80:   PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat");
 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;

 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

 89:   PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL);
 90:   c->final_pack = (int)chol->pack;

 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:     PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL);
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:     PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);
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:     PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);
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:   PetscOptionsEnd();
127:   return(0);
128: }

132: static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool  *aijalloc)
133: {
134:   Mat_SeqSBAIJ   *sbaij = (Mat_SeqSBAIJ*)A->data;

138:   PetscMemzero(C,sizeof(*C));
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:   return(0);
154: }

158: static PetscErrorCode VecWrapCholmodRead(Vec X,cholmod_dense *Y)
159: {
160:   PetscErrorCode    ierr;
161:   const PetscScalar *x;
162:   PetscInt          n;

165:   PetscMemzero(Y,sizeof(*Y));
166:   VecGetArrayRead(X,&x);
167:   VecGetSize(X,&n);

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:   return(0);
178: }

182: static PetscErrorCode VecUnWrapCholmodRead(Vec X,cholmod_dense *Y)
183: {
184:   PetscErrorCode    ierr;

187:   VecRestoreArrayRead(X,PETSC_NULL);
188:   return(0);
189: }

193: PETSC_INTERN PetscErrorCode  MatDestroy_CHOLMOD(Mat F)
194: {
196:   Mat_CHOLMOD    *chol=(Mat_CHOLMOD*)F->spptr;

199:   if (chol) {
200:     !cholmod_X_free_factor(&chol->factor,chol->common);
201:     !cholmod_X_finish(chol->common);
202:     PetscFree(chol->common);
203:     PetscFree(chol->matrix);
204:     (*chol->Destroy)(F);
205:   }
206:   PetscFree(F->spptr);
207:   return(0);
208: }

210: static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec);

212: /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/

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;

224:   if (F->ops->solve != MatSolve_CHOLMOD) return(0);
225:   PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");
226:   PetscViewerASCIIPushTab(viewer);
227:   PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");
228:   PetscViewerASCIIPrintf(viewer,"Common.dbound            %g  (Smallest absolute value of diagonal entries of D)\n",c->dbound);
229:   PetscViewerASCIIPrintf(viewer,"Common.grow0             %g\n",c->grow0);
230:   PetscViewerASCIIPrintf(viewer,"Common.grow1             %g\n",c->grow1);
231:   PetscViewerASCIIPrintf(viewer,"Common.grow2             %u\n",(unsigned)c->grow2);
232:   PetscViewerASCIIPrintf(viewer,"Common.maxrank           %u\n",(unsigned)c->maxrank);
233:   PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);
234:   PetscViewerASCIIPrintf(viewer,"Common.supernodal        %d\n",c->supernodal);
235:   PetscViewerASCIIPrintf(viewer,"Common.final_asis        %d\n",c->final_asis);
236:   PetscViewerASCIIPrintf(viewer,"Common.final_super       %d\n",c->final_super);
237:   PetscViewerASCIIPrintf(viewer,"Common.final_ll          %d\n",c->final_ll);
238:   PetscViewerASCIIPrintf(viewer,"Common.final_pack        %d\n",c->final_pack);
239:   PetscViewerASCIIPrintf(viewer,"Common.final_monotonic   %d\n",c->final_monotonic);
240:   PetscViewerASCIIPrintf(viewer,"Common.final_resymbol    %d\n",c->final_resymbol);
241:   PetscViewerASCIIPrintf(viewer,"Common.zrelax            [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);
242:   PetscViewerASCIIPrintf(viewer,"Common.nrelax            [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);
243:   PetscViewerASCIIPrintf(viewer,"Common.prefer_upper      %d\n",c->prefer_upper);
244:   PetscViewerASCIIPrintf(viewer,"Common.print             %d\n",c->print);
245:   for (i=0; i<c->nmethods; i++) {
246:     PetscViewerASCIIPrintf(viewer,"Ordering method %D%s:\n",i,i==c->selected ? " [SELECTED]" : "");
247:     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);
249:   }
250:   PetscViewerASCIIPrintf(viewer,"Common.postorder         %d\n",c->postorder);
251:   PetscViewerASCIIPrintf(viewer,"Common.default_nesdis    %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);
252:   /* Statistics */
253:   PetscViewerASCIIPrintf(viewer,"Common.fl                %g (flop count from most recent analysis)\n",c->fl);
254:   PetscViewerASCIIPrintf(viewer,"Common.lnz               %g (fundamental nz in L)\n",c->lnz);
255:   PetscViewerASCIIPrintf(viewer,"Common.anz               %g\n",c->anz);
256:   PetscViewerASCIIPrintf(viewer,"Common.modfl             %g (flop count from most recent update)\n",c->modfl);
257:   PetscViewerASCIIPrintf(viewer,"Common.malloc_count      %g (number of live objects)\n",(double)c->malloc_count);
258:   PetscViewerASCIIPrintf(viewer,"Common.memory_usage      %g (peak memory usage in bytes)\n",(double)c->memory_usage);
259:   PetscViewerASCIIPrintf(viewer,"Common.memory_inuse      %g (current memory usage in bytes)\n",(double)c->memory_inuse);
260:   PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col      %g (number of column reallocations)\n",c->nrealloc_col);
261:   PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor   %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);
262:   PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit      %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);
263:   PetscViewerASCIIPrintf(viewer,"Common.rowfacfl          %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);
264:   PetscViewerASCIIPrintf(viewer,"Common.aatfl             %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);
265: #if defined(PETSC_USE_SUITESPARSE_GPU)
266:   PetscViewerASCIIPrintf(viewer,"Common.useGPU            %d\n",c->useGPU);
267: #endif
268:   PetscViewerASCIIPopTab(viewer);
269:   return(0);
270: }

274: PETSC_INTERN PetscErrorCode  MatView_CHOLMOD(Mat F,PetscViewer viewer)
275: {
276:   PetscErrorCode    ierr;
277:   PetscBool         iascii;
278:   PetscViewerFormat format;

281:   MatView_SeqSBAIJ(F,viewer);
282:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
283:   if (iascii) {
284:     PetscViewerGetFormat(viewer,&format);
285:     if (format == PETSC_VIEWER_ASCII_INFO) {
286:       MatFactorInfo_CHOLMOD(F,viewer);
287:     }
288:   }
289:   return(0);
290: }

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;

302:   VecWrapCholmodRead(B,&cholB);
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:   VecUnWrapCholmodRead(B,&cholB);
307:   VecGetArray(X,&x);
308:   PetscMemcpy(x,cholX->x,cholX->nrow*sizeof(*x));
309:   !cholmod_X_free_dense(&cholX,chol->common);
310:   VecRestoreArray(X,&x);
311:   return(0);
312: }

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;

324:   (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc);
325:   static_F = F;
326:   !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);

330:   if (aijalloc) {PetscFree3(cholA.p,cholA.i,cholA.x);}

332:   F->ops->solve          = MatSolve_CHOLMOD;
333:   F->ops->solvetranspose = MatSolve_CHOLMOD;
334:   return(0);
335: }

339: PETSC_INTERN PetscErrorCode  MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info)
340: {
341:   Mat_CHOLMOD    *chol = (Mat_CHOLMOD*)F->spptr;
343:   cholmod_sparse cholA;
344:   PetscBool      aijalloc;
345:   PetscInt       *fset = 0;
346:   size_t         fsize = 0;

349:   (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc);
350:   static_F = F;
351:   if (chol->factor) {
352:     !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:     ISGetIndices(perm,&ip);
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:     ISRestoreIndices(perm,&ip);
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:   }

365:   if (aijalloc) {PetscFree3(cholA.p,cholA.i,cholA.x);}

367:   F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD;
368:   return(0);
369: }

373: static PetscErrorCode MatFactorGetSolverPackage_seqsbaij_cholmod(Mat A,const MatSolverPackage *type)
374: {
376:   *type = MATSOLVERCHOLMOD;
377:   return(0);
378: }

380: /*MC
381:   MATSOLVERCHOLMOD = "cholmod" - A matrix type providing direct solvers (Cholesky) for sequential matrices
382:   via the external package CHOLMOD.

384:   Use ./configure --download-suitesparse to install PETSc to use CHOLMOD

386:   Use -pc_type lu -pc_factor_mat_solver_package cholmod to use this direct solver

388:   Consult CHOLMOD documentation for more information about the Common parameters
389:   which correspond to the options database keys below.

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)

406:    Level: beginner

408:    Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

410: .seealso: PCCHOLESKY, PCFactorSetMatSolverPackage(), MatSolverPackage
411: M*/

415: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F)
416: {
417:   Mat            B;
418:   Mat_CHOLMOD    *chol;
420:   PetscInt       m=A->rmap->n,n=A->cmap->n,bs;

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:   MatGetBlockSize(A,&bs);
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:   MatCreate(PetscObjectComm((PetscObject)A),&B);
429:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
430:   MatSetType(B,((PetscObject)A)->type_name);
431:   MatSeqSBAIJSetPreallocation(B,1,0,NULL);
432:   PetscNewLog(B,&chol);

434:   chol->Wrap    = MatWrapCholmod_seqsbaij;
435:   chol->Destroy = MatDestroy_SeqSBAIJ;
436:   B->spptr      = chol;

438:   B->ops->view                   = MatView_CHOLMOD;
439:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
440:   B->ops->destroy                = MatDestroy_CHOLMOD;
441:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_seqsbaij_cholmod);
442:   B->factortype                  = MAT_FACTOR_CHOLESKY;
443:   B->assembled                   = PETSC_TRUE; /* required by -ksp_view */
444:   B->preallocated                = PETSC_TRUE;

446:   CholmodStart(B);

448:   PetscFree(B->solvertype);
449:   PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);

451:   *F   = B;
452:   return(0);
453: }