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