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