Actual source code: sbaijcholmod.c
2: /*
3: Provides an interface to the CHOLMOD sparse solver available through SuiteSparse version 4.2.1
5: When built 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;
21: static void CholmodErrorHandler(int status,const char *file,int line,const char *message)
22: {
23: if (status > CHOLMOD_OK) {
24: PetscInfo(static_F,"CHOLMOD warning %d at %s:%d: %s\n",status,file,line,message);
25: } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */
26: PetscInfo(static_F,"CHOLMOD OK at %s:%d: %s\n",file,line,message);
27: } else {
28: PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n",status,file,line,message);
29: }
30: return;
31: }
33: PetscErrorCode CholmodStart(Mat F)
34: {
36: Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data;
37: cholmod_common *c;
38: PetscBool flg;
40: if (chol->common) return 0;
41: PetscMalloc1(1,&chol->common);
42: !cholmod_X_start(chol->common);
44: c = chol->common;
45: c->error_handler = CholmodErrorHandler;
47: #define CHOLMOD_OPTION_DOUBLE(name,help) do { \
48: PetscReal tmp = (PetscReal)c->name; \
49: PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
50: c->name = (double)tmp; \
51: } while (0)
53: #define CHOLMOD_OPTION_INT(name,help) do { \
54: PetscInt tmp = (PetscInt)c->name; \
55: PetscOptionsInt("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
56: c->name = (int)tmp; \
57: } while (0)
59: #define CHOLMOD_OPTION_SIZE_T(name,help) do { \
60: PetscReal tmp = (PetscInt)c->name; \
61: PetscOptionsReal("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
63: c->name = (size_t)tmp; \
64: } while (0)
66: #define CHOLMOD_OPTION_BOOL(name,help) do { \
67: PetscBool tmp = (PetscBool) !!c->name; \
68: PetscOptionsBool("-mat_cholmod_" #name,help,"None",tmp,&tmp,NULL); \
69: c->name = (int)tmp; \
70: } while (0)
72: PetscOptionsBegin(PetscObjectComm((PetscObject)F),((PetscObject)F)->prefix,"CHOLMOD Options","Mat");
73: CHOLMOD_OPTION_INT(nmethods,"Number of different ordering methods to try");
75: #if defined(PETSC_USE_SUITESPARSE_GPU)
76: c->useGPU = 1;
77: CHOLMOD_OPTION_INT(useGPU,"Use GPU for BLAS 1, otherwise 0");
78: CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes,"Maximum memory to allocate on the GPU");
79: CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction,"Fraction of available GPU memory to allocate");
80: #endif
82: /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */
83: chol->pack = (PetscBool)c->final_pack;
84: PetscOptionsBool("-mat_cholmod_pack","Pack factors after factorization [disable for frequent repeat factorization]","None",chol->pack,&chol->pack,NULL);
85: c->final_pack = (int)chol->pack;
87: CHOLMOD_OPTION_DOUBLE(dbound,"Minimum absolute value of diagonal entries of D");
88: CHOLMOD_OPTION_DOUBLE(grow0,"Global growth ratio when factors are modified");
89: CHOLMOD_OPTION_DOUBLE(grow1,"Column growth ratio when factors are modified");
90: CHOLMOD_OPTION_SIZE_T(grow2,"Affine column growth constant when factors are modified");
91: CHOLMOD_OPTION_SIZE_T(maxrank,"Max rank of update, larger values are faster but use more memory [2,4,8]");
92: {
93: static const char *const list[] = {"SIMPLICIAL","AUTO","SUPERNODAL","MatCholmodFactorType","MAT_CHOLMOD_FACTOR_",0};
94: PetscOptionsEnum("-mat_cholmod_factor","Factorization method","None",list,(PetscEnum)c->supernodal,(PetscEnum*)&c->supernodal,NULL);
95: }
96: if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch,"flop/nnz_L threshold for switching to supernodal factorization");
97: CHOLMOD_OPTION_BOOL(final_asis,"Leave factors \"as is\"");
98: CHOLMOD_OPTION_BOOL(final_pack,"Pack the columns when finished (use FALSE if the factors will be updated later)");
99: if (!c->final_asis) {
100: CHOLMOD_OPTION_BOOL(final_super,"Leave supernodal factors instead of converting to simplicial");
101: CHOLMOD_OPTION_BOOL(final_ll,"Turn LDL' factorization into LL'");
102: CHOLMOD_OPTION_BOOL(final_monotonic,"Ensure columns are monotonic when done");
103: CHOLMOD_OPTION_BOOL(final_resymbol,"Remove numerically zero values resulting from relaxed supernodal amalgamation");
104: }
105: {
106: PetscReal tmp[] = {(PetscReal)c->zrelax[0],(PetscReal)c->zrelax[1],(PetscReal)c->zrelax[2]};
107: PetscInt n = 3;
108: PetscOptionsRealArray("-mat_cholmod_zrelax","3 real supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);
110: if (flg) while (n--) c->zrelax[n] = (double)tmp[n];
111: }
112: {
113: PetscInt n,tmp[] = {(PetscInt)c->nrelax[0],(PetscInt)c->nrelax[1],(PetscInt)c->nrelax[2]};
114: PetscOptionsIntArray("-mat_cholmod_nrelax","3 size_t supernodal relaxed amalgamation parameters","None",tmp,&n,&flg);
116: if (flg) while (n--) c->nrelax[n] = (size_t)tmp[n];
117: }
118: CHOLMOD_OPTION_BOOL(prefer_upper,"Work with upper triangular form [faster when using fill-reducing ordering, slower in natural ordering]");
119: CHOLMOD_OPTION_BOOL(default_nesdis,"Use NESDIS instead of METIS for nested dissection");
120: CHOLMOD_OPTION_INT(print,"Verbosity level");
121: PetscOptionsEnd();
122: return 0;
123: }
125: static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A,PetscBool values,cholmod_sparse *C,PetscBool *aijalloc,PetscBool *valloc)
126: {
127: Mat_SeqSBAIJ *sbaij = (Mat_SeqSBAIJ*)A->data;
128: PetscBool vallocin = PETSC_FALSE;
130: PetscMemzero(C,sizeof(*C));
131: /* 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 */
132: C->nrow = (size_t)A->cmap->n;
133: C->ncol = (size_t)A->rmap->n;
134: C->nzmax = (size_t)sbaij->maxnz;
135: C->p = sbaij->i;
136: C->i = sbaij->j;
137: if (values) {
138: #if defined(PETSC_USE_COMPLEX)
139: /* we need to pass CHOLMOD the conjugate matrix */
140: PetscScalar *v;
141: PetscInt i;
143: PetscMalloc1(sbaij->maxnz,&v);
144: for (i = 0; i < sbaij->maxnz; i++) v[i] = PetscConj(sbaij->a[i]);
145: C->x = v;
146: vallocin = PETSC_TRUE;
147: #else
148: C->x = sbaij->a;
149: #endif
150: }
151: C->stype = -1;
152: C->itype = CHOLMOD_INT_TYPE;
153: C->xtype = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
154: C->dtype = CHOLMOD_DOUBLE;
155: C->sorted = 1;
156: C->packed = 1;
157: *aijalloc = PETSC_FALSE;
158: *valloc = vallocin;
159: return 0;
160: }
162: #define GET_ARRAY_READ 0
163: #define GET_ARRAY_WRITE 1
165: PetscErrorCode VecWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y)
166: {
167: PetscScalar *x;
168: PetscInt n;
170: PetscMemzero(Y,sizeof(*Y));
171: switch (rw) {
172: case GET_ARRAY_READ:
173: VecGetArrayRead(X,(const PetscScalar**)&x);
174: break;
175: case GET_ARRAY_WRITE:
176: VecGetArrayWrite(X,&x);
177: break;
178: default:
179: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw);
180: break;
181: }
182: VecGetSize(X,&n);
184: Y->x = x;
185: Y->nrow = n;
186: Y->ncol = 1;
187: Y->nzmax = n;
188: Y->d = n;
189: Y->xtype = CHOLMOD_SCALAR_TYPE;
190: Y->dtype = CHOLMOD_DOUBLE;
191: return 0;
192: }
194: PetscErrorCode VecUnWrapCholmod(Vec X,PetscInt rw,cholmod_dense *Y)
195: {
196: switch (rw) {
197: case GET_ARRAY_READ:
198: VecRestoreArrayRead(X,(const PetscScalar**)&Y->x);
199: break;
200: case GET_ARRAY_WRITE:
201: VecRestoreArrayWrite(X,(PetscScalar**)&Y->x);
202: break;
203: default:
204: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw);
205: break;
206: }
207: return 0;
208: }
210: PetscErrorCode MatDenseWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y)
211: {
212: PetscScalar *x;
213: PetscInt m,n,lda;
215: PetscMemzero(Y,sizeof(*Y));
216: switch (rw) {
217: case GET_ARRAY_READ:
218: MatDenseGetArrayRead(X,(const PetscScalar**)&x);
219: break;
220: case GET_ARRAY_WRITE:
221: MatDenseGetArrayWrite(X,&x);
222: break;
223: default:
224: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw);
225: break;
226: }
227: MatDenseGetLDA(X,&lda);
228: MatGetLocalSize(X,&m,&n);
230: Y->x = x;
231: Y->nrow = m;
232: Y->ncol = n;
233: Y->nzmax = lda*n;
234: Y->d = lda;
235: Y->xtype = CHOLMOD_SCALAR_TYPE;
236: Y->dtype = CHOLMOD_DOUBLE;
237: return 0;
238: }
240: PetscErrorCode MatDenseUnWrapCholmod(Mat X,PetscInt rw,cholmod_dense *Y)
241: {
242: switch (rw) {
243: case GET_ARRAY_READ:
244: MatDenseRestoreArrayRead(X,(const PetscScalar**)&Y->x);
245: break;
246: case GET_ARRAY_WRITE:
247: /* we don't have MatDenseRestoreArrayWrite */
248: MatDenseRestoreArray(X,(PetscScalar**)&Y->x);
249: break;
250: default:
251: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Case %" PetscInt_FMT " not handled",rw);
252: break;
253: }
254: return 0;
255: }
257: PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F)
258: {
259: Mat_CHOLMOD *chol=(Mat_CHOLMOD*)F->data;
261: if (chol->spqrfact) {
262: !SuiteSparseQR_C_free(&chol->spqrfact, chol->common);
263: }
264: if (chol->factor) {
265: !cholmod_X_free_factor(&chol->factor,chol->common);
266: }
267: if (chol->common->itype == CHOLMOD_INT) {
268: !cholmod_finish(chol->common);
269: } else {
270: !cholmod_l_finish(chol->common);
271: }
272: PetscFree(chol->common);
273: PetscFree(chol->matrix);
274: PetscObjectComposeFunction((PetscObject)F,"MatFactorGetSolverType_C",NULL);
275: PetscFree(F->data);
276: return 0;
277: }
279: static PetscErrorCode MatSolve_CHOLMOD(Mat,Vec,Vec);
280: static PetscErrorCode MatMatSolve_CHOLMOD(Mat,Mat,Mat);
282: /*static const char *const CholmodOrderingMethods[] = {"User","AMD","METIS","NESDIS(default)","Natural","NESDIS(small=20000)","NESDIS(small=4,no constrained)","NESDIS()"};*/
284: static PetscErrorCode MatView_Info_CHOLMOD(Mat F,PetscViewer viewer)
285: {
286: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
287: const cholmod_common *c = chol->common;
288: PetscErrorCode ierr;
289: PetscInt i;
291: if (F->ops->solve != MatSolve_CHOLMOD) return 0;
292: PetscViewerASCIIPrintf(viewer,"CHOLMOD run parameters:\n");
293: PetscViewerASCIIPushTab(viewer);
294: PetscViewerASCIIPrintf(viewer,"Pack factors after symbolic factorization: %s\n",chol->pack ? "TRUE" : "FALSE");
295: PetscViewerASCIIPrintf(viewer,"Common.dbound %g (Smallest absolute value of diagonal entries of D)\n",c->dbound);
296: PetscViewerASCIIPrintf(viewer,"Common.grow0 %g\n",c->grow0);
297: PetscViewerASCIIPrintf(viewer,"Common.grow1 %g\n",c->grow1);
298: PetscViewerASCIIPrintf(viewer,"Common.grow2 %u\n",(unsigned)c->grow2);
299: PetscViewerASCIIPrintf(viewer,"Common.maxrank %u\n",(unsigned)c->maxrank);
300: PetscViewerASCIIPrintf(viewer,"Common.supernodal_switch %g\n",c->supernodal_switch);
301: PetscViewerASCIIPrintf(viewer,"Common.supernodal %d\n",c->supernodal);
302: PetscViewerASCIIPrintf(viewer,"Common.final_asis %d\n",c->final_asis);
303: PetscViewerASCIIPrintf(viewer,"Common.final_super %d\n",c->final_super);
304: PetscViewerASCIIPrintf(viewer,"Common.final_ll %d\n",c->final_ll);
305: PetscViewerASCIIPrintf(viewer,"Common.final_pack %d\n",c->final_pack);
306: PetscViewerASCIIPrintf(viewer,"Common.final_monotonic %d\n",c->final_monotonic);
307: PetscViewerASCIIPrintf(viewer,"Common.final_resymbol %d\n",c->final_resymbol);
308: PetscViewerASCIIPrintf(viewer,"Common.zrelax [%g,%g,%g]\n",c->zrelax[0],c->zrelax[1],c->zrelax[2]);
309: PetscViewerASCIIPrintf(viewer,"Common.nrelax [%u,%u,%u]\n",(unsigned)c->nrelax[0],(unsigned)c->nrelax[1],(unsigned)c->nrelax[2]);
310: PetscViewerASCIIPrintf(viewer,"Common.prefer_upper %d\n",c->prefer_upper);
311: PetscViewerASCIIPrintf(viewer,"Common.print %d\n",c->print);
312: for (i=0; i<c->nmethods; i++) {
313: PetscViewerASCIIPrintf(viewer,"Ordering method %" PetscInt_FMT "%s:\n",i,i==c->selected ? " [SELECTED]" : "");
314: PetscViewerASCIIPrintf(viewer," lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n",
315: c->method[i].lnz,c->method[i].fl,c->method[i].prune_dense,c->method[i].prune_dense2);
316: }
317: PetscViewerASCIIPrintf(viewer,"Common.postorder %d\n",c->postorder);
318: PetscViewerASCIIPrintf(viewer,"Common.default_nesdis %d (use NESDIS instead of METIS for nested dissection)\n",c->default_nesdis);
319: /* Statistics */
320: PetscViewerASCIIPrintf(viewer,"Common.fl %g (flop count from most recent analysis)\n",c->fl);
321: PetscViewerASCIIPrintf(viewer,"Common.lnz %g (fundamental nz in L)\n",c->lnz);
322: PetscViewerASCIIPrintf(viewer,"Common.anz %g\n",c->anz);
323: PetscViewerASCIIPrintf(viewer,"Common.modfl %g (flop count from most recent update)\n",c->modfl);
324: PetscViewerASCIIPrintf(viewer,"Common.malloc_count %g (number of live objects)\n",(double)c->malloc_count);
325: PetscViewerASCIIPrintf(viewer,"Common.memory_usage %g (peak memory usage in bytes)\n",(double)c->memory_usage);
326: PetscViewerASCIIPrintf(viewer,"Common.memory_inuse %g (current memory usage in bytes)\n",(double)c->memory_inuse);
327: PetscViewerASCIIPrintf(viewer,"Common.nrealloc_col %g (number of column reallocations)\n",c->nrealloc_col);
328: PetscViewerASCIIPrintf(viewer,"Common.nrealloc_factor %g (number of factor reallocations due to column reallocations)\n",c->nrealloc_factor);
329: PetscViewerASCIIPrintf(viewer,"Common.ndbounds_hit %g (number of times diagonal was modified by dbound)\n",c->ndbounds_hit);
330: PetscViewerASCIIPrintf(viewer,"Common.rowfacfl %g (number of flops in last call to cholmod_rowfac)\n",c->rowfacfl);
331: PetscViewerASCIIPrintf(viewer,"Common.aatfl %g (number of flops to compute A(:,f)*A(:,f)')\n",c->aatfl);
332: #if defined(PETSC_USE_SUITESPARSE_GPU)
333: PetscViewerASCIIPrintf(viewer,"Common.useGPU %d\n",c->useGPU);
334: #endif
335: PetscViewerASCIIPopTab(viewer);
336: return 0;
337: }
339: PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F,PetscViewer viewer)
340: {
341: PetscBool iascii;
342: PetscViewerFormat format;
344: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
345: if (iascii) {
346: PetscViewerGetFormat(viewer,&format);
347: if (format == PETSC_VIEWER_ASCII_INFO) {
348: MatView_Info_CHOLMOD(F,viewer);
349: }
350: }
351: return 0;
352: }
354: static PetscErrorCode MatSolve_CHOLMOD(Mat F,Vec B,Vec X)
355: {
356: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
357: cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL;
359: static_F = F;
360: VecWrapCholmod(B,GET_ARRAY_READ,&cholB);
361: VecWrapCholmod(X,GET_ARRAY_WRITE,&cholX);
362: X_handle = &cholX;
363: !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);
364: !cholmod_X_free_dense(&Y_handle,chol->common);
365: !cholmod_X_free_dense(&E_handle,chol->common);
366: VecUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
367: VecUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);
368: PetscLogFlops(4.0*chol->common->lnz);
369: return 0;
370: }
372: static PetscErrorCode MatMatSolve_CHOLMOD(Mat F,Mat B,Mat X)
373: {
374: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
375: cholmod_dense cholB,cholX,*X_handle,*Y_handle = NULL,*E_handle = NULL;
377: static_F = F;
378: MatDenseWrapCholmod(B,GET_ARRAY_READ,&cholB);
379: MatDenseWrapCholmod(X,GET_ARRAY_WRITE,&cholX);
380: X_handle = &cholX;
381: !cholmod_X_solve2(CHOLMOD_A,chol->factor,&cholB,NULL,&X_handle,NULL,&Y_handle,&E_handle,chol->common);
382: !cholmod_X_free_dense(&Y_handle,chol->common);
383: !cholmod_X_free_dense(&E_handle,chol->common);
384: MatDenseUnWrapCholmod(B,GET_ARRAY_READ,&cholB);
385: MatDenseUnWrapCholmod(X,GET_ARRAY_WRITE,&cholX);
386: PetscLogFlops(4.0*B->cmap->n*chol->common->lnz);
387: return 0;
388: }
390: static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F,Mat A,const MatFactorInfo *info)
391: {
392: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
393: cholmod_sparse cholA;
394: PetscBool aijalloc,valloc;
397: (*chol->Wrap)(A,PETSC_TRUE,&cholA,&aijalloc,&valloc);
398: static_F = F;
399: !cholmod_X_factorize(&cholA,chol->factor,chol->common);
403: PetscLogFlops(chol->common->fl);
404: if (aijalloc) PetscFree2(cholA.p,cholA.i);
405: if (valloc) PetscFree(cholA.x);
406: #if defined(PETSC_USE_SUITESPARSE_GPU)
407: PetscLogGpuTimeAdd(chol->common->CHOLMOD_GPU_GEMM_TIME + chol->common->CHOLMOD_GPU_SYRK_TIME + chol->common->CHOLMOD_GPU_TRSM_TIME + chol->common->CHOLMOD_GPU_POTRF_TIME);
408: #endif
410: F->ops->solve = MatSolve_CHOLMOD;
411: F->ops->solvetranspose = MatSolve_CHOLMOD;
412: F->ops->matsolve = MatMatSolve_CHOLMOD;
413: F->ops->matsolvetranspose = MatMatSolve_CHOLMOD;
414: return 0;
415: }
417: PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F,Mat A,IS perm,const MatFactorInfo *info)
418: {
419: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
421: cholmod_sparse cholA;
422: PetscBool aijalloc,valloc;
423: PetscInt *fset = 0;
424: size_t fsize = 0;
426: (*chol->Wrap)(A,PETSC_FALSE,&cholA,&aijalloc,&valloc);
427: static_F = F;
428: if (chol->factor) {
429: !cholmod_X_resymbol(&cholA,fset,fsize,(int)chol->pack,chol->factor,chol->common);
431: } else if (perm) {
432: const PetscInt *ip;
433: ISGetIndices(perm,&ip);
434: chol->factor = cholmod_X_analyze_p(&cholA,(PetscInt*)ip,fset,fsize,chol->common);
436: ISRestoreIndices(perm,&ip);
437: } else {
438: chol->factor = cholmod_X_analyze(&cholA,chol->common);
440: }
442: if (aijalloc) PetscFree2(cholA.p,cholA.i);
443: if (valloc) PetscFree(cholA.x);
445: F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD;
446: return 0;
447: }
449: static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A,MatSolverType *type)
450: {
451: *type = MATSOLVERCHOLMOD;
452: return 0;
453: }
455: PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F,MatInfoType flag,MatInfo *info)
456: {
457: Mat_CHOLMOD *chol = (Mat_CHOLMOD*)F->data;
459: info->block_size = 1.0;
460: info->nz_allocated = chol->common->lnz;
461: info->nz_used = chol->common->lnz;
462: info->nz_unneeded = 0.0;
463: info->assemblies = 0.0;
464: info->mallocs = 0.0;
465: info->memory = chol->common->memory_inuse;
466: info->fill_ratio_given = 0;
467: info->fill_ratio_needed = 0;
468: info->factor_mallocs = chol->common->malloc_count;
469: return 0;
470: }
472: /*MC
473: MATSOLVERCHOLMOD
475: A matrix type providing direct solvers (Cholesky) for sequential matrices
476: via the external package CHOLMOD.
478: Use ./configure --download-suitesparse to install PETSc to use CHOLMOD
480: Use -pc_type cholesky -pc_factor_mat_solver_type cholmod to use this direct solver
482: Consult CHOLMOD documentation for more information about the Common parameters
483: which correspond to the options database keys below.
485: Options Database Keys:
486: + -mat_cholmod_dbound <0> - Minimum absolute value of diagonal entries of D (None)
487: . -mat_cholmod_grow0 <1.2> - Global growth ratio when factors are modified (None)
488: . -mat_cholmod_grow1 <1.2> - Column growth ratio when factors are modified (None)
489: . -mat_cholmod_grow2 <5> - Affine column growth constant when factors are modified (None)
490: . -mat_cholmod_maxrank <8> - Max rank of update, larger values are faster but use more memory [2,4,8] (None)
491: . -mat_cholmod_factor <AUTO> - (choose one of) SIMPLICIAL AUTO SUPERNODAL
492: . -mat_cholmod_supernodal_switch <40> - flop/nnz_L threshold for switching to supernodal factorization (None)
493: . -mat_cholmod_final_asis <TRUE> - Leave factors "as is" (None)
494: . -mat_cholmod_final_pack <TRUE> - Pack the columns when finished (use FALSE if the factors will be updated later) (None)
495: . -mat_cholmod_zrelax <0.8> - 3 real supernodal relaxed amalgamation parameters (None)
496: . -mat_cholmod_nrelax <4> - 3 size_t supernodal relaxed amalgamation parameters (None)
497: . -mat_cholmod_prefer_upper <TRUE> - Work with upper triangular form (faster when using fill-reducing ordering, slower in natural ordering) (None)
498: . -mat_cholmod_print <3> - Verbosity level (None)
499: - -mat_ordering_type internal - Use the ordering provided by Cholmod
501: Level: beginner
503: Note: CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html
505: .seealso: PCCHOLESKY, PCFactorSetMatSolverType(), MatSolverType
506: M*/
508: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A,MatFactorType ftype,Mat *F)
509: {
510: Mat B;
511: Mat_CHOLMOD *chol;
512: PetscInt m=A->rmap->n,n=A->cmap->n,bs;
513: const char *prefix;
515: MatGetBlockSize(A,&bs);
517: #if defined(PETSC_USE_COMPLEX)
519: #endif
520: /* Create the factorization matrix F */
521: MatCreate(PetscObjectComm((PetscObject)A),&B);
522: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
523: PetscStrallocpy("cholmod",&((PetscObject)B)->type_name);
524: MatGetOptionsPrefix(A,&prefix);
525: MatSetOptionsPrefix(B,prefix);
526: MatSetUp(B);
527: PetscNewLog(B,&chol);
529: chol->Wrap = MatWrapCholmod_seqsbaij;
530: B->data = chol;
532: B->ops->getinfo = MatGetInfo_CHOLMOD;
533: B->ops->view = MatView_CHOLMOD;
534: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
535: B->ops->destroy = MatDestroy_CHOLMOD;
536: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_seqsbaij_cholmod);
537: B->factortype = MAT_FACTOR_CHOLESKY;
538: B->assembled = PETSC_TRUE;
539: B->preallocated = PETSC_TRUE;
541: CholmodStart(B);
543: PetscFree(B->solvertype);
544: PetscStrallocpy(MATSOLVERCHOLMOD,&B->solvertype);
545: B->canuseordering = PETSC_TRUE;
546: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_CHOLESKY]);
547: *F = B;
548: return 0;
549: }