Actual source code: lusol.c
1: /*
2: Provides an interface to the LUSOL package of ....
4: */
5: #include <../src/mat/impls/aij/seq/aij.h>
7: #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
8: #define LU1FAC lu1fac_
9: #define LU6SOL lu6sol_
10: #define M1PAGE m1page_
11: #define M5SETX m5setx_
12: #define M6RDEL m6rdel_
13: #elif !defined(PETSC_HAVE_FORTRAN_CAPS)
14: #define LU1FAC lu1fac
15: #define LU6SOL lu6sol
16: #define M1PAGE m1page
17: #define M5SETX m5setx
18: #define M6RDEL m6rdel
19: #endif
21: /*
22: Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require
23: */
24: PETSC_EXTERN void M1PAGE()
25: {
26: ;
27: }
28: PETSC_EXTERN void M5SETX()
29: {
30: ;
31: }
33: PETSC_EXTERN void M6RDEL()
34: {
35: ;
36: }
38: PETSC_EXTERN void LU1FAC(int *m, int *n, int *nnz, int *size, int *luparm, double *parmlu, double *data, int *indc, int *indr, int *rowperm, int *colperm, int *collen, int *rowlen, int *colstart, int *rowstart, int *rploc, int *cploc, int *rpinv, int *cpinv, double *w, int *inform);
40: PETSC_EXTERN void LU6SOL(int *mode, int *m, int *n, double *rhs, double *x, int *size, int *luparm, double *parmlu, double *data, int *indc, int *indr, int *rowperm, int *colperm, int *collen, int *rowlen, int *colstart, int *rowstart, int *inform);
42: typedef struct {
43: double *data;
44: int *indc;
45: int *indr;
47: int *ip;
48: int *iq;
49: int *lenc;
50: int *lenr;
51: int *locc;
52: int *locr;
53: int *iploc;
54: int *iqloc;
55: int *ipinv;
56: int *iqinv;
57: double *mnsw;
58: double *mnsv;
60: double elbowroom;
61: double luroom; /* Extra space allocated when factor fails */
62: double parmlu[30]; /* Input/output to LUSOL */
64: int n; /* Number of rows/columns in matrix */
65: int nz; /* Number of nonzeros */
66: int nnz; /* Number of nonzeros allocated for factors */
67: int luparm[30]; /* Input/output to LUSOL */
69: PetscBool CleanUpLUSOL;
71: } Mat_LUSOL;
73: /*
74: LUSOL input/Output Parameters (Description uses C-style indexes
76: Input parameters Typical value
77: luparm(0) = nout File number for printed messages. 6
78: luparm(1) = lprint Print level. 0
79: < 0 suppresses output.
80: = 0 gives error messages.
81: = 1 gives debug output from some of the
82: other routines in LUSOL.
83: >= 2 gives the pivot row and column and the
84: no. of rows and columns involved at
85: each elimination step in lu1fac.
86: luparm(2) = maxcol lu1fac: maximum number of columns 5
87: searched allowed in a Markowitz-type
88: search for the next pivot element.
89: For some of the factorization, the
90: number of rows searched is
91: maxrow = maxcol - 1.
93: Output parameters:
95: luparm(9) = inform Return code from last call to any LU routine.
96: luparm(10) = nsing No. of singularities marked in the
97: output array w(*).
98: luparm(11) = jsing Column index of last singularity.
99: luparm(12) = minlen Minimum recommended value for lena.
100: luparm(13) = maxlen ?
101: luparm(14) = nupdat No. of updates performed by the lu8 routines.
102: luparm(15) = nrank No. of nonempty rows of U.
103: luparm(16) = ndens1 No. of columns remaining when the density of
104: the matrix being factorized reached dens1.
105: luparm(17) = ndens2 No. of columns remaining when the density of
106: the matrix being factorized reached dens2.
107: luparm(18) = jumin The column index associated with dumin.
108: luparm(19) = numl0 No. of columns in initial L.
109: luparm(20) = lenl0 Size of initial L (no. of nonzeros).
110: luparm(21) = lenu0 Size of initial U.
111: luparm(22) = lenl Size of current L.
112: luparm(23) = lenu Size of current U.
113: luparm(24) = lrow Length of row file.
114: luparm(25) = ncp No. of compressions of LU data structures.
115: luparm(26) = mersum lu1fac: sum of Markowitz merit counts.
116: luparm(27) = nutri lu1fac: triangular rows in U.
117: luparm(28) = nltri lu1fac: triangular rows in L.
118: luparm(29) =
120: Input parameters Typical value
121: parmlu(0) = elmax1 Max multiplier allowed in L 10.0
122: during factor.
123: parmlu(1) = elmax2 Max multiplier allowed in L 10.0
124: during updates.
125: parmlu(2) = small Absolute tolerance for eps**0.8
126: treating reals as zero. IBM double: 3.0d-13
127: parmlu(3) = utol1 Absolute tol for flagging eps**0.66667
128: small diagonals of U. IBM double: 3.7d-11
129: parmlu(4) = utol2 Relative tol for flagging eps**0.66667
130: small diagonals of U. IBM double: 3.7d-11
131: parmlu(5) = uspace Factor limiting waste space in U. 3.0
132: In lu1fac, the row or column lists
133: are compressed if their length
134: exceeds uspace times the length of
135: either file after the last compression.
136: parmlu(6) = dens1 The density at which the Markowitz 0.3
137: strategy should search maxcol columns
138: and no rows.
139: parmlu(7) = dens2 the density at which the Markowitz 0.6
140: strategy should search only 1 column
141: or (preferably) use a dense LU for
142: all the remaining rows and columns.
144: Output parameters:
145: parmlu(9) = amax Maximum element in A.
146: parmlu(10) = elmax Maximum multiplier in current L.
147: parmlu(11) = umax Maximum element in current U.
148: parmlu(12) = dumax Maximum diagonal in U.
149: parmlu(13) = dumin Minimum diagonal in U.
150: parmlu(14) =
151: parmlu(15) =
152: parmlu(16) =
153: parmlu(17) =
154: parmlu(18) =
155: parmlu(19) = resid lu6sol: residual after solve with U or U'.
156: ...
157: parmlu(29) =
158: */
160: #define Factorization_Tolerance 1e-1
161: #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0)
162: #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */
164: static PetscErrorCode MatDestroy_LUSOL(Mat A)
165: {
166: Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;
168: PetscFunctionBegin;
169: if (lusol && lusol->CleanUpLUSOL) {
170: PetscCall(PetscFree(lusol->ip));
171: PetscCall(PetscFree(lusol->iq));
172: PetscCall(PetscFree(lusol->lenc));
173: PetscCall(PetscFree(lusol->lenr));
174: PetscCall(PetscFree(lusol->locc));
175: PetscCall(PetscFree(lusol->locr));
176: PetscCall(PetscFree(lusol->iploc));
177: PetscCall(PetscFree(lusol->iqloc));
178: PetscCall(PetscFree(lusol->ipinv));
179: PetscCall(PetscFree(lusol->iqinv));
180: PetscCall(PetscFree(lusol->mnsw));
181: PetscCall(PetscFree(lusol->mnsv));
182: PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
183: }
184: PetscCall(PetscFree(A->spptr));
185: PetscCall(MatDestroy_SeqAIJ(A));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode MatSolve_LUSOL(Mat A, Vec b, Vec x)
190: {
191: Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;
192: double *xx;
193: const double *bb;
194: int mode = 5;
195: int i, m, n, nnz, status;
197: PetscFunctionBegin;
198: PetscCall(VecGetArray(x, &xx));
199: PetscCall(VecGetArrayRead(b, &bb));
201: m = n = lusol->n;
202: nnz = lusol->nnz;
204: for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i];
206: LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz, lusol->luparm, lusol->parmlu, lusol->data, lusol->indc, lusol->indr, lusol->ip, lusol->iq, lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status);
208: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "solve failed, error code %d", status);
210: PetscCall(VecRestoreArray(x, &xx));
211: PetscCall(VecRestoreArrayRead(b, &bb));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: static PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F, Mat A, const MatFactorInfo *info)
216: {
217: Mat_SeqAIJ *a;
218: Mat_LUSOL *lusol = (Mat_LUSOL *)F->spptr;
219: int m, n, nz, nnz, status;
220: int i, rs, re;
221: int factorizations;
223: PetscFunctionBegin;
224: PetscCall(MatGetSize(A, &m, &n));
225: a = (Mat_SeqAIJ *)A->data;
227: PetscCheck(m == lusol->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "factorization struct inconsistent");
229: factorizations = 0;
230: do {
231: /*******************************************************************/
232: /* Check the workspace allocation. */
233: /*******************************************************************/
235: nz = a->nz;
236: nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom * nz));
237: nnz = PetscMax(nnz, 5 * n);
239: if (nnz < lusol->luparm[12]) {
240: nnz = (int)(lusol->luroom * lusol->luparm[12]);
241: } else if ((factorizations > 0) && (lusol->luroom < 6)) {
242: lusol->luroom += 0.1;
243: }
245: nnz = PetscMax(nnz, (int)(lusol->luroom * (lusol->luparm[22] + lusol->luparm[23])));
247: if (nnz > lusol->nnz) {
248: PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
249: PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));
250: lusol->nnz = nnz;
251: }
253: /* Fill in the data for the problem. (1-based Fortran style) */
254: nz = 0;
255: for (i = 0; i < n; i++) {
256: rs = a->i[i];
257: re = a->i[i + 1];
259: while (rs < re) {
260: if (a->a[rs] != 0.0) {
261: lusol->indc[nz] = i + 1;
262: lusol->indr[nz] = a->j[rs] + 1;
263: lusol->data[nz] = a->a[rs];
264: nz++;
265: }
266: rs++;
267: }
268: }
270: /* Do the factorization. */
271: LU1FAC(&m, &n, &nz, &nnz, lusol->luparm, lusol->parmlu, lusol->data, lusol->indc, lusol->indr, lusol->ip, lusol->iq, lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, lusol->iploc, lusol->iqloc, lusol->ipinv, lusol->iqinv, lusol->mnsw, &status);
273: switch (status) {
274: case 0: /* factored */
275: break;
277: case 7: /* insufficient memory */
278: break;
280: case 1:
281: case -1: /* singular */
282: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Singular matrix");
284: case 3:
285: case 4: /* error conditions */
286: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix error");
288: default: /* unknown condition */
289: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix unknown return code");
290: }
292: factorizations++;
293: } while (status == 7);
294: F->ops->solve = MatSolve_LUSOL;
295: F->assembled = PETSC_TRUE;
296: F->preallocated = PETSC_TRUE;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
301: {
302: /*
303: Input
304: A - matrix to factor
305: r - row permutation (ignored)
306: c - column permutation (ignored)
308: Output
309: F - matrix storing the factorization;
310: */
311: Mat_LUSOL *lusol;
312: int i, m, n, nz, nnz;
314: PetscFunctionBegin;
315: /* Check the arguments. */
316: PetscCall(MatGetSize(A, &m, &n));
317: nz = ((Mat_SeqAIJ *)A->data)->nz;
319: /* Create the factorization. */
320: F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
321: lusol = (Mat_LUSOL *)F->spptr;
323: /* Initialize parameters */
324: for (i = 0; i < 30; i++) {
325: lusol->luparm[i] = 0;
326: lusol->parmlu[i] = 0;
327: }
329: lusol->luparm[1] = -1;
330: lusol->luparm[2] = 5;
331: lusol->luparm[7] = 1;
333: lusol->parmlu[0] = 1 / Factorization_Tolerance;
334: lusol->parmlu[1] = 1 / Factorization_Tolerance;
335: lusol->parmlu[2] = Factorization_Small_Tolerance;
336: lusol->parmlu[3] = Factorization_Pivot_Tolerance;
337: lusol->parmlu[4] = Factorization_Pivot_Tolerance;
338: lusol->parmlu[5] = 3.0;
339: lusol->parmlu[6] = 0.3;
340: lusol->parmlu[7] = 0.6;
342: /* Allocate the workspace needed by LUSOL. */
343: lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
344: nnz = PetscMax((int)(lusol->elbowroom * nz), 5 * n);
346: lusol->n = n;
347: lusol->nz = nz;
348: lusol->nnz = nnz;
349: lusol->luroom = 1.75;
351: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ip));
352: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iq));
353: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenc));
354: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenr));
355: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locc));
356: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locr));
357: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iploc));
358: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqloc));
359: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ipinv));
360: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqinv));
361: PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsw));
362: PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsv));
363: PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));
365: lusol->CleanUpLUSOL = PETSC_TRUE;
366: F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A, MatSolverType *type)
371: {
372: PetscFunctionBegin;
373: *type = MATSOLVERLUSOL;
374: PetscFunctionReturn(PETSC_SUCCESS);
375: }
377: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A, MatFactorType ftype, Mat *F)
378: {
379: Mat B;
380: Mat_LUSOL *lusol;
381: int m, n;
383: PetscFunctionBegin;
384: PetscCall(MatGetSize(A, &m, &n));
385: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
386: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
387: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
388: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
390: PetscCall(PetscNew(&lusol));
391: B->spptr = lusol;
393: B->trivialsymbolic = PETSC_TRUE;
394: B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL;
395: B->ops->destroy = MatDestroy_LUSOL;
397: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_lusol));
399: B->factortype = MAT_FACTOR_LU;
400: PetscCall(PetscFree(B->solvertype));
401: PetscCall(PetscStrallocpy(MATSOLVERLUSOL, &B->solvertype));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Lusol(void)
406: {
407: PetscFunctionBegin;
408: PetscCall(MatSolverTypeRegister(MATSOLVERLUSOL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_lusol));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*MC
413: MATSOLVERLUSOL - "lusol" - Provides direct solvers, LU, for sequential matrices
414: via the external package LUSOL.
416: Works with `MATSEQAIJ` matrices
418: Level: beginner
420: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
421: M*/