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 32-bit pointer systems [as Suitesparse_long is 32-bit 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:   PetscFunctionBegin;
 24:   if (status > CHOLMOD_OK) {
 25:     PetscCallVoid(PetscInfo(static_F, "CHOLMOD warning %d at %s:%d: %s\n", status, file, line, message));
 26:   } else if (status == CHOLMOD_OK) { /* Documentation says this can happen, but why? */
 27:     PetscCallVoid(PetscInfo(static_F, "CHOLMOD OK at %s:%d: %s\n", file, line, message));
 28:   } else {
 29:     PetscCallVoid(PetscErrorPrintf("CHOLMOD error %d at %s:%d: %s\n", status, file, line, message));
 30:   }
 31:   PetscFunctionReturnVoid();
 32: }

 34: #define CHOLMOD_OPTION_DOUBLE(name, help) \
 35:   do { \
 36:     PetscReal tmp = (PetscReal)c->name; \
 37:     PetscCall(PetscOptionsReal("-mat_cholmod_" #name, help, "None", tmp, &tmp, NULL)); \
 38:     c->name = (double)tmp; \
 39:   } while (0)

 41: #define CHOLMOD_OPTION_INT(name, help) \
 42:   do { \
 43:     PetscInt tmp = (PetscInt)c->name; \
 44:     PetscCall(PetscOptionsInt("-mat_cholmod_" #name, help, "None", tmp, &tmp, NULL)); \
 45:     c->name = (int)tmp; \
 46:   } while (0)

 48: #define CHOLMOD_OPTION_SIZE_T(name, help) \
 49:   do { \
 50:     PetscReal tmp = (PetscInt)c->name; \
 51:     PetscCall(PetscOptionsReal("-mat_cholmod_" #name, help, "None", tmp, &tmp, NULL)); \
 52:     PetscCheck(tmp >= 0, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "value must be positive"); \
 53:     c->name = (size_t)tmp; \
 54:   } while (0)

 56: #define CHOLMOD_OPTION_BOOL(name, help) \
 57:   do { \
 58:     PetscBool tmp = (PetscBool) !!c->name; \
 59:     PetscCall(PetscOptionsBool("-mat_cholmod_" #name, help, "None", tmp, &tmp, NULL)); \
 60:     c->name = (int)tmp; \
 61:   } while (0)

 63: static PetscErrorCode CholmodSetOptions(Mat F)
 64: {
 65:   Mat_CHOLMOD    *chol = (Mat_CHOLMOD *)F->data;
 66:   cholmod_common *c    = chol->common;
 67:   PetscBool       flg;

 69:   PetscFunctionBegin;
 70:   PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "CHOLMOD Options", "Mat");
 71:   CHOLMOD_OPTION_INT(nmethods, "Number of different ordering methods to try");

 73: #if defined(PETSC_USE_SUITESPARSE_GPU)
 74:   c->useGPU = 1;
 75:   CHOLMOD_OPTION_INT(useGPU, "Use GPU for BLAS 1, otherwise 0");
 76:   CHOLMOD_OPTION_SIZE_T(maxGpuMemBytes, "Maximum memory to allocate on the GPU");
 77:   CHOLMOD_OPTION_DOUBLE(maxGpuMemFraction, "Fraction of available GPU memory to allocate");
 78: #endif

 80:   /* CHOLMOD handles first-time packing and refactor-packing separately, but we usually want them to be the same. */
 81:   chol->pack = (PetscBool)c->final_pack;
 82:   PetscCall(PetscOptionsBool("-mat_cholmod_pack", "Pack factors after factorization [disable for frequent repeat factorization]", "None", chol->pack, &chol->pack, NULL));
 83:   c->final_pack = (int)chol->pack;

 85:   CHOLMOD_OPTION_DOUBLE(dbound, "Minimum absolute value of diagonal entries of D");
 86:   CHOLMOD_OPTION_DOUBLE(grow0, "Global growth ratio when factors are modified");
 87:   CHOLMOD_OPTION_DOUBLE(grow1, "Column growth ratio when factors are modified");
 88:   CHOLMOD_OPTION_SIZE_T(grow2, "Affine column growth constant when factors are modified");
 89:   CHOLMOD_OPTION_SIZE_T(maxrank, "Max rank of update, larger values are faster but use more memory [2,4,8]");
 90:   {
 91:     static const char *const list[] = {"SIMPLICIAL", "AUTO", "SUPERNODAL", "MatCholmodFactorType", "MAT_CHOLMOD_FACTOR_", 0};
 92:     PetscCall(PetscOptionsEnum("-mat_cholmod_factor", "Factorization method", "None", list, (PetscEnum)c->supernodal, (PetscEnum *)&c->supernodal, NULL));
 93:   }
 94:   if (c->supernodal) CHOLMOD_OPTION_DOUBLE(supernodal_switch, "flop/nnz_L threshold for switching to supernodal factorization");
 95:   CHOLMOD_OPTION_BOOL(final_asis, "Leave factors \"as is\"");
 96:   CHOLMOD_OPTION_BOOL(final_pack, "Pack the columns when finished (use FALSE if the factors will be updated later)");
 97:   if (!c->final_asis) {
 98:     CHOLMOD_OPTION_BOOL(final_super, "Leave supernodal factors instead of converting to simplicial");
 99:     CHOLMOD_OPTION_BOOL(final_ll, "Turn LDL' factorization into LL'");
100:     CHOLMOD_OPTION_BOOL(final_monotonic, "Ensure columns are monotonic when done");
101:     CHOLMOD_OPTION_BOOL(final_resymbol, "Remove numerically zero values resulting from relaxed supernodal amalgamation");
102:   }
103:   {
104:     PetscReal tmp[] = {(PetscReal)c->zrelax[0], (PetscReal)c->zrelax[1], (PetscReal)c->zrelax[2]};
105:     PetscInt  n     = 3;
106:     PetscCall(PetscOptionsRealArray("-mat_cholmod_zrelax", "3 real supernodal relaxed amalgamation parameters", "None", tmp, &n, &flg));
107:     PetscCheck(!flg || n == 3, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "must provide exactly 3 parameters to -mat_cholmod_zrelax");
108:     if (flg)
109:       while (n--) c->zrelax[n] = (double)tmp[n];
110:   }
111:   {
112:     PetscInt n, tmp[] = {(PetscInt)c->nrelax[0], (PetscInt)c->nrelax[1], (PetscInt)c->nrelax[2]};
113:     PetscCall(PetscOptionsIntArray("-mat_cholmod_nrelax", "3 size_t supernodal relaxed amalgamation parameters", "None", tmp, &n, &flg));
114:     PetscCheck(!flg || n == 3, PetscObjectComm((PetscObject)F), PETSC_ERR_ARG_OUTOFRANGE, "must provide exactly 3 parameters to -mat_cholmod_nrelax");
115:     if (flg)
116:       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:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: PetscErrorCode CholmodStart(Mat F)
126: {
127:   Mat_CHOLMOD    *chol = (Mat_CHOLMOD *)F->data;
128:   cholmod_common *c;

130:   PetscFunctionBegin;
131:   if (chol->common) PetscFunctionReturn(PETSC_SUCCESS);
132:   PetscCall(PetscMalloc1(1, &chol->common));
133:   PetscCallExternal(!cholmod_X_start, chol->common);

135:   c                = chol->common;
136:   c->error_handler = CholmodErrorHandler;
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: static PetscErrorCode MatWrapCholmod_seqsbaij(Mat A, PetscBool values, cholmod_sparse *C, PetscBool *aijalloc, PetscBool *valloc)
141: {
142:   Mat_SeqSBAIJ *sbaij    = (Mat_SeqSBAIJ *)A->data;
143:   PetscBool     vallocin = PETSC_FALSE;

145:   PetscFunctionBegin;
146:   PetscCall(PetscMemzero(C, sizeof(*C)));
147:   /* 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 */
148:   C->nrow  = (size_t)A->cmap->n;
149:   C->ncol  = (size_t)A->rmap->n;
150:   C->nzmax = (size_t)sbaij->maxnz;
151:   C->p     = sbaij->i;
152:   C->i     = sbaij->j;
153:   if (values) {
154: #if defined(PETSC_USE_COMPLEX)
155:     /* we need to pass CHOLMOD the conjugate matrix */
156:     PetscScalar *v;
157:     PetscInt     i;

159:     PetscCall(PetscMalloc1(sbaij->maxnz, &v));
160:     for (i = 0; i < sbaij->maxnz; i++) v[i] = PetscConj(sbaij->a[i]);
161:     C->x     = v;
162:     vallocin = PETSC_TRUE;
163: #else
164:     C->x = sbaij->a;
165: #endif
166:   }
167:   C->stype  = -1;
168:   C->itype  = CHOLMOD_INT_TYPE;
169:   C->xtype  = values ? CHOLMOD_SCALAR_TYPE : CHOLMOD_PATTERN;
170:   C->dtype  = CHOLMOD_DOUBLE;
171:   C->sorted = 1;
172:   C->packed = 1;
173:   *aijalloc = PETSC_FALSE;
174:   *valloc   = vallocin;
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: #define GET_ARRAY_READ  0
179: #define GET_ARRAY_WRITE 1

181: PetscErrorCode VecWrapCholmod(Vec X, PetscInt rw, cholmod_dense *Y)
182: {
183:   PetscScalar *x;
184:   PetscInt     n;

186:   PetscFunctionBegin;
187:   PetscCall(PetscMemzero(Y, sizeof(*Y)));
188:   switch (rw) {
189:   case GET_ARRAY_READ:
190:     PetscCall(VecGetArrayRead(X, (const PetscScalar **)&x));
191:     break;
192:   case GET_ARRAY_WRITE:
193:     PetscCall(VecGetArrayWrite(X, &x));
194:     break;
195:   default:
196:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT " not handled", rw);
197:     break;
198:   }
199:   PetscCall(VecGetSize(X, &n));

201:   Y->x     = x;
202:   Y->nrow  = n;
203:   Y->ncol  = 1;
204:   Y->nzmax = n;
205:   Y->d     = n;
206:   Y->xtype = CHOLMOD_SCALAR_TYPE;
207:   Y->dtype = CHOLMOD_DOUBLE;
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: PetscErrorCode VecUnWrapCholmod(Vec X, PetscInt rw, cholmod_dense *Y)
212: {
213:   PetscFunctionBegin;
214:   switch (rw) {
215:   case GET_ARRAY_READ:
216:     PetscCall(VecRestoreArrayRead(X, (const PetscScalar **)&Y->x));
217:     break;
218:   case GET_ARRAY_WRITE:
219:     PetscCall(VecRestoreArrayWrite(X, (PetscScalar **)&Y->x));
220:     break;
221:   default:
222:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT " not handled", rw);
223:     break;
224:   }
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: PetscErrorCode MatDenseWrapCholmod(Mat X, PetscInt rw, cholmod_dense *Y)
229: {
230:   PetscScalar *x;
231:   PetscInt     m, n, lda;

233:   PetscFunctionBegin;
234:   PetscCall(PetscMemzero(Y, sizeof(*Y)));
235:   switch (rw) {
236:   case GET_ARRAY_READ:
237:     PetscCall(MatDenseGetArrayRead(X, (const PetscScalar **)&x));
238:     break;
239:   case GET_ARRAY_WRITE:
240:     PetscCall(MatDenseGetArrayWrite(X, &x));
241:     break;
242:   default:
243:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT " not handled", rw);
244:     break;
245:   }
246:   PetscCall(MatDenseGetLDA(X, &lda));
247:   PetscCall(MatGetLocalSize(X, &m, &n));

249:   Y->x     = x;
250:   Y->nrow  = m;
251:   Y->ncol  = n;
252:   Y->nzmax = lda * n;
253:   Y->d     = lda;
254:   Y->xtype = CHOLMOD_SCALAR_TYPE;
255:   Y->dtype = CHOLMOD_DOUBLE;
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: PetscErrorCode MatDenseUnWrapCholmod(Mat X, PetscInt rw, cholmod_dense *Y)
260: {
261:   PetscFunctionBegin;
262:   switch (rw) {
263:   case GET_ARRAY_READ:
264:     PetscCall(MatDenseRestoreArrayRead(X, (const PetscScalar **)&Y->x));
265:     break;
266:   case GET_ARRAY_WRITE:
267:     /* we don't have MatDenseRestoreArrayWrite */
268:     PetscCall(MatDenseRestoreArray(X, (PetscScalar **)&Y->x));
269:     break;
270:   default:
271:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Case %" PetscInt_FMT " not handled", rw);
272:     break;
273:   }
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: PETSC_INTERN PetscErrorCode MatDestroy_CHOLMOD(Mat F)
278: {
279:   Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;

281:   PetscFunctionBegin;
282:   if (chol->spqrfact) PetscCallExternal(!SuiteSparseQR_C_free, &chol->spqrfact, chol->common);
283:   if (chol->factor) PetscCallExternal(!cholmod_X_free_factor, &chol->factor, chol->common);
284:   if (chol->common->itype == CHOLMOD_INT) {
285:     PetscCallExternal(!cholmod_finish, chol->common);
286:   } else {
287:     PetscCallExternal(!cholmod_l_finish, chol->common);
288:   }
289:   PetscCall(PetscFree(chol->common));
290:   PetscCall(PetscFree(chol->matrix));
291:   PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatFactorGetSolverType_C", NULL));
292:   PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorSymbolic_C", NULL));
293:   PetscCall(PetscObjectComposeFunction((PetscObject)F, "MatQRFactorNumeric_C", NULL));
294:   PetscCall(PetscFree(F->data));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: static PetscErrorCode MatSolve_CHOLMOD(Mat, Vec, Vec);
299: static PetscErrorCode MatMatSolve_CHOLMOD(Mat, Mat, Mat);

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

303: static PetscErrorCode MatView_Info_CHOLMOD(Mat F, PetscViewer viewer)
304: {
305:   Mat_CHOLMOD          *chol = (Mat_CHOLMOD *)F->data;
306:   const cholmod_common *c    = chol->common;
307:   PetscInt              i;

309:   PetscFunctionBegin;
310:   if (F->ops->solve != MatSolve_CHOLMOD) PetscFunctionReturn(PETSC_SUCCESS);
311:   PetscCall(PetscViewerASCIIPrintf(viewer, "CHOLMOD run parameters:\n"));
312:   PetscCall(PetscViewerASCIIPushTab(viewer));
313:   PetscCall(PetscViewerASCIIPrintf(viewer, "Pack factors after symbolic factorization: %s\n", chol->pack ? "TRUE" : "FALSE"));
314:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.dbound            %g  (Smallest absolute value of diagonal entries of D)\n", c->dbound));
315:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.grow0             %g\n", c->grow0));
316:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.grow1             %g\n", c->grow1));
317:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.grow2             %u\n", (unsigned)c->grow2));
318:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.maxrank           %u\n", (unsigned)c->maxrank));
319:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.supernodal_switch %g\n", c->supernodal_switch));
320:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.supernodal        %d\n", c->supernodal));
321:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_asis        %d\n", c->final_asis));
322:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_super       %d\n", c->final_super));
323:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_ll          %d\n", c->final_ll));
324:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_pack        %d\n", c->final_pack));
325:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_monotonic   %d\n", c->final_monotonic));
326:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.final_resymbol    %d\n", c->final_resymbol));
327:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.zrelax            [%g,%g,%g]\n", c->zrelax[0], c->zrelax[1], c->zrelax[2]));
328:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.nrelax            [%u,%u,%u]\n", (unsigned)c->nrelax[0], (unsigned)c->nrelax[1], (unsigned)c->nrelax[2]));
329:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.prefer_upper      %d\n", c->prefer_upper));
330:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.print             %d\n", c->print));
331:   for (i = 0; i < c->nmethods; i++) {
332:     PetscCall(PetscViewerASCIIPrintf(viewer, "Ordering method %" PetscInt_FMT "%s:\n", i, i == c->selected ? " [SELECTED]" : ""));
333:     PetscCall(PetscViewerASCIIPrintf(viewer, "  lnz %g, fl %g, prune_dense %g, prune_dense2 %g\n", c->method[i].lnz, c->method[i].fl, c->method[i].prune_dense, c->method[i].prune_dense2));
334:   }
335:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.postorder         %d\n", c->postorder));
336:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.default_nesdis    %d (use NESDIS instead of METIS for nested dissection)\n", c->default_nesdis));
337:   /* Statistics */
338:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.fl                %g (flop count from most recent analysis)\n", c->fl));
339:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.lnz               %g (fundamental nz in L)\n", c->lnz));
340:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.anz               %g\n", c->anz));
341:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.modfl             %g (flop count from most recent update)\n", c->modfl));
342:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.malloc_count      %g (number of live objects)\n", (double)c->malloc_count));
343:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.memory_usage      %g (peak memory usage in bytes)\n", (double)c->memory_usage));
344:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.memory_inuse      %g (current memory usage in bytes)\n", (double)c->memory_inuse));
345:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.nrealloc_col      %g (number of column reallocations)\n", c->nrealloc_col));
346:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.nrealloc_factor   %g (number of factor reallocations due to column reallocations)\n", c->nrealloc_factor));
347:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.ndbounds_hit      %g (number of times diagonal was modified by dbound)\n", c->ndbounds_hit));
348:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.rowfacfl          %g (number of flops in last call to cholmod_rowfac)\n", c->rowfacfl));
349:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.aatfl             %g (number of flops to compute A(:,f)*A(:,f)')\n", c->aatfl));
350: #if defined(PETSC_USE_SUITESPARSE_GPU)
351:   PetscCall(PetscViewerASCIIPrintf(viewer, "Common.useGPU            %d\n", c->useGPU));
352: #endif
353:   PetscCall(PetscViewerASCIIPopTab(viewer));
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: PETSC_INTERN PetscErrorCode MatView_CHOLMOD(Mat F, PetscViewer viewer)
358: {
359:   PetscBool         iascii;
360:   PetscViewerFormat format;

362:   PetscFunctionBegin;
363:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
364:   if (iascii) {
365:     PetscCall(PetscViewerGetFormat(viewer, &format));
366:     if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_CHOLMOD(F, viewer));
367:   }
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }

371: static PetscErrorCode MatSolve_CHOLMOD(Mat F, Vec B, Vec X)
372: {
373:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
374:   cholmod_dense cholB, cholX, *X_handle, *Y_handle = NULL, *E_handle = NULL;

376:   PetscFunctionBegin;
377:   static_F = F;
378:   PetscCall(VecWrapCholmod(B, GET_ARRAY_READ, &cholB));
379:   PetscCall(VecWrapCholmod(X, GET_ARRAY_WRITE, &cholX));
380:   X_handle = &cholX;
381:   PetscCallExternal(!cholmod_X_solve2, CHOLMOD_A, chol->factor, &cholB, NULL, &X_handle, NULL, &Y_handle, &E_handle, chol->common);
382:   PetscCallExternal(!cholmod_X_free_dense, &Y_handle, chol->common);
383:   PetscCallExternal(!cholmod_X_free_dense, &E_handle, chol->common);
384:   PetscCall(VecUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
385:   PetscCall(VecUnWrapCholmod(X, GET_ARRAY_WRITE, &cholX));
386:   PetscCall(PetscLogFlops(4.0 * chol->common->lnz));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: static PetscErrorCode MatMatSolve_CHOLMOD(Mat F, Mat B, Mat X)
391: {
392:   Mat_CHOLMOD  *chol = (Mat_CHOLMOD *)F->data;
393:   cholmod_dense cholB, cholX, *X_handle, *Y_handle = NULL, *E_handle = NULL;

395:   PetscFunctionBegin;
396:   static_F = F;
397:   PetscCall(MatDenseWrapCholmod(B, GET_ARRAY_READ, &cholB));
398:   PetscCall(MatDenseWrapCholmod(X, GET_ARRAY_WRITE, &cholX));
399:   X_handle = &cholX;
400:   PetscCallExternal(!cholmod_X_solve2, CHOLMOD_A, chol->factor, &cholB, NULL, &X_handle, NULL, &Y_handle, &E_handle, chol->common);
401:   PetscCallExternal(!cholmod_X_free_dense, &Y_handle, chol->common);
402:   PetscCallExternal(!cholmod_X_free_dense, &E_handle, chol->common);
403:   PetscCall(MatDenseUnWrapCholmod(B, GET_ARRAY_READ, &cholB));
404:   PetscCall(MatDenseUnWrapCholmod(X, GET_ARRAY_WRITE, &cholX));
405:   PetscCall(PetscLogFlops(4.0 * B->cmap->n * chol->common->lnz));
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }

409: static PetscErrorCode MatCholeskyFactorNumeric_CHOLMOD(Mat F, Mat A, const MatFactorInfo *info)
410: {
411:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
412:   cholmod_sparse cholA;
413:   PetscBool      aijalloc, valloc;
414:   int            err;

416:   PetscFunctionBegin;
417:   PetscCall((*chol->Wrap)(A, PETSC_TRUE, &cholA, &aijalloc, &valloc));
418:   static_F = F;
419:   err      = !cholmod_X_factorize(&cholA, chol->factor, chol->common);
420:   PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD factorization failed with status %d", chol->common->status);
421:   PetscCheck(chol->common->status != CHOLMOD_NOT_POSDEF, 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);

423:   PetscCall(PetscLogFlops(chol->common->fl));
424:   if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
425:   if (valloc) PetscCall(PetscFree(cholA.x));
426: #if defined(PETSC_USE_SUITESPARSE_GPU)
427:   PetscCall(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));
428: #endif

430:   F->ops->solve             = MatSolve_CHOLMOD;
431:   F->ops->solvetranspose    = MatSolve_CHOLMOD;
432:   F->ops->matsolve          = MatMatSolve_CHOLMOD;
433:   F->ops->matsolvetranspose = MatMatSolve_CHOLMOD;
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: PETSC_INTERN PetscErrorCode MatCholeskyFactorSymbolic_CHOLMOD(Mat F, Mat A, IS perm, const MatFactorInfo *info)
438: {
439:   Mat_CHOLMOD   *chol = (Mat_CHOLMOD *)F->data;
440:   int            err;
441:   cholmod_sparse cholA;
442:   PetscBool      aijalloc, valloc;
443:   PetscInt      *fset  = 0;
444:   size_t         fsize = 0;

446:   PetscFunctionBegin;
447:   /* Set options to F */
448:   PetscCall(CholmodSetOptions(F));

450:   PetscCall((*chol->Wrap)(A, PETSC_FALSE, &cholA, &aijalloc, &valloc));
451:   static_F = F;
452:   if (chol->factor) {
453:     err = !cholmod_X_resymbol(&cholA, fset, fsize, (int)chol->pack, chol->factor, chol->common);
454:     PetscCheck(!err, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed with status %d", chol->common->status);
455:   } else if (perm) {
456:     const PetscInt *ip;
457:     PetscCall(ISGetIndices(perm, &ip));
458:     chol->factor = cholmod_X_analyze_p(&cholA, (PetscInt *)ip, fset, fsize, chol->common);
459:     PetscCheck(chol->factor, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using PETSc ordering with status %d", chol->common->status);
460:     PetscCall(ISRestoreIndices(perm, &ip));
461:   } else {
462:     chol->factor = cholmod_X_analyze(&cholA, chol->common);
463:     PetscCheck(chol->factor, PetscObjectComm((PetscObject)F), PETSC_ERR_LIB, "CHOLMOD analysis failed using internal ordering with status %d", chol->common->status);
464:   }

466:   if (aijalloc) PetscCall(PetscFree2(cholA.p, cholA.i));
467:   if (valloc) PetscCall(PetscFree(cholA.x));

469:   F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_CHOLMOD;
470:   PetscFunctionReturn(PETSC_SUCCESS);
471: }

473: static PetscErrorCode MatFactorGetSolverType_seqsbaij_cholmod(Mat A, MatSolverType *type)
474: {
475:   PetscFunctionBegin;
476:   *type = MATSOLVERCHOLMOD;
477:   PetscFunctionReturn(PETSC_SUCCESS);
478: }

480: PETSC_INTERN PetscErrorCode MatGetInfo_CHOLMOD(Mat F, MatInfoType flag, MatInfo *info)
481: {
482:   Mat_CHOLMOD *chol = (Mat_CHOLMOD *)F->data;

484:   PetscFunctionBegin;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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:
530:    CHOLMOD is part of SuiteSparse http://faculty.cse.tamu.edu/davis/suitesparse.html

532: .seealso: [](ch_matrices), `Mat`, `PCCHOLESKY`, `PCFactorSetMatSolverType()`, `MatSolverType`
533: M*/

535: PETSC_INTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat A, MatFactorType ftype, Mat *F)
536: {
537:   Mat          B;
538:   Mat_CHOLMOD *chol;
539:   PetscInt     m = A->rmap->n, n = A->cmap->n, bs;

541:   PetscFunctionBegin;
542:   PetscCall(MatGetBlockSize(A, &bs));
543:   PetscCheck(bs == 1, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "CHOLMOD only supports block size=1, given %" PetscInt_FMT, bs);
544: #if defined(PETSC_USE_COMPLEX)
545:   PetscCheck(A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for Hermitian matrices");
546: #endif
547:   /* Create the factorization matrix F */
548:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
549:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
550:   PetscCall(PetscStrallocpy("cholmod", &((PetscObject)B)->type_name));
551:   PetscCall(MatSetUp(B));
552:   PetscCall(PetscNew(&chol));

554:   chol->Wrap = MatWrapCholmod_seqsbaij;
555:   B->data    = chol;

557:   B->ops->getinfo                = MatGetInfo_CHOLMOD;
558:   B->ops->view                   = MatView_CHOLMOD;
559:   B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_CHOLMOD;
560:   B->ops->destroy                = MatDestroy_CHOLMOD;
561:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqsbaij_cholmod));
562:   B->factortype   = MAT_FACTOR_CHOLESKY;
563:   B->assembled    = PETSC_TRUE;
564:   B->preallocated = PETSC_TRUE;

566:   PetscCall(CholmodStart(B));

568:   PetscCall(PetscFree(B->solvertype));
569:   PetscCall(PetscStrallocpy(MATSOLVERCHOLMOD, &B->solvertype));
570:   B->canuseordering = PETSC_TRUE;
571:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_CHOLESKY]));
572:   *F = B;
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }