1 /*
2 Provides an interface to the LUSOL package of ....
3
4 */
5 #include <../src/mat/impls/aij/seq/aij.h>
6
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
20
21 /*
22 Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require
23 */
M1PAGE()24 PETSC_EXTERN void M1PAGE()
25 {
26 ;
27 }
M5SETX()28 PETSC_EXTERN void M5SETX()
29 {
30 ;
31 }
32
M6RDEL()33 PETSC_EXTERN void M6RDEL()
34 {
35 ;
36 }
37
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);
39
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);
41
42 typedef struct {
43 double *data;
44 int *indc;
45 int *indr;
46
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;
59
60 double elbowroom;
61 double luroom; /* Extra space allocated when factor fails */
62 double parmlu[30]; /* Input/output to LUSOL */
63
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 */
68
69 PetscBool CleanUpLUSOL;
70
71 } Mat_LUSOL;
72
73 /*
74 LUSOL input/Output Parameters (Description uses C-style indexes
75
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.
92
93 Output parameters:
94
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) =
119
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.
143
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 */
159
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) */
163
MatDestroy_LUSOL(Mat A)164 static PetscErrorCode MatDestroy_LUSOL(Mat A)
165 {
166 Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;
167
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 }
188
MatSolve_LUSOL(Mat A,Vec b,Vec x)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;
196
197 PetscFunctionBegin;
198 PetscCall(VecGetArray(x, &xx));
199 PetscCall(VecGetArrayRead(b, &bb));
200
201 m = n = lusol->n;
202 nnz = lusol->nnz;
203
204 for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i];
205
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);
207
208 PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "solve failed, error code %d", status);
209
210 PetscCall(VecRestoreArray(x, &xx));
211 PetscCall(VecRestoreArrayRead(b, &bb));
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
MatLUFactorNumeric_LUSOL(Mat F,Mat A,const MatFactorInfo * info)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;
222
223 PetscFunctionBegin;
224 PetscCall(MatGetSize(A, &m, &n));
225 a = (Mat_SeqAIJ *)A->data;
226
227 PetscCheck(m == lusol->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "factorization struct inconsistent");
228
229 factorizations = 0;
230 do {
231 /*******************************************************************/
232 /* Check the workspace allocation. */
233 /*******************************************************************/
234
235 nz = a->nz;
236 nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom * nz));
237 nnz = PetscMax(nnz, 5 * n);
238
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 }
244
245 nnz = PetscMax(nnz, (int)(lusol->luroom * (lusol->luparm[22] + lusol->luparm[23])));
246
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 }
252
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];
258
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 }
269
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);
272
273 switch (status) {
274 case 0: /* factored */
275 break;
276
277 case 7: /* insufficient memory */
278 break;
279
280 case 1:
281 case -1: /* singular */
282 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Singular matrix");
283
284 case 3:
285 case 4: /* error conditions */
286 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix error");
287
288 default: /* unknown condition */
289 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix unknown return code");
290 }
291
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 }
299
MatLUFactorSymbolic_LUSOL(Mat F,Mat A,IS r,IS c,const MatFactorInfo * info)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)
307
308 Output
309 F - matrix storing the factorization;
310 */
311 Mat_LUSOL *lusol;
312 int i, m, n, nz, nnz;
313
314 PetscFunctionBegin;
315 /* Check the arguments. */
316 PetscCall(MatGetSize(A, &m, &n));
317 nz = ((Mat_SeqAIJ *)A->data)->nz;
318
319 /* Create the factorization. */
320 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
321 lusol = (Mat_LUSOL *)F->spptr;
322
323 /* Initialize parameters */
324 for (i = 0; i < 30; i++) {
325 lusol->luparm[i] = 0;
326 lusol->parmlu[i] = 0;
327 }
328
329 lusol->luparm[1] = -1;
330 lusol->luparm[2] = 5;
331 lusol->luparm[7] = 1;
332
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;
341
342 /* Allocate the workspace needed by LUSOL. */
343 lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
344 nnz = PetscMax((int)(lusol->elbowroom * nz), 5 * n);
345
346 lusol->n = n;
347 lusol->nz = nz;
348 lusol->nnz = nnz;
349 lusol->luroom = 1.75;
350
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));
364
365 lusol->CleanUpLUSOL = PETSC_TRUE;
366 F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
MatFactorGetSolverType_seqaij_lusol(Mat A,MatSolverType * type)370 static PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A, MatSolverType *type)
371 {
372 PetscFunctionBegin;
373 *type = MATSOLVERLUSOL;
374 PetscFunctionReturn(PETSC_SUCCESS);
375 }
376
MatGetFactor_seqaij_lusol(Mat A,MatFactorType ftype,Mat * F)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;
382
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));
389
390 PetscCall(PetscNew(&lusol));
391 B->spptr = lusol;
392
393 B->trivialsymbolic = PETSC_TRUE;
394 B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL;
395 B->ops->destroy = MatDestroy_LUSOL;
396
397 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_lusol));
398
399 B->factortype = MAT_FACTOR_LU;
400 PetscCall(PetscFree(B->solvertype));
401 PetscCall(PetscStrallocpy(MATSOLVERLUSOL, &B->solvertype));
402 PetscFunctionReturn(PETSC_SUCCESS);
403 }
404
MatSolverTypeRegister_Lusol(void)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 }
411
412 /*MC
413 MATSOLVERLUSOL - "lusol" - Provides direct solvers, LU, for sequential matrices
414 via the external package LUSOL.
415
416 Works with `MATSEQAIJ` matrices
417
418 Level: beginner
419
420 .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
421 M*/
422