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: /*
     76:     LUSOL input/Output Parameters (Description uses C-style indexes

     78:     Input parameters                                        Typical value
     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.

     95:     Output parameters:

     97:     luparm(9) = inform   Return code from last call to any LU routine.
     98:     luparm(10) = nsing    No. of singularities marked in the
     99:                           output array w(*).
    100:     luparm(11) = jsing    Column index of last singularity.
    101:     luparm(12) = minlen   Minimum recommended value for  lena.
    102:     luparm(13) = maxlen   ?
    103:     luparm(14) = nupdat   No. of updates performed by the lu8 routines.
    104:     luparm(15) = nrank    No. of nonempty rows of U.
    105:     luparm(16) = ndens1   No. of columns remaining when the density of
    106:                           the matrix being factorized reached dens1.
    107:     luparm(17) = ndens2   No. of columns remaining when the density of
    108:                           the matrix being factorized reached dens2.
    109:     luparm(18) = jumin    The column index associated with dumin.
    110:     luparm(19) = numl0    No. of columns in initial  L.
    111:     luparm(20) = lenl0    Size of initial  L  (no. of nonzeros).
    112:     luparm(21) = lenu0    Size of initial  U.
    113:     luparm(22) = lenl     Size of current  L.
    114:     luparm(23) = lenu     Size of current  U.
    115:     luparm(24) = lrow     Length of row file.
    116:     luparm(25) = ncp      No. of compressions of LU data structures.
    117:     luparm(26) = mersum   lu1fac: sum of Markowitz merit counts.
    118:     luparm(27) = nutri    lu1fac: triangular rows in U.
    119:     luparm(28) = nltri    lu1fac: triangular rows in L.
    120:     luparm(29) =

    122:     Input parameters                                        Typical value
    123:     parmlu(0) = elmax1   Max multiplier allowed in  L           10.0
    124:                           during factor.
    125:     parmlu(1) = elmax2   Max multiplier allowed in  L           10.0
    126:                           during updates.
    127:     parmlu(2) = small    Absolute tolerance for             eps**0.8
    128:                           treating reals as zero.     IBM double: 3.0d-13
    129:     parmlu(3) = utol1    Absolute tol for flagging          eps**0.66667
    130:                           small diagonals of U.       IBM double: 3.7d-11
    131:     parmlu(4) = utol2    Relative tol for flagging          eps**0.66667
    132:                           small diagonals of U.       IBM double: 3.7d-11
    133:     parmlu(5) = uspace   Factor limiting waste space in  U.      3.0
    134:                           In lu1fac, the row or column lists
    135:                           are compressed if their length
    136:                           exceeds uspace times the length of
    137:                           either file after the last compression.
    138:     parmlu(6) = dens1    The density at which the Markowitz      0.3
    139:                           strategy should search maxcol columns
    140:                           and no rows.
    141:     parmlu(7) = dens2    the density at which the Markowitz      0.6
    142:                           strategy should search only 1 column
    143:                           or (preferably) use a dense LU for
    144:                           all the remaining rows and columns.

    146:     Output parameters:
    147:     parmlu(9) = amax     Maximum element in  A.
    148:     parmlu(10) = elmax    Maximum multiplier in current  L.
    149:     parmlu(11) = umax     Maximum element in current  U.
    150:     parmlu(12) = dumax    Maximum diagonal in  U.
    151:     parmlu(13) = dumin    Minimum diagonal in  U.
    152:     parmlu(14) =
    153:     parmlu(15) =
    154:     parmlu(16) =
    155:     parmlu(17) =
    156:     parmlu(18) =
    157:     parmlu(19) = resid    lu6sol: residual after solve with U or U'.
    158:     ...
    159:     parmlu(29) =
    160:   */

    162: #define Factorization_Tolerance       1e-1
    163: #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0)
    164: #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */

    166: static PetscErrorCode MatDestroy_LUSOL(Mat A)
    167: {
    168:   Mat_LUSOL *lusol = (Mat_LUSOL *)A->spptr;

    170:   PetscFunctionBegin;
    171:   if (lusol && lusol->CleanUpLUSOL) {
    172:     PetscCall(PetscFree(lusol->ip));
    173:     PetscCall(PetscFree(lusol->iq));
    174:     PetscCall(PetscFree(lusol->lenc));
    175:     PetscCall(PetscFree(lusol->lenr));
    176:     PetscCall(PetscFree(lusol->locc));
    177:     PetscCall(PetscFree(lusol->locr));
    178:     PetscCall(PetscFree(lusol->iploc));
    179:     PetscCall(PetscFree(lusol->iqloc));
    180:     PetscCall(PetscFree(lusol->ipinv));
    181:     PetscCall(PetscFree(lusol->iqinv));
    182:     PetscCall(PetscFree(lusol->mnsw));
    183:     PetscCall(PetscFree(lusol->mnsv));
    184:     PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
    185:   }
    186:   PetscCall(PetscFree(A->spptr));
    187:   PetscCall(MatDestroy_SeqAIJ(A));
    188:   PetscFunctionReturn(PETSC_SUCCESS);
    189: }

    191: static PetscErrorCode MatSolve_LUSOL(Mat A, Vec b, Vec x)
    192: {
    193:   Mat_LUSOL    *lusol = (Mat_LUSOL *)A->spptr;
    194:   double       *xx;
    195:   const double *bb;
    196:   int           mode = 5;
    197:   int           i, m, n, nnz, status;

    199:   PetscFunctionBegin;
    200:   PetscCall(VecGetArray(x, &xx));
    201:   PetscCall(VecGetArrayRead(b, &bb));

    203:   m = n = lusol->n;
    204:   nnz   = lusol->nnz;

    206:   for (i = 0; i < m; i++) lusol->mnsv[i] = bb[i];

    208:   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);

    210:   PetscCheck(!status, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "solve failed, error code %d", status);

    212:   PetscCall(VecRestoreArray(x, &xx));
    213:   PetscCall(VecRestoreArrayRead(b, &bb));
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: static PetscErrorCode MatLUFactorNumeric_LUSOL(Mat F, Mat A, const MatFactorInfo *info)
    218: {
    219:   Mat_SeqAIJ *a;
    220:   Mat_LUSOL  *lusol = (Mat_LUSOL *)F->spptr;
    221:   int         m, n, nz, nnz, status;
    222:   int         i, rs, re;
    223:   int         factorizations;

    225:   PetscFunctionBegin;
    226:   PetscCall(MatGetSize(A, &m, &n));
    227:   a = (Mat_SeqAIJ *)A->data;

    229:   PetscCheck(m == lusol->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "factorization struct inconsistent");

    231:   factorizations = 0;
    232:   do {
    233:     /*******************************************************************/
    234:     /* Check the workspace allocation.                                 */
    235:     /*******************************************************************/

    237:     nz  = a->nz;
    238:     nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom * nz));
    239:     nnz = PetscMax(nnz, 5 * n);

    241:     if (nnz < lusol->luparm[12]) {
    242:       nnz = (int)(lusol->luroom * lusol->luparm[12]);
    243:     } else if ((factorizations > 0) && (lusol->luroom < 6)) {
    244:       lusol->luroom += 0.1;
    245:     }

    247:     nnz = PetscMax(nnz, (int)(lusol->luroom * (lusol->luparm[22] + lusol->luparm[23])));

    249:     if (nnz > lusol->nnz) {
    250:       PetscCall(PetscFree3(lusol->data, lusol->indc, lusol->indr));
    251:       PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));
    252:       lusol->nnz = nnz;
    253:     }

    255:     /* Fill in the data for the problem.      (1-based Fortran style)  */
    256:     nz = 0;
    257:     for (i = 0; i < n; i++) {
    258:       rs = a->i[i];
    259:       re = a->i[i + 1];

    261:       while (rs < re) {
    262:         if (a->a[rs] != 0.0) {
    263:           lusol->indc[nz] = i + 1;
    264:           lusol->indr[nz] = a->j[rs] + 1;
    265:           lusol->data[nz] = a->a[rs];
    266:           nz++;
    267:         }
    268:         rs++;
    269:       }
    270:     }

    272:     /* Do the factorization.                                           */
    273:     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);

    275:     switch (status) {
    276:     case 0: /* factored */
    277:       break;

    279:     case 7: /* insufficient memory */
    280:       break;

    282:     case 1:
    283:     case -1: /* singular */
    284:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Singular matrix");

    286:     case 3:
    287:     case 4: /* error conditions */
    288:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix error");

    290:     default: /* unknown condition */
    291:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "matrix unknown return code");
    292:     }

    294:     factorizations++;
    295:   } while (status == 7);
    296:   F->ops->solve   = MatSolve_LUSOL;
    297:   F->assembled    = PETSC_TRUE;
    298:   F->preallocated = PETSC_TRUE;
    299:   PetscFunctionReturn(PETSC_SUCCESS);
    300: }

    302: static PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
    303: {
    304:   /*
    305:      Input
    306:          A  - matrix to factor
    307:          r  - row permutation (ignored)
    308:          c  - column permutation (ignored)

    310:      Output
    311:          F  - matrix storing the factorization;
    312:   */
    313:   Mat_LUSOL *lusol;
    314:   int        i, m, n, nz, nnz;

    316:   PetscFunctionBegin;
    317:   /* Check the arguments.                                                 */
    318:   PetscCall(MatGetSize(A, &m, &n));
    319:   nz = ((Mat_SeqAIJ *)A->data)->nz;

    321:   /* Create the factorization.                                            */
    322:   F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
    323:   lusol                   = (Mat_LUSOL *)F->spptr;

    325:   /* Initialize parameters                                                */
    326:   for (i = 0; i < 30; i++) {
    327:     lusol->luparm[i] = 0;
    328:     lusol->parmlu[i] = 0;
    329:   }

    331:   lusol->luparm[1] = -1;
    332:   lusol->luparm[2] = 5;
    333:   lusol->luparm[7] = 1;

    335:   lusol->parmlu[0] = 1 / Factorization_Tolerance;
    336:   lusol->parmlu[1] = 1 / Factorization_Tolerance;
    337:   lusol->parmlu[2] = Factorization_Small_Tolerance;
    338:   lusol->parmlu[3] = Factorization_Pivot_Tolerance;
    339:   lusol->parmlu[4] = Factorization_Pivot_Tolerance;
    340:   lusol->parmlu[5] = 3.0;
    341:   lusol->parmlu[6] = 0.3;
    342:   lusol->parmlu[7] = 0.6;

    344:   /* Allocate the workspace needed by LUSOL.                              */
    345:   lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
    346:   nnz              = PetscMax((int)(lusol->elbowroom * nz), 5 * n);

    348:   lusol->n      = n;
    349:   lusol->nz     = nz;
    350:   lusol->nnz    = nnz;
    351:   lusol->luroom = 1.75;

    353:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ip));
    354:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iq));
    355:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenc));
    356:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->lenr));
    357:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locc));
    358:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->locr));
    359:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iploc));
    360:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqloc));
    361:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->ipinv));
    362:   PetscCall(PetscMalloc(sizeof(int) * n, &lusol->iqinv));
    363:   PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsw));
    364:   PetscCall(PetscMalloc(sizeof(double) * n, &lusol->mnsv));
    365:   PetscCall(PetscMalloc3(nnz, &lusol->data, nnz, &lusol->indc, nnz, &lusol->indr));

    367:   lusol->CleanUpLUSOL     = PETSC_TRUE;
    368:   F->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
    369:   PetscFunctionReturn(PETSC_SUCCESS);
    370: }

    372: static PetscErrorCode MatFactorGetSolverType_seqaij_lusol(Mat A, MatSolverType *type)
    373: {
    374:   PetscFunctionBegin;
    375:   *type = MATSOLVERLUSOL;
    376:   PetscFunctionReturn(PETSC_SUCCESS);
    377: }

    379: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_lusol(Mat A, MatFactorType ftype, Mat *F)
    380: {
    381:   Mat        B;
    382:   Mat_LUSOL *lusol;
    383:   int        m, n;

    385:   PetscFunctionBegin;
    386:   PetscCall(MatGetSize(A, &m, &n));
    387:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    388:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, m, n));
    389:   PetscCall(MatSetType(B, ((PetscObject)A)->type_name));
    390:   PetscCall(MatSeqAIJSetPreallocation(B, 0, NULL));

    392:   PetscCall(PetscNew(&lusol));
    393:   B->spptr = lusol;

    395:   B->trivialsymbolic       = PETSC_TRUE;
    396:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_LUSOL;
    397:   B->ops->destroy          = MatDestroy_LUSOL;

    399:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_seqaij_lusol));

    401:   B->factortype = MAT_FACTOR_LU;
    402:   PetscCall(PetscFree(B->solvertype));
    403:   PetscCall(PetscStrallocpy(MATSOLVERLUSOL, &B->solvertype));
    404:   PetscFunctionReturn(PETSC_SUCCESS);
    405: }

    407: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Lusol(void)
    408: {
    409:   PetscFunctionBegin;
    410:   PetscCall(MatSolverTypeRegister(MATSOLVERLUSOL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_lusol));
    411:   PetscFunctionReturn(PETSC_SUCCESS);
    412: }

    414: /*MC
    415:   MATSOLVERLUSOL - "lusol" - Provides direct solvers, LU, for sequential matrices
    416:                    via the external package LUSOL.

    418:   Works with `MATSEQAIJ` matrices

    420:   Level: beginner

    422: .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
    423: M*/