Actual source code: sbaijcholmod.c
petsc-3.6.1 2015-08-06
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: 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: 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: 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: 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_EXTERN 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);
447: *F = B;
448: return(0);
449: }