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: extern PetscErrorCode MatDuplicate_LUSOL(Mat, MatDuplicateOption, Mat *);
44: typedef struct {
45: double *data;
46: int *indc;
47: int *indr;
49: int *ip;
50: int *iq;
51: int *lenc;
52: int *lenr;
53: int *locc;
54: int *locr;
55: int *iploc;
56: int *iqloc;
57: int *ipinv;
58: int *iqinv;
59: double *mnsw;
60: double *mnsv;
62: double elbowroom;
63: double luroom; /* Extra space allocated when factor fails */
64: double parmlu[30]; /* Input/output to LUSOL */
66: int n; /* Number of rows/columns in matrix */
67: int nz; /* Number of nonzeros */
68: int nnz; /* Number of nonzeros allocated for factors */
69: int luparm[30]; /* Input/output to LUSOL */
71: PetscBool CleanUpLUSOL;
73: } Mat_LUSOL;
75: /* LUSOL input/Output Parameters (Description uses C-style indexes
76: *
77: * Input parameters Typical value
78: *
79: * luparm(0) = nout File number for printed messages. 6
80: * luparm(1) = lprint Print level. 0
81: * < 0 suppresses output.
82: * = 0 gives error messages.
83: * = 1 gives debug output from some of the
84: * other routines in LUSOL.
85: * >= 2 gives the pivot row and column and the
86: * no. of rows and columns involved at
87: * each elimination step in lu1fac.
88: * luparm(2) = maxcol lu1fac: maximum number of columns 5
89: * searched allowed in a Markowitz-type
90: * search for the next pivot element.
91: * For some of the factorization, the
92: * number of rows searched is
93: * maxrow = maxcol - 1.
94: *
95: *
96: * Output parameters:
97: *
98: * luparm(9) = inform Return code from last call to any LU routine.
99: * luparm(10) = nsing No. of singularities marked in the
100: * output array w(*).
101: * luparm(11) = jsing Column index of last singularity.
102: * luparm(12) = minlen Minimum recommended value for lena.
103: * luparm(13) = maxlen ?
104: * luparm(14) = nupdat No. of updates performed by the lu8 routines.
105: * luparm(15) = nrank No. of nonempty rows of U.
106: * luparm(16) = ndens1 No. of columns remaining when the density of
107: * the matrix being factorized reached dens1.
108: * luparm(17) = ndens2 No. of columns remaining when the density of
109: * the matrix being factorized reached dens2.
110: * luparm(18) = jumin The column index associated with dumin.
111: * luparm(19) = numl0 No. of columns in initial L.
112: * luparm(20) = lenl0 Size of initial L (no. of nonzeros).
113: * luparm(21) = lenu0 Size of initial U.
114: * luparm(22) = lenl Size of current L.
115: * luparm(23) = lenu Size of current U.
116: * luparm(24) = lrow Length of row file.
117: * luparm(25) = ncp No. of compressions of LU data structures.
118: * luparm(26) = mersum lu1fac: sum of Markowitz merit counts.
119: * luparm(27) = nutri lu1fac: triangular rows in U.
120: * luparm(28) = nltri lu1fac: triangular rows in L.
121: * luparm(29) =
122: *
123: *
124: * Input parameters Typical value
125: *
126: * parmlu(0) = elmax1 Max multiplier allowed in L 10.0
127: * during factor.
128: * parmlu(1) = elmax2 Max multiplier allowed in L 10.0
129: * during updates.
130: * parmlu(2) = small Absolute tolerance for eps**0.8
131: * treating reals as zero. IBM double: 3.0d-13
132: * parmlu(3) = utol1 Absolute tol for flagging eps**0.66667
133: * small diagonals of U. IBM double: 3.7d-11
134: * parmlu(4) = utol2 Relative tol for flagging eps**0.66667
135: * small diagonals of U. IBM double: 3.7d-11
136: * parmlu(5) = uspace Factor limiting waste space in U. 3.0
137: * In lu1fac, the row or column lists
138: * are compressed if their length
139: * exceeds uspace times the length of
140: * either file after the last compression.
141: * parmlu(6) = dens1 The density at which the Markowitz 0.3
142: * strategy should search maxcol columns
143: * and no rows.
144: * parmlu(7) = dens2 the density at which the Markowitz 0.6
145: * strategy should search only 1 column
146: * or (preferably) use a dense LU for
147: * all the remaining rows and columns.
148: *
149: *
150: * Output parameters:
151: *
152: * parmlu(9) = amax Maximum element in A.
153: * parmlu(10) = elmax Maximum multiplier in current L.
154: * parmlu(11) = umax Maximum element in current U.
155: * parmlu(12) = dumax Maximum diagonal in U.
156: * parmlu(13) = dumin Minimum diagonal in U.
157: * parmlu(14) =
158: * parmlu(15) =
159: * parmlu(16) =
160: * parmlu(17) =
161: * parmlu(18) =
162: * parmlu(19) = resid lu6sol: residual after solve with U or U'.
163: * ...
164: * parmlu(29) =
165: */
167: #define Factorization_Tolerance 1e-1
168: #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0)
169: #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */
171: static PetscErrorCode MatDestroy_LUSOL(Mat A)
172: {
173: Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;
175: PetscFunctionBegin;
176: if (lusol && lusol->CleanUpLUSOL) {
177: PetscCall(PetscFree(lusol->ip));
178: PetscCall(PetscFree(lusol->iq));
179: PetscCall(PetscFree(lusol->lenc));
180: PetscCall(PetscFree(lusol->lenr));
181: PetscCall(PetscFree(lusol->locc));
182: PetscCall(PetscFree(lusol->locr));
183: PetscCall(PetscFree(lusol->iploc));
184: PetscCall(PetscFree(lusol->iqloc));
185: PetscCall(PetscFree(lusol->ipinv));
186: PetscCall(PetscFree(lusol->iqinv));
187: PetscCall(PetscFree(lusol->mnsw));
188: PetscCall(PetscFree(lusol->mnsv));
189: PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
190: }
191: PetscCall(PetscFree(A->spptr));
192: PetscCall(MatDestroy_SeqAIJ(A));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode MatSolve_LUSOL(Mat A, Vec b, Vec x)
197: {
198: Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;
199: double *xx;
200: const double *bb;
201: int mode = 5;
202: int i, m, n, nnz, status;
204: PetscFunctionBegin;
205: PetscCall(VecGetArray(x, &xx));
206: PetscCall(VecGetArrayRead(b, &bb));
208: m = n = lusol->n;
209: nnz = lusol->nnz;
211: for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i];
213: 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);
215: PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "solve failed, error code %d", status);
217: PetscCall(VecRestoreArray(x, &xx));
218: PetscCall(VecRestoreArrayRead(b, &bb));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F, Mat A, const MatFactorInfo *info)
223: {
224: Mat_SeqAIJ *a;
225: Mat_LUSOL *lusol = (Mat_LUSOL *)F->spptr;
226: int m, n, nz, nnz, status;
227: int i, rs, re;
228: int factorizations;
230: PetscFunctionBegin;
231: PetscCall(MatGetSize(A, &m, &n));
232: a = (Mat_SeqAIJ *)A->data;
234: PetscCheck(m == lusol->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "factorization struct inconsistent");
236: factorizations = 0;
237: do {
238: /*******************************************************************/
239: /* Check the workspace allocation. */
240: /*******************************************************************/
242: nz = a->nz;
243: nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom * nz));
244: nnz = PetscMax(nnz, 5 * n);
246: if (nnz < lusol->luparm[12]) {
247: nnz = (int)(lusol->luroom * lusol->luparm[12]);
248: } else if ((factorizations > 0) && (lusol->luroom < 6)) {
249: lusol->luroom += 0.1;
250: }
252: nnz = PetscMax(nnz, (int)(lusol->luroom * (lusol->luparm[22] + lusol->luparm[23])));
254: if (nnz > lusol->nnz) {
255: PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
256: PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));
257: lusol->nnz = nnz;
258: }
260: /*******************************************************************/
261: /* Fill in the data for the problem. (1-based Fortran style) */
262: /*******************************************************************/
264: nz = 0;
265: for (i = 0; i < n; i++) {
266: rs = a->i[i];
267: re = a->i[i + 1];
269: while (rs < re) {
270: if (a->a[rs] != 0.0) {
271: lusol->indc[nz] = i + 1;
272: lusol->indr[nz] = a->j[rs] + 1;
273: lusol->data[nz] = a->a[rs];
274: nz++;
275: }
276: rs++;
277: }
278: }
280: /*******************************************************************/
281: /* Do the factorization. */
282: /*******************************************************************/
284: 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);
286: switch (status) {
287: case 0: /* factored */
288: break;
290: case 7: /* insufficient memory */
291: break;
293: case 1:
294: case -1: /* singular */
295: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Singular matrix");
297: case 3:
298: case 4: /* error conditions */
299: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix error");
301: default: /* unknown condition */
302: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix unknown return code");
303: }
305: factorizations++;
306: } while (status == 7);
307: F->ops->solve = MatSolve_LUSOL;
308: F->assembled = PETSC_TRUE;
309: F->preallocated = PETSC_TRUE;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
314: {
315: /************************************************************************/
316: /* Input */
317: /* A - matrix to factor */
318: /* r - row permutation (ignored) */
319: /* c - column permutation (ignored) */
320: /* */
321: /* Output */
322: /* F - matrix storing the factorization; */
323: /************************************************************************/
324: Mat_LUSOL *lusol;
325: int i, m, n, nz, nnz;
327: PetscFunctionBegin;
328: /************************************************************************/
329: /* Check the arguments. */
330: /************************************************************************/
332: PetscCall(MatGetSize(A, &m, &n));
333: nz = ((Mat_SeqAIJ *)A->data)->nz;
335: /************************************************************************/
336: /* Create the factorization. */
337: /************************************************************************/
339: F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
340: lusol = (Mat_LUSOL *)(F->spptr);
342: /************************************************************************/
343: /* Initialize parameters */
344: /************************************************************************/
346: for (i = 0; i < 30; i++) {
347: lusol->luparm[i] = 0;
348: lusol->parmlu[i] = 0;
349: }
351: lusol->luparm[1] = -1;
352: lusol->luparm[2] = 5;
353: lusol->luparm[7] = 1;
355: lusol->parmlu[0] = 1 / Factorization_Tolerance;
356: lusol->parmlu[1] = 1 / Factorization_Tolerance;
357: lusol->parmlu[2] = Factorization_Small_Tolerance;
358: lusol->parmlu[3] = Factorization_Pivot_Tolerance;
359: lusol->parmlu[4] = Factorization_Pivot_Tolerance;
360: lusol->parmlu[5] = 3.0;
361: lusol->parmlu[6] = 0.3;
362: lusol->parmlu[7] = 0.6;
364: /************************************************************************/
365: /* Allocate the workspace needed by LUSOL. */
366: /************************************************************************/
368: lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
369: nnz = PetscMax((int)(lusol->elbowroom * nz), 5 * n);
371: lusol->n = n;
372: lusol->nz = nz;
373: lusol->nnz = nnz;
374: lusol->luroom = 1.75;
376: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ip));
377: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iq));
378: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenc));
379: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenr));
380: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locc));
381: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locr));
382: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iploc));
383: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqloc));
384: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ipinv));
385: PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqinv));
386: PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsw));
387: PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsv));
388: PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));
390: lusol->CleanUpLUSOL = PETSC_TRUE;
391: F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: static PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A, MatSolverType *type)
396: {
397: PetscFunctionBegin;
398: *type = MATSOLVERLUSOL;
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A, MatFactorType ftype, Mat *F)
403: {
404: Mat B;
405: Mat_LUSOL *lusol;
406: int m, n;
408: PetscFunctionBegin;
409: PetscCall(MatGetSize(A, &m, &n));
410: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
411: PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
412: PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
413: PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));
415: PetscCall(PetscNew(&lusol));
416: B->spptr = lusol;
418: B->trivialsymbolic = PETSC_TRUE;
419: B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL;
420: B->ops->destroy = MatDestroy_LUSOL;
422: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_lusol));
424: B->factortype = MAT_FACTOR_LU;
425: PetscCall(PetscFree(B->solvertype));
426: PetscCall(PetscStrallocpy(MATSOLVERLUSOL, &B->solvertype));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Lusol(void)
432: {
433: PetscFunctionBegin;
434: PetscCall(MatSolverTypeRegister(MATSOLVERLUSOL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_lusol));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*MC
439: MATSOLVERLUSOL - "lusol" - Provides direct solvers, LU, for sequential matrices
440: via the external package LUSOL.
442: Works with `MATSEQAIJ` matrices
444: Level: beginner
446: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
447: M*/