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