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