Actual source code: superlu_dist.c
1: /*
2: Provides an interface to the SuperLU_DIST sparse solver
3: */
5: #include <../src/mat/impls/aij/seq/aij.h>
6: #include <../src/mat/impls/aij/mpi/mpiaij.h>
7: #include <petscpkg_version.h>
9: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Wundef")
10: EXTERN_C_BEGIN
11: // To suppress compiler warnings
12: #define DISABLE_CUSPARSE_DEPRECATED
13: #if PetscDefined(USE_COMPLEX)
14: #define CASTDOUBLECOMPLEX (doublecomplex *)
15: #define CASTDOUBLECOMPLEXSTAR (doublecomplex **)
16: #include <superlu_zdefs.h>
17: #define LUstructInit zLUstructInit
18: #define ScalePermstructInit zScalePermstructInit
19: #define ScalePermstructFree zScalePermstructFree
20: #define LUstructFree zLUstructFree
21: #define Destroy_LU zDestroy_LU
22: #define ScalePermstruct_t zScalePermstruct_t
23: #define LUstruct_t zLUstruct_t
24: #define SOLVEstruct_t zSOLVEstruct_t
25: #define SolveFinalize zSolveFinalize
26: #define pGetDiagU pzGetDiagU
27: #define pgssvx pzgssvx
28: #define allocateA_dist zallocateA_dist
29: #define Create_CompRowLoc_Matrix_dist zCreate_CompRowLoc_Matrix_dist
30: #define SLU SLU_Z
31: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
32: #define DeAllocLlu_3d zDeAllocLlu_3d
33: #define DeAllocGlu_3d zDeAllocGlu_3d
34: #define Destroy_A3d_gathered_on_2d zDestroy_A3d_gathered_on_2d
35: #define pgssvx3d pzgssvx3d
36: #endif
37: #elif PetscDefined(USE_REAL_SINGLE)
38: #define CASTDOUBLECOMPLEX
39: #define CASTDOUBLECOMPLEXSTAR
40: #include <superlu_sdefs.h>
41: #define LUstructInit sLUstructInit
42: #define ScalePermstructInit sScalePermstructInit
43: #define ScalePermstructFree sScalePermstructFree
44: #define LUstructFree sLUstructFree
45: #define Destroy_LU sDestroy_LU
46: #define ScalePermstruct_t sScalePermstruct_t
47: #define LUstruct_t sLUstruct_t
48: #define SOLVEstruct_t sSOLVEstruct_t
49: #define SolveFinalize sSolveFinalize
50: #define pGetDiagU psGetDiagU
51: #define pgssvx psgssvx
52: #define allocateA_dist sallocateA_dist
53: #define Create_CompRowLoc_Matrix_dist sCreate_CompRowLoc_Matrix_dist
54: #define SLU SLU_S
55: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
56: #define DeAllocLlu_3d sDeAllocLlu_3d
57: #define DeAllocGlu_3d sDeAllocGlu_3d
58: #define Destroy_A3d_gathered_on_2d sDestroy_A3d_gathered_on_2d
59: #define pgssvx3d psgssvx3d
60: #endif
61: #else
62: #define CASTDOUBLECOMPLEX
63: #define CASTDOUBLECOMPLEXSTAR
64: #include <superlu_ddefs.h>
65: #define LUstructInit dLUstructInit
66: #define ScalePermstructInit dScalePermstructInit
67: #define ScalePermstructFree dScalePermstructFree
68: #define LUstructFree dLUstructFree
69: #define Destroy_LU dDestroy_LU
70: #define ScalePermstruct_t dScalePermstruct_t
71: #define LUstruct_t dLUstruct_t
72: #define SOLVEstruct_t dSOLVEstruct_t
73: #define SolveFinalize dSolveFinalize
74: #define pGetDiagU pdGetDiagU
75: #define pgssvx pdgssvx
76: #define allocateA_dist dallocateA_dist
77: #define Create_CompRowLoc_Matrix_dist dCreate_CompRowLoc_Matrix_dist
78: #define SLU SLU_D
79: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
80: #define DeAllocLlu_3d dDeAllocLlu_3d
81: #define DeAllocGlu_3d dDeAllocGlu_3d
82: #define Destroy_A3d_gathered_on_2d dDestroy_A3d_gathered_on_2d
83: #define pgssvx3d pdgssvx3d
84: #endif
85: #endif
86: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
87: #include <superlu_sdefs.h>
88: #endif
89: EXTERN_C_END
90: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
92: typedef struct {
93: int_t nprow, npcol, *row, *col;
94: gridinfo_t grid;
95: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
96: PetscBool use3d;
97: int_t npdep; /* replication factor, must be power of two */
98: gridinfo3d_t grid3d;
99: #endif
100: superlu_dist_options_t options;
101: SuperMatrix A_sup;
102: ScalePermstruct_t ScalePermstruct;
103: LUstruct_t LUstruct;
104: int StatPrint;
105: SOLVEstruct_t SOLVEstruct;
106: fact_t FactPattern;
107: MPI_Comm comm_superlu;
108: PetscScalar *val;
109: PetscBool matsolve_iscalled, matmatsolve_iscalled;
110: PetscBool CleanUpSuperLU_Dist; /* Flag to clean up (non-global) SuperLU objects during Destroy */
111: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
112: sScalePermstruct_t sScalePermstruct;
113: sLUstruct_t sLUstruct;
114: sSOLVEstruct_t sSOLVEstruct;
115: float *sval;
116: PetscBool singleprecision; /* use single precision SuperLU_DIST from double precision PETSc */
117: float *sbptr;
118: #endif
119: } Mat_SuperLU_DIST;
121: static PetscErrorCode MatSuperluDistGetDiagU_SuperLU_DIST(Mat F, PetscScalar *diagU)
122: {
123: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)F->data;
125: PetscFunctionBegin;
126: PetscCallExternalVoid("SuperLU_DIST:pGetDiagU", pGetDiagU(F->rmap->N, &lu->LUstruct, &lu->grid, CASTDOUBLECOMPLEX diagU));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: /*@
131: MatSuperluDistGetDiagU - Get the diagonal of the U factor of a `MATSOLVERSUPERLU_DIST` factored matrix.
133: Collective
135: Input Parameter:
136: . F - the factored matrix returned by `MatGetFactor()` with `MATSOLVERSUPERLU_DIST`
138: Output Parameter:
139: . diagU - array of the U factor diagonal; must be preallocated by the caller
141: Level: advanced
143: .seealso: `Mat`, `MATSOLVERSUPERLU_DIST`, `MatGetFactor()`, `MatLUFactorSymbolic()`, `MatLUFactorNumeric()`
144: @*/
145: PetscErrorCode MatSuperluDistGetDiagU(Mat F, PetscScalar *diagU)
146: {
147: PetscFunctionBegin;
149: PetscTryMethod(F, "MatSuperluDistGetDiagU_C", (Mat, PetscScalar *), (F, diagU));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: /* This allows reusing the Superlu_DIST communicator and grid when only a single SuperLU_DIST matrix is used at a time */
154: typedef struct {
155: MPI_Comm comm;
156: PetscBool busy;
157: gridinfo_t grid;
158: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
159: PetscBool use3d;
160: gridinfo3d_t grid3d;
161: #endif
162: } PetscSuperLU_DIST;
164: static PetscMPIInt Petsc_Superlu_dist_keyval = MPI_KEYVAL_INVALID;
166: PETSC_EXTERN PetscMPIInt MPIAPI Petsc_Superlu_dist_keyval_DeleteFn(MPI_Comm comm, PetscMPIInt keyval, void *attr_val, void *extra_state)
167: {
168: PetscSuperLU_DIST *context = (PetscSuperLU_DIST *)attr_val;
170: PetscFunctionBegin;
171: PetscCheckReturnMPI(keyval == Petsc_Superlu_dist_keyval, PETSC_COMM_SELF, PETSC_ERR_ARG_CORRUPT, "Unexpected keyval");
172: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
173: if (context->use3d) {
174: PetscCallExternalVoid("SuperLU_DIST:superlu_gridexit3d", superlu_gridexit3d(&context->grid3d));
175: } else
176: #endif
177: PetscCallExternalVoid("SuperLU_DIST:superlu_gridexit", superlu_gridexit(&context->grid));
178: PetscCallMPIReturnMPI(MPI_Comm_free(&context->comm));
179: PetscCallReturnMPI(PetscFree(context));
180: PetscFunctionReturn(MPI_SUCCESS);
181: }
183: /*
184: Performs MPI_Comm_free_keyval() on Petsc_Superlu_dist_keyval but keeps the global variable for
185: users who do not destroy all PETSc objects before PetscFinalize().
187: The value Petsc_Superlu_dist_keyval is retained so that Petsc_Superlu_dist_keyval_DeleteFn()
188: can still check that the keyval associated with the MPI communicator is correct when the MPI
189: communicator is destroyed.
191: This is called in PetscFinalize()
192: */
193: static PetscErrorCode Petsc_Superlu_dist_keyval_free(void)
194: {
195: PetscMPIInt Petsc_Superlu_dist_keyval_temp = Petsc_Superlu_dist_keyval;
197: PetscFunctionBegin;
198: PetscCallMPI(MPI_Comm_free_keyval(&Petsc_Superlu_dist_keyval_temp));
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode MatDestroy_SuperLU_DIST(Mat A)
203: {
204: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)A->data;
206: PetscFunctionBegin;
207: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
208: PetscCall(PetscFree(lu->sbptr));
209: #endif
210: if (lu->CleanUpSuperLU_Dist) {
211: /* Deallocate SuperLU_DIST storage */
212: PetscCallExternalVoid("SuperLU_DIST:Destroy_CompRowLoc_Matrix_dist", Destroy_CompRowLoc_Matrix_dist(&lu->A_sup));
213: if (lu->options.SolveInitialized) {
214: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
215: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", sSolveFinalize(&lu->options, &lu->sSOLVEstruct));
216: else
217: #endif
218: PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", SolveFinalize(&lu->options, &lu->SOLVEstruct));
219: }
220: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
221: if (lu->use3d) {
222: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
223: if (lu->singleprecision) {
224: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", sDestroy_LU(A->cmap->N, &lu->grid3d.grid2d, &lu->sLUstruct));
225: PetscCallExternalVoid("SuperLU_DIST:Destroy_A3d_gathered_on_2d", sDestroy_A3d_gathered_on_2d(&lu->sSOLVEstruct, &lu->grid3d));
226: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructFree", sScalePermstructFree(&lu->sScalePermstruct));
227: PetscCallExternalVoid("SuperLU_DIST:LUstructFree", sLUstructFree(&lu->sLUstruct));
228: } else
229: #endif
230: {
231: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", Destroy_LU(A->cmap->N, &lu->grid3d.grid2d, &lu->LUstruct));
232: PetscCallExternalVoid("SuperLU_DIST:Destroy_A3d_gathered_on_2d", Destroy_A3d_gathered_on_2d(&lu->SOLVEstruct, &lu->grid3d));
233: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructFree", ScalePermstructFree(&lu->ScalePermstruct));
234: PetscCallExternalVoid("SuperLU_DIST:LUstructFree", LUstructFree(&lu->LUstruct));
235: }
236: } else
237: #endif
238: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
239: if (lu->singleprecision) {
240: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", sDestroy_LU(A->cmap->N, &lu->grid, &lu->sLUstruct));
241: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructFree", sScalePermstructFree(&lu->sScalePermstruct));
242: PetscCallExternalVoid("SuperLU_DIST:LUstructFree", sLUstructFree(&lu->sLUstruct));
243: } else
244: #endif
245: {
246: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", Destroy_LU(A->cmap->N, &lu->grid, &lu->LUstruct));
247: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructFree", ScalePermstructFree(&lu->ScalePermstruct));
248: PetscCallExternalVoid("SuperLU_DIST:LUstructFree", LUstructFree(&lu->LUstruct));
249: }
250: /* Release the SuperLU_DIST process grid only if the matrix has its own copy, that is it is not in the communicator context */
251: if (lu->comm_superlu) {
252: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
253: if (lu->use3d) {
254: PetscCallExternalVoid("SuperLU_DIST:superlu_gridexit3d", superlu_gridexit3d(&lu->grid3d));
255: } else
256: #endif
257: PetscCallExternalVoid("SuperLU_DIST:superlu_gridexit", superlu_gridexit(&lu->grid));
258: }
259: }
260: /*
261: * We always need to release the communicator that was created in MatGetFactor_aij_superlu_dist.
262: * lu->CleanUpSuperLU_Dist was turned on in MatLUFactorSymbolic_SuperLU_DIST. There are some use
263: * cases where we only create a matrix but do not solve mat. In these cases, lu->CleanUpSuperLU_Dist
264: * is off, and the communicator was not released or marked as "not busy " in the old code.
265: * Here we try to release comm regardless.
266: */
267: if (lu->comm_superlu) {
268: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &lu->comm_superlu));
269: } else {
270: PetscSuperLU_DIST *context;
271: MPI_Comm comm;
272: PetscMPIInt iflg;
274: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
275: PetscCallMPI(MPI_Comm_get_attr(comm, Petsc_Superlu_dist_keyval, &context, &iflg));
276: if (iflg) context->busy = PETSC_FALSE;
277: }
279: PetscCall(PetscFree(A->data));
280: /* clear composed functions */
281: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
282: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSuperluDistGetDiagU_C", NULL));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode MatSolve_SuperLU_DIST(Mat A, Vec b_mpi, Vec x)
287: {
288: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)A->data;
289: PetscInt m = A->rmap->n;
290: SuperLUStat_t stat;
291: PetscReal berr[1];
292: PetscScalar *bptr = NULL;
293: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
294: float sberr[1];
295: #endif
296: int info; /* SuperLU_Dist info code is ALWAYS an int, even with long long indices */
297: static PetscBool cite = PETSC_FALSE;
299: PetscFunctionBegin;
300: PetscCheck(lu->options.Fact == FACTORED, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "SuperLU_DIST options.Fact must equal FACTORED");
301: PetscCall(PetscCitationsRegister("@article{lidemmel03,\n author = {Xiaoye S. Li and James W. Demmel},\n title = {{SuperLU_DIST}: A Scalable Distributed-Memory Sparse Direct\n Solver for Unsymmetric Linear Systems},\n journal = {ACM "
302: "Trans. Mathematical Software},\n volume = {29},\n number = {2},\n pages = {110-140},\n year = 2003\n}\n",
303: &cite));
305: if (lu->options.SolveInitialized && !lu->matsolve_iscalled) {
306: /* see comments in MatMatSolve() */
307: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
308: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", sSolveFinalize(&lu->options, &lu->sSOLVEstruct));
309: else
310: #endif
311: PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", SolveFinalize(&lu->options, &lu->SOLVEstruct));
312: lu->options.SolveInitialized = NO;
313: }
314: PetscCall(VecCopy(b_mpi, x));
315: PetscCall(VecGetArray(x, &bptr));
316: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
317: if (lu->singleprecision) {
318: PetscInt n;
319: PetscCall(VecGetLocalSize(x, &n));
320: if (!lu->sbptr) PetscCall(PetscMalloc1(n, &lu->sbptr));
321: for (PetscInt i = 0; i < n; i++) lu->sbptr[i] = PetscRealPart(bptr[i]); /* PetscRealPart() is a no-op to allow compiling with complex */
322: }
323: #endif
325: PetscCallExternalVoid("SuperLU_DIST:PStatInit", PStatInit(&stat)); /* Initialize the statistics variables. */
326: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
327: if (lu->use3d) {
328: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
329: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:pgssvx3d", psgssvx3d(&lu->options, &lu->A_sup, &lu->sScalePermstruct, lu->sbptr, (int)m, 1, &lu->grid3d, &lu->sLUstruct, &lu->sSOLVEstruct, sberr, &stat, &info));
330: else
331: #endif
332: PetscCallExternalVoid("SuperLU_DIST:pgssvx3d", pgssvx3d(&lu->options, &lu->A_sup, &lu->ScalePermstruct, CASTDOUBLECOMPLEX bptr, (int)m, 1, &lu->grid3d, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info));
333: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "pgssvx3d fails, info: %d", info);
334: } else
335: #endif
336: {
337: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
338: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:pgssvx", psgssvx(&lu->options, &lu->A_sup, &lu->sScalePermstruct, lu->sbptr, (int)m, 1, &lu->grid, &lu->sLUstruct, &lu->sSOLVEstruct, sberr, &stat, &info));
339: else
340: #endif
341: PetscCallExternalVoid("SuperLU_DIST:pgssvx", pgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, CASTDOUBLECOMPLEX bptr, (int)m, 1, &lu->grid, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info));
342: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "pgssvx fails, info: %d", info);
343: }
344: if (lu->options.PrintStat) PetscCallExternalVoid("SuperLU_DIST:PStatPrint", PStatPrint(&lu->options, &stat, &lu->grid)); /* Print the statistics. */
345: PetscCallExternalVoid("SuperLU_DIST:PStatFree", PStatFree(&stat));
346: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
347: if (lu->singleprecision) {
348: PetscInt n;
349: PetscCall(VecGetLocalSize(x, &n));
350: for (PetscInt i = 0; i < n; i++) bptr[i] = lu->sbptr[i];
351: }
352: #endif
353: PetscCall(VecRestoreArray(x, &bptr));
354: lu->matsolve_iscalled = PETSC_TRUE;
355: lu->matmatsolve_iscalled = PETSC_FALSE;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode MatMatSolve_SuperLU_DIST(Mat A, Mat B_mpi, Mat X)
360: {
361: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)A->data;
362: PetscInt m = A->rmap->n, nrhs;
363: SuperLUStat_t stat;
364: PetscReal berr[1];
365: PetscScalar *bptr;
366: int info; /* SuperLU_Dist info code is ALWAYS an int, even with long long indices */
367: PetscBool flg;
369: PetscFunctionBegin;
370: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
371: PetscCheck(!lu->singleprecision, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for -pc_precision single");
372: #endif
373: PetscCheck(lu->options.Fact == FACTORED, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "SuperLU_DIST options.Fact must equal FACTORED");
374: PetscCall(PetscObjectTypeCompareAny((PetscObject)B_mpi, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
375: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix B must be MATDENSE matrix");
376: if (X != B_mpi) {
377: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSEQDENSE, MATMPIDENSE, NULL));
378: PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix X must be MATDENSE matrix");
379: }
381: if (lu->options.SolveInitialized && !lu->matmatsolve_iscalled) {
382: /* communication pattern of SOLVEstruct is unlikely created for matmatsolve,
383: thus destroy it and create a new SOLVEstruct.
384: Otherwise it may result in memory corruption or incorrect solution
385: See src/mat/tests/ex125.c */
386: PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", SolveFinalize(&lu->options, &lu->SOLVEstruct));
387: lu->options.SolveInitialized = NO;
388: }
389: if (X != B_mpi) PetscCall(MatCopy(B_mpi, X, SAME_NONZERO_PATTERN));
391: PetscCall(MatGetSize(B_mpi, NULL, &nrhs));
393: PetscCallExternalVoid("SuperLU_DIST:PStatInit", PStatInit(&stat)); /* Initialize the statistics variables. */
394: PetscCall(MatDenseGetArray(X, &bptr));
396: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
397: if (lu->use3d) PetscCallExternalVoid("SuperLU_DIST:pgssvx3d", pgssvx3d(&lu->options, &lu->A_sup, &lu->ScalePermstruct, CASTDOUBLECOMPLEX bptr, (int)m, (int)nrhs, &lu->grid3d, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info));
398: else
399: #endif
400: PetscCallExternalVoid("SuperLU_DIST:pgssvx", pgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, CASTDOUBLECOMPLEX bptr, (int)m, (int)nrhs, &lu->grid, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &info));
401: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "pdgssvx fails, info: %d", info);
402: PetscCall(MatDenseRestoreArray(X, &bptr));
404: if (lu->options.PrintStat) PetscCallExternalVoid("SuperLU_DIST:PStatPrint", PStatPrint(&lu->options, &stat, &lu->grid)); /* Print the statistics. */
405: PetscCallExternalVoid("SuperLU_DIST:PStatFree", PStatFree(&stat));
406: lu->matsolve_iscalled = PETSC_FALSE;
407: lu->matmatsolve_iscalled = PETSC_TRUE;
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: /*
412: input:
413: F: numeric Cholesky factor
414: output:
415: nneg: total number of negative pivots
416: nzero: total number of zero pivots
417: npos: (global dimension of F) - nneg - nzero
418: */
419: static PetscErrorCode MatGetInertia_SuperLU_DIST(Mat F, PetscInt *nneg, PetscInt *nzero, PetscInt *npos)
420: {
421: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)F->data;
422: PetscScalar *diagU = NULL;
423: PetscInt M, i, neg = 0, zero = 0, pos = 0;
424: PetscReal r;
426: PetscFunctionBegin;
427: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
428: PetscCheck(!lu->singleprecision, PetscObjectComm((PetscObject)F), PETSC_ERR_SUP, "Not for -pc_precision single");
429: #endif
430: PetscCheck(F->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix factor F is not assembled");
431: PetscCheck(lu->options.RowPerm == NOROWPERM, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Must set NOROWPERM");
432: PetscCall(MatGetSize(F, &M, NULL));
433: PetscCall(PetscMalloc1(M, &diagU));
434: PetscCall(MatSuperluDistGetDiagU(F, diagU));
435: for (i = 0; i < M; i++) {
436: #if PetscDefined(USE_COMPLEX)
437: r = PetscAbsReal(PetscImaginaryPart(diagU[i]));
438: PetscCheck(r < 1000 * PETSC_MACHINE_EPSILON, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "diagU[%" PetscInt_FMT "]=%g + i %g is non-real", i, (double)PetscRealPart(diagU[i]), (double)PetscImaginaryPart(diagU[i]));
439: r = PetscRealPart(diagU[i]);
440: #else
441: r = diagU[i];
442: #endif
443: if (r > 0) {
444: pos++;
445: } else if (r < 0) {
446: neg++;
447: } else zero++;
448: }
450: PetscCall(PetscFree(diagU));
451: if (nneg) *nneg = neg;
452: if (nzero) *nzero = zero;
453: if (npos) *npos = pos;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode MatLUFactorNumeric_SuperLU_DIST(Mat F, Mat A, const MatFactorInfo *info)
458: {
459: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)F->data;
460: Mat Aloc;
461: const PetscScalar *av;
462: const PetscInt *ai = NULL, *aj = NULL;
463: PetscInt nz, unused;
464: int sinfo; /* SuperLU_Dist info flag is always an int even with long long indices */
465: SuperLUStat_t stat;
466: PetscReal *berr = 0;
467: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
468: float *sberr = 0;
469: #endif
470: PetscBool ismpiaij, isseqaij, flg;
472: PetscFunctionBegin;
473: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
474: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
475: if (ismpiaij) PetscCall(MatMPIAIJGetLocalMat(A, MAT_INITIAL_MATRIX, &Aloc));
476: else if (isseqaij) {
477: PetscCall(PetscObjectReference((PetscObject)A));
478: Aloc = A;
479: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for type %s", ((PetscObject)A)->type_name);
481: PetscCall(MatGetRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &unused, &ai, &aj, &flg));
482: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "GetRowIJ failed");
483: PetscCall(MatSeqAIJGetArrayRead(Aloc, &av));
484: nz = ai[Aloc->rmap->n];
486: /* Allocations for A_sup */
487: if (lu->options.Fact == DOFACT) { /* first numeric factorization */
488: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
489: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:allocateA_dist", sallocateA_dist(Aloc->rmap->n, nz, &lu->sval, &lu->col, &lu->row));
490: else
491: #endif
492: PetscCallExternalVoid("SuperLU_DIST:allocateA_dist", allocateA_dist(Aloc->rmap->n, nz, CASTDOUBLECOMPLEXSTAR & lu->val, &lu->col, &lu->row));
493: } else { /* successive numeric factorization, sparsity pattern and perm_c are reused. */
494: if (lu->FactPattern == SamePattern_SameRowPerm) {
495: lu->options.Fact = SamePattern_SameRowPerm; /* matrix has similar numerical values */
496: } else if (lu->FactPattern == SamePattern) {
497: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
498: if (lu->use3d) {
499: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
500: if (lu->singleprecision) {
501: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", sDestroy_LU(A->cmap->N, &lu->grid3d.grid2d, &lu->sLUstruct));
502: PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", sSolveFinalize(&lu->options, &lu->sSOLVEstruct));
503: } else
504: #endif
505: {
506: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", Destroy_LU(A->cmap->N, &lu->grid3d.grid2d, &lu->LUstruct));
507: PetscCallExternalVoid("SuperLU_DIST:SolveFinalize", SolveFinalize(&lu->options, &lu->SOLVEstruct));
508: }
509: } else
510: #endif
511: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
512: if (lu->singleprecision)
513: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", sDestroy_LU(A->rmap->N, &lu->grid, &lu->sLUstruct));
514: else
515: #endif
516: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", Destroy_LU(A->rmap->N, &lu->grid, &lu->LUstruct));
517: lu->options.Fact = SamePattern;
518: } else if (lu->FactPattern == DOFACT) {
519: PetscCallExternalVoid("SuperLU_DIST:Destroy_CompRowLoc_Matrix_dist", Destroy_CompRowLoc_Matrix_dist(&lu->A_sup));
520: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
521: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", sDestroy_LU(A->rmap->N, &lu->grid, &lu->sLUstruct));
522: else
523: #endif
524: PetscCallExternalVoid("SuperLU_DIST:Destroy_LU", Destroy_LU(A->rmap->N, &lu->grid, &lu->LUstruct));
525: lu->options.Fact = DOFACT;
526: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
527: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:allocateA_dist", sallocateA_dist(Aloc->rmap->n, nz, &lu->sval, &lu->col, &lu->row));
528: else
529: #endif
530: PetscCallExternalVoid("SuperLU_DIST:allocateA_dist", allocateA_dist(Aloc->rmap->n, nz, CASTDOUBLECOMPLEXSTAR & lu->val, &lu->col, &lu->row));
531: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "options.Fact must be one of SamePattern SamePattern_SameRowPerm DOFACT");
532: }
534: /* Copy AIJ matrix to superlu_dist matrix */
535: PetscCall(PetscArraycpy(lu->row, ai, Aloc->rmap->n + 1));
536: PetscCall(PetscArraycpy(lu->col, aj, nz));
537: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
538: if (lu->singleprecision)
539: for (PetscInt i = 0; i < nz; i++) lu->sval[i] = PetscRealPart(av[i]); /* PetscRealPart() is a no-op to allow compiling with complex */
540: else
541: #endif
542: PetscCall(PetscArraycpy(lu->val, av, nz));
543: PetscCall(MatRestoreRowIJ(Aloc, 0, PETSC_FALSE, PETSC_FALSE, &unused, &ai, &aj, &flg));
544: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_SUP, "RestoreRowIJ failed");
545: PetscCall(MatSeqAIJRestoreArrayRead(Aloc, &av));
546: PetscCall(MatDestroy(&Aloc));
548: /* Create and setup A_sup */
549: if (lu->options.Fact == DOFACT) {
550: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
551: if (lu->singleprecision)
552: PetscCallExternalVoid("SuperLU_DIST:Create_CompRowLoc_Matrix_dist", sCreate_CompRowLoc_Matrix_dist(&lu->A_sup, A->rmap->N, A->cmap->N, nz, A->rmap->n, A->rmap->rstart, lu->sval, lu->col, lu->row, SLU_NR_loc, SLU_S, SLU_GE));
553: else
554: #endif
555: PetscCallExternalVoid("SuperLU_DIST:Create_CompRowLoc_Matrix_dist", Create_CompRowLoc_Matrix_dist(&lu->A_sup, A->rmap->N, A->cmap->N, nz, A->rmap->n, A->rmap->rstart, CASTDOUBLECOMPLEX lu->val, lu->col, lu->row, SLU_NR_loc, SLU, SLU_GE));
556: }
558: /* Factor the matrix. */
559: PetscCallExternalVoid("SuperLU_DIST:PStatInit", PStatInit(&stat)); /* Initialize the statistics variables. */
560: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
561: if (lu->use3d) {
562: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
563: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:pgssvx3d", psgssvx3d(&lu->options, &lu->A_sup, &lu->sScalePermstruct, 0, (int)A->rmap->n, 0, &lu->grid3d, &lu->sLUstruct, &lu->sSOLVEstruct, sberr, &stat, &sinfo));
564: else
565: #endif
566: PetscCallExternalVoid("SuperLU_DIST:pgssvx3d", pgssvx3d(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, (int)A->rmap->n, 0, &lu->grid3d, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &sinfo));
567: } else
568: #endif
569: {
570: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
571: if (lu->singleprecision) PetscCallExternalVoid("SuperLU_DIST:pgssvx", psgssvx(&lu->options, &lu->A_sup, &lu->sScalePermstruct, 0, (int)A->rmap->n, 0, &lu->grid, &lu->sLUstruct, &lu->sSOLVEstruct, sberr, &stat, &sinfo));
572: else
573: #endif
574: PetscCallExternalVoid("SuperLU_DIST:pgssvx", pgssvx(&lu->options, &lu->A_sup, &lu->ScalePermstruct, 0, (int)A->rmap->n, 0, &lu->grid, &lu->LUstruct, &lu->SOLVEstruct, berr, &stat, &sinfo));
575: }
576: if (sinfo > 0) {
577: PetscCheck(!A->erroriffailure, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot in row %d", sinfo);
578: if (sinfo <= lu->A_sup.ncol) {
579: F->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
580: PetscCall(PetscInfo(F, "U(i,i) is exactly zero, i= %d\n", sinfo));
581: } else if (sinfo > lu->A_sup.ncol) {
582: /*
583: number of bytes allocated when memory allocation
584: failure occurred, plus A->ncol.
585: */
586: F->factorerrortype = MAT_FACTOR_OUTMEMORY;
587: PetscCall(PetscInfo(F, "Number of bytes allocated when memory allocation fails %d\n", sinfo));
588: }
589: } else PetscCheck(sinfo >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "info = %d, argument in p*gssvx() had an illegal value", sinfo);
591: if (lu->options.PrintStat) PetscCallExternalVoid("SuperLU_DIST:PStatPrint", PStatPrint(&lu->options, &stat, &lu->grid)); /* Print the statistics. */
592: PetscCallExternalVoid("SuperLU_DIST:PStatFree", PStatFree(&stat));
593: F->assembled = PETSC_TRUE;
594: F->preallocated = PETSC_TRUE;
595: lu->options.Fact = FACTORED; /* The factored form of A is supplied. Local option used by this func. only */
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: /* Note the PETSc r and c permutations are ignored */
600: static PetscErrorCode MatLUFactorSymbolic_SuperLU_DIST(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
601: {
602: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)F->data;
603: PetscInt M = A->rmap->N, N = A->cmap->N, indx;
604: PetscMPIInt size, iflg;
605: PetscBool flg, set;
606: const char *colperm[] = {"NATURAL", "MMD_AT_PLUS_A", "MMD_ATA", "METIS_AT_PLUS_A", "PARMETIS"};
607: const char *rowperm[] = {"NOROWPERM", "LargeDiag_MC64", "LargeDiag_AWPM", "MY_PERMR"};
608: const char *factPattern[] = {"SamePattern", "SamePattern_SameRowPerm", "DOFACT"};
609: MPI_Comm comm;
610: PetscSuperLU_DIST *context = NULL;
611: PetscPrecision precision = PETSC_PRECISION_INVALID;
613: PetscFunctionBegin;
614: /* Set options to F */
615: PetscCall(PetscObjectGetComm((PetscObject)F, &comm));
616: PetscCallMPI(MPI_Comm_size(comm, &size));
618: PetscOptionsBegin(PetscObjectComm((PetscObject)F), ((PetscObject)F)->prefix, "SuperLU_Dist Options", "Mat");
619: PetscCall(PetscOptionsBool("-mat_superlu_dist_equil", "Equilibrate matrix", "None", lu->options.Equil ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
620: if (set && !flg) lu->options.Equil = NO;
622: PetscCall(PetscOptionsEnum("-pc_precision", "Precision used by SuperLU_DIST", "MATSOLVERSUPERLU_DIST", PetscPrecisionTypes, (PetscEnum)precision, (PetscEnum *)&precision, &flg));
623: if (flg) {
624: PetscCheck(precision == PETSC_PRECISION_SINGLE || precision == PETSC_PRECISION_DOUBLE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER_INPUT, "-pc_precision only accepts single or double as option for SuperLU_DIST");
625: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
626: lu->singleprecision = (PetscBool)(precision == PETSC_PRECISION_SINGLE); // It also implies PetscReal is not single; not merely SuperLU_DIST is running in single
627: #endif
628: }
630: PetscCall(PetscOptionsEList("-mat_superlu_dist_rowperm", "Row permutation", "None", rowperm, 4, rowperm[1], &indx, &flg));
631: if (flg) {
632: switch (indx) {
633: case 0:
634: lu->options.RowPerm = NOROWPERM;
635: break;
636: case 1:
637: lu->options.RowPerm = LargeDiag_MC64;
638: break;
639: case 2:
640: lu->options.RowPerm = LargeDiag_AWPM;
641: break;
642: case 3:
643: lu->options.RowPerm = MY_PERMR;
644: break;
645: default:
646: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown row permutation");
647: }
648: }
650: PetscCall(PetscOptionsEList("-mat_superlu_dist_colperm", "Column permutation", "None", colperm, 5, colperm[3], &indx, &flg));
651: if (flg) {
652: switch (indx) {
653: case 0:
654: lu->options.ColPerm = NATURAL;
655: break;
656: case 1:
657: lu->options.ColPerm = MMD_AT_PLUS_A;
658: break;
659: case 2:
660: lu->options.ColPerm = MMD_ATA;
661: break;
662: case 3:
663: lu->options.ColPerm = METIS_AT_PLUS_A;
664: break;
665: case 4:
666: lu->options.ColPerm = PARMETIS; /* only works for np>1 */
667: break;
668: default:
669: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown column permutation");
670: }
671: }
673: lu->options.ReplaceTinyPivot = NO;
674: PetscCall(PetscOptionsBool("-mat_superlu_dist_replacetinypivot", "Replace tiny pivots", "None", lu->options.ReplaceTinyPivot ? PETSC_TRUE : PETSC_FALSE, &flg, &set));
675: if (set && flg) lu->options.ReplaceTinyPivot = YES;
677: lu->options.ParSymbFact = NO;
678: PetscCall(PetscOptionsBool("-mat_superlu_dist_parsymbfact", "Parallel symbolic factorization", "None", PETSC_FALSE, &flg, &set));
679: if (set && flg && size > 1) {
680: #if PetscDefined(HAVE_PARMETIS)
681: lu->options.ParSymbFact = YES;
682: lu->options.ColPerm = PARMETIS; /* in v2.2, PARMETIS is forced for ParSymbFact regardless of user ordering setting */
683: #else
684: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "parsymbfact needs PARMETIS");
685: #endif
686: }
688: lu->FactPattern = SamePattern;
689: PetscCall(PetscOptionsEList("-mat_superlu_dist_fact", "Sparsity pattern for repeated matrix factorization", "None", factPattern, 3, factPattern[0], &indx, &flg));
690: if (flg) {
691: switch (indx) {
692: case 0:
693: lu->FactPattern = SamePattern;
694: break;
695: case 1:
696: lu->FactPattern = SamePattern_SameRowPerm;
697: break;
698: case 2:
699: lu->FactPattern = DOFACT;
700: break;
701: }
702: }
704: lu->options.IterRefine = NOREFINE;
705: PetscCall(PetscOptionsBool("-mat_superlu_dist_iterrefine", "Use iterative refinement", "None", lu->options.IterRefine == NOREFINE ? PETSC_FALSE : PETSC_TRUE, &flg, &set));
706: if (set && flg) lu->options.IterRefine = SLU_DOUBLE;
708: if (PetscLogPrintInfo) lu->options.PrintStat = YES;
709: else lu->options.PrintStat = NO;
710: PetscCall(PetscOptionsDeprecated("-mat_superlu_dist_statprint", "-mat_superlu_dist_printstat", "3.19", NULL));
711: PetscCall(PetscOptionsBool("-mat_superlu_dist_printstat", "Print factorization information", "None", (PetscBool)lu->options.PrintStat, (PetscBool *)&lu->options.PrintStat, NULL));
713: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(8, 0, 0)
714: lu->options.superlu_acc_offload = 1;
715: PetscCall(PetscOptionsBool("-mat_superlu_dist_gpuoffload", "Offload factorization onto the GPUs", "None", (PetscBool)lu->options.superlu_acc_offload, (PetscBool *)&lu->options.superlu_acc_offload, NULL));
716: #endif
718: PetscCallMPI(MPI_Comm_get_attr(comm, Petsc_Superlu_dist_keyval, &context, &iflg));
719: if (!iflg || context->busy) { /* additional options */
720: if (!iflg) {
721: PetscCall(PetscNew(&context));
722: context->busy = PETSC_TRUE;
723: PetscCallMPI(MPI_Comm_dup(comm, &context->comm));
724: PetscCallMPI(MPI_Comm_set_attr(comm, Petsc_Superlu_dist_keyval, context));
725: } else {
726: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)A), &lu->comm_superlu));
727: }
729: /* Default number of process columns and rows */
730: lu->nprow = (int_t)(0.5 + PetscSqrtReal((PetscReal)size));
731: if (!lu->nprow) lu->nprow = 1;
732: while (lu->nprow > 0) {
733: lu->npcol = (int_t)(size / lu->nprow);
734: if (size == lu->nprow * lu->npcol) break;
735: lu->nprow--;
736: }
737: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
738: lu->use3d = PETSC_FALSE;
739: lu->npdep = 1;
740: #endif
742: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
743: PetscCall(PetscOptionsBool("-mat_superlu_dist_3d", "Use SuperLU_DIST 3D distribution", "None", lu->use3d, &lu->use3d, NULL));
744: if (lu->use3d) {
745: PetscInt t;
746: PetscCall(PetscOptionsInt("-mat_superlu_dist_d", "Number of z entries in processor partition", "None", lu->npdep, (PetscInt *)&lu->npdep, NULL));
747: t = (PetscInt)PetscLog2Real((PetscReal)lu->npdep);
748: PetscCheck(PetscPowInt(2, t) == lu->npdep, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "-mat_superlu_dist_d %lld must be a power of 2", (long long)lu->npdep);
749: if (lu->npdep > 1) {
750: lu->nprow = (int_t)(0.5 + PetscSqrtReal((PetscReal)(size / lu->npdep)));
751: if (!lu->nprow) lu->nprow = 1;
752: while (lu->nprow > 0) {
753: lu->npcol = (int_t)(size / (lu->npdep * lu->nprow));
754: if (size == lu->nprow * lu->npcol * lu->npdep) break;
755: lu->nprow--;
756: }
757: }
758: }
759: #endif
760: PetscCall(PetscOptionsInt("-mat_superlu_dist_r", "Number rows in processor partition", "None", lu->nprow, (PetscInt *)&lu->nprow, NULL));
761: PetscCall(PetscOptionsInt("-mat_superlu_dist_c", "Number columns in processor partition", "None", lu->npcol, (PetscInt *)&lu->npcol, NULL));
762: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
763: PetscCheck(size == lu->nprow * lu->npcol * lu->npdep, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of processes %d must equal to nprow %lld * npcol %lld * npdep %lld", size, (long long)lu->nprow, (long long)lu->npcol, (long long)lu->npdep);
764: #else
765: PetscCheck(size == lu->nprow * lu->npcol, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of processes %d must equal to nprow %lld * npcol %lld", size, (long long)lu->nprow, (long long)lu->npcol);
766: #endif
767: /* end of adding additional options */
769: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
770: if (lu->use3d) {
771: PetscCallExternalVoid("SuperLU_DIST:superlu_gridinit3d", superlu_gridinit3d(context ? context->comm : lu->comm_superlu, (int)lu->nprow, (int)lu->npcol, (int)lu->npdep, &lu->grid3d));
772: if (context) {
773: context->grid3d = lu->grid3d;
774: context->use3d = lu->use3d;
775: }
776: } else {
777: #endif
778: PetscCallExternalVoid("SuperLU_DIST:superlu_gridinit", superlu_gridinit(context ? context->comm : lu->comm_superlu, (int)lu->nprow, (int)lu->npcol, &lu->grid));
779: if (context) context->grid = lu->grid;
780: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
781: }
782: #endif
783: PetscCall(PetscInfo(NULL, "Duplicating a communicator for SuperLU_DIST and calling superlu_gridinit()\n"));
784: if (iflg) {
785: PetscCall(PetscInfo(NULL, "Communicator attribute already in use so not saving communicator and SuperLU_DIST grid in communicator attribute \n"));
786: } else {
787: PetscCall(PetscInfo(NULL, "Storing communicator and SuperLU_DIST grid in communicator attribute\n"));
788: }
789: } else { /* (iflg && !context->busy) */
790: PetscCall(PetscInfo(NULL, "Reusing communicator and superlu_gridinit() for SuperLU_DIST from communicator attribute.\n"));
791: context->busy = PETSC_TRUE;
792: lu->grid = context->grid;
793: }
794: PetscOptionsEnd();
796: /* Initialize ScalePermstruct and LUstruct. */
797: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
798: if (lu->singleprecision) {
799: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructInit", sScalePermstructInit(M, N, &lu->sScalePermstruct));
800: PetscCallExternalVoid("SuperLU_DIST:LUstructInit", sLUstructInit(N, &lu->sLUstruct));
801: } else
802: #endif
803: {
804: PetscCallExternalVoid("SuperLU_DIST:ScalePermstructInit", ScalePermstructInit(M, N, &lu->ScalePermstruct));
805: PetscCallExternalVoid("SuperLU_DIST:LUstructInit", LUstructInit(N, &lu->LUstruct));
806: }
807: F->ops->lufactornumeric = MatLUFactorNumeric_SuperLU_DIST;
808: F->ops->solve = MatSolve_SuperLU_DIST;
809: F->ops->matsolve = MatMatSolve_SuperLU_DIST;
810: F->ops->getinertia = NULL;
812: if (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE) F->ops->getinertia = MatGetInertia_SuperLU_DIST;
813: lu->CleanUpSuperLU_Dist = PETSC_TRUE;
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: static PetscErrorCode MatCholeskyFactorSymbolic_SuperLU_DIST(Mat F, Mat A, IS r, const MatFactorInfo *info)
818: {
819: PetscFunctionBegin;
820: PetscCall(MatLUFactorSymbolic_SuperLU_DIST(F, A, r, r, info));
821: F->ops->choleskyfactornumeric = MatLUFactorNumeric_SuperLU_DIST;
822: PetscFunctionReturn(PETSC_SUCCESS);
823: }
825: static PetscErrorCode MatFactorGetSolverType_aij_superlu_dist(Mat A, MatSolverType *type)
826: {
827: PetscFunctionBegin;
828: *type = MATSOLVERSUPERLU_DIST;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: static PetscErrorCode MatView_Info_SuperLU_DIST(Mat A, PetscViewer viewer)
833: {
834: Mat_SuperLU_DIST *lu = (Mat_SuperLU_DIST *)A->data;
835: superlu_dist_options_t options;
837: PetscFunctionBegin;
838: /* check if matrix is superlu_dist type */
839: if (A->ops->solve != MatSolve_SuperLU_DIST) PetscFunctionReturn(PETSC_SUCCESS);
841: options = lu->options;
842: PetscCall(PetscViewerASCIIPrintf(viewer, "SuperLU_DIST run parameters:\n"));
843: /* would love to use superlu 'IFMT' macro but it looks like it's inconsistently applied, the
844: * format spec for int64_t is set to %d for whatever reason */
845: PetscCall(PetscViewerASCIIPrintf(viewer, " Process grid nprow %lld x npcol %lld \n", (long long)lu->nprow, (long long)lu->npcol));
846: #if PETSC_PKG_SUPERLU_DIST_VERSION_GE(9, 0, 0)
847: if (lu->use3d) PetscCall(PetscViewerASCIIPrintf(viewer, " Using 3d decomposition with npdep %lld \n", (long long)lu->npdep));
848: #endif
850: PetscCall(PetscViewerASCIIPrintf(viewer, " Equilibrate matrix %s \n", PetscBools[options.Equil != NO]));
851: PetscCall(PetscViewerASCIIPrintf(viewer, " Replace tiny pivots %s \n", PetscBools[options.ReplaceTinyPivot != NO]));
852: PetscCall(PetscViewerASCIIPrintf(viewer, " Use iterative refinement %s \n", PetscBools[options.IterRefine == SLU_DOUBLE]));
853: PetscCall(PetscViewerASCIIPrintf(viewer, " Processors in row %lld col partition %lld \n", (long long)lu->nprow, (long long)lu->npcol));
855: switch (options.RowPerm) {
856: case NOROWPERM:
857: PetscCall(PetscViewerASCIIPrintf(viewer, " Row permutation NOROWPERM\n"));
858: break;
859: case LargeDiag_MC64:
860: PetscCall(PetscViewerASCIIPrintf(viewer, " Row permutation LargeDiag_MC64\n"));
861: break;
862: case LargeDiag_AWPM:
863: PetscCall(PetscViewerASCIIPrintf(viewer, " Row permutation LargeDiag_AWPM\n"));
864: break;
865: case MY_PERMR:
866: PetscCall(PetscViewerASCIIPrintf(viewer, " Row permutation MY_PERMR\n"));
867: break;
868: default:
869: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown column permutation");
870: }
872: switch (options.ColPerm) {
873: case NATURAL:
874: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation NATURAL\n"));
875: break;
876: case MMD_AT_PLUS_A:
877: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation MMD_AT_PLUS_A\n"));
878: break;
879: case MMD_ATA:
880: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation MMD_ATA\n"));
881: break;
882: /* Even though this is called METIS, the SuperLU_DIST code sets this by default if PARMETIS is defined, not METIS */
883: case METIS_AT_PLUS_A:
884: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation METIS_AT_PLUS_A\n"));
885: break;
886: case PARMETIS:
887: PetscCall(PetscViewerASCIIPrintf(viewer, " Column permutation PARMETIS\n"));
888: break;
889: default:
890: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown column permutation");
891: }
893: PetscCall(PetscViewerASCIIPrintf(viewer, " Parallel symbolic factorization %s \n", PetscBools[options.ParSymbFact != NO]));
895: if (lu->FactPattern == SamePattern) {
896: PetscCall(PetscViewerASCIIPrintf(viewer, " Repeated factorization SamePattern\n"));
897: } else if (lu->FactPattern == SamePattern_SameRowPerm) {
898: PetscCall(PetscViewerASCIIPrintf(viewer, " Repeated factorization SamePattern_SameRowPerm\n"));
899: } else if (lu->FactPattern == DOFACT) {
900: PetscCall(PetscViewerASCIIPrintf(viewer, " Repeated factorization DOFACT\n"));
901: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unknown factorization pattern");
902: #if PetscDefined(HAVE_SUPERLU_DIST_SINGLE)
903: if (lu->singleprecision) PetscCall(PetscViewerASCIIPrintf(viewer, " Using SuperLU_DIST in single precision\n"));
904: #endif
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: static PetscErrorCode MatView_SuperLU_DIST(Mat A, PetscViewer viewer)
909: {
910: PetscBool isascii;
911: PetscViewerFormat format;
913: PetscFunctionBegin;
914: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
915: if (isascii) {
916: PetscCall(PetscViewerGetFormat(viewer, &format));
917: if (format == PETSC_VIEWER_ASCII_INFO) PetscCall(MatView_Info_SuperLU_DIST(A, viewer));
918: }
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: static PetscErrorCode MatGetFactor_aij_superlu_dist(Mat A, MatFactorType ftype, Mat *F)
923: {
924: Mat B;
925: Mat_SuperLU_DIST *lu;
926: PetscInt M = A->rmap->N, N = A->cmap->N;
927: PetscMPIInt size;
928: superlu_dist_options_t options;
930: PetscFunctionBegin;
931: /* Create the factorization matrix */
932: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
933: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, M, N));
934: PetscCall(PetscStrallocpy("superlu_dist", &((PetscObject)B)->type_name));
935: PetscCall(MatSetUp(B));
936: B->ops->getinfo = MatGetInfo_External;
937: B->ops->view = MatView_SuperLU_DIST;
938: B->ops->destroy = MatDestroy_SuperLU_DIST;
940: /* set_default_options_dist() sets the default input options to the following values:
941: options.Fact = DOFACT;
942: options.Equil = YES;
943: options.ParSymbFact = NO;
944: options.ColPerm = METIS_AT_PLUS_A;
945: options.RowPerm = LargeDiag_MC64;
946: options.ReplaceTinyPivot = NO;
947: options.IterRefine = SLU_DOUBLE;
948: options.Trans = NOTRANS;
949: options.SolveInitialized = NO; -hold the communication pattern used MatSolve() and MatMatSolve()
950: options.RefineInitialized = NO;
951: options.PrintStat = YES;
952: options.superlu_acc_offload = 1;
953: options.SymPattern = NO;
954: */
955: set_default_options_dist(&options);
957: B->trivialsymbolic = PETSC_TRUE;
958: if (ftype == MAT_FACTOR_LU) {
959: B->factortype = MAT_FACTOR_LU;
960: B->ops->lufactorsymbolic = MatLUFactorSymbolic_SuperLU_DIST;
961: } else {
962: B->factortype = MAT_FACTOR_CHOLESKY;
963: B->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SuperLU_DIST;
964: options.SymPattern = YES;
965: }
967: /* set solvertype */
968: PetscCall(PetscFree(B->solvertype));
969: PetscCall(PetscStrallocpy(MATSOLVERSUPERLU_DIST, &B->solvertype));
971: PetscCall(PetscNew(&lu));
972: B->data = lu;
973: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
975: lu->options = options;
976: lu->options.Fact = DOFACT;
977: lu->matsolve_iscalled = PETSC_FALSE;
978: lu->matmatsolve_iscalled = PETSC_FALSE;
980: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_aij_superlu_dist));
981: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSuperluDistGetDiagU_C", MatSuperluDistGetDiagU_SuperLU_DIST));
983: *F = B;
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_SuperLU_DIST(void)
988: {
989: PetscFunctionBegin;
990: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU_DIST, MATMPIAIJ, MAT_FACTOR_LU, MatGetFactor_aij_superlu_dist));
991: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU_DIST, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_aij_superlu_dist));
992: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU_DIST, MATMPIAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_aij_superlu_dist));
993: PetscCall(MatSolverTypeRegister(MATSOLVERSUPERLU_DIST, MATSEQAIJ, MAT_FACTOR_CHOLESKY, MatGetFactor_aij_superlu_dist));
994: if (Petsc_Superlu_dist_keyval == MPI_KEYVAL_INVALID) {
995: PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, Petsc_Superlu_dist_keyval_DeleteFn, &Petsc_Superlu_dist_keyval, NULL));
996: PetscCall(PetscRegisterFinalize(Petsc_Superlu_dist_keyval_free));
997: }
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: /*MC
1002: MATSOLVERSUPERLU_DIST - Parallel direct solver package for LU factorization
1004: Use `./configure --download-superlu_dist --download-parmetis --download-metis --download-ptscotch` to have PETSc installed with SuperLU_DIST
1006: Use `-pc_type lu` `-pc_factor_mat_solver_type superlu_dist` to use this direct solver
1008: Works with `MATAIJ` matrices
1010: Options Database Keys:
1011: + -mat_superlu_dist_r <n> - number of rows in processor partition
1012: . -mat_superlu_dist_c <n> - number of columns in processor partition
1013: . -mat_superlu_dist_3d - use 3d partition, requires SuperLU_DIST 7.2 or later
1014: . -mat_superlu_dist_d <n> - depth in 3d partition (valid only if `-mat_superlu_dist_3d`) is provided
1015: . -mat_superlu_dist_equil - equilibrate the matrix
1016: . -mat_superlu_dist_rowperm <NOROWPERM,LargeDiag_MC64,LargeDiag_AWPM,MY_PERMR> - row permutation
1017: . -mat_superlu_dist_colperm <NATURAL,MMD_AT_PLUS_A,MMD_ATA,METIS_AT_PLUS_A,PARMETIS> - column permutation
1018: . -mat_superlu_dist_replacetinypivot - replace tiny pivots
1019: . -mat_superlu_dist_fact <SamePattern> - (choose one of) `SamePattern`, `SamePattern_SameRowPerm`, `DOFACT`
1020: . -mat_superlu_dist_iterrefine - use iterative refinement
1021: . -mat_superlu_dist_printstat - print factorization information
1022: . -mat_superlu_dist_gpuoffload - offload factorization onto the GPUs, requires SuperLU_DIST 8.0.0 or later
1023: - -pc_precision single - use SuperLU_DIST single precision with PETSc double precision
1025: Level: beginner
1027: Note:
1028: If PETSc was configured with `--with-cuda` then this solver will use the GPUs by default.
1030: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`, `MatGetFactor()`
1031: M*/