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*/