Actual source code: lusol.c

      1: #define PETSCMAT_DLL

      3: /*
      4:         Provides an interface to the LUSOL package of ....

      6: */
     7:  #include src/mat/impls/aij/seq/aij.h

      9: #if defined(PETSC_HAVE_FORTRAN_UNDERSCORE)
     10: #define LU1FAC   lu1fac_
     11: #define LU6SOL   lu6sol_
     12: #define M1PAGE   m1page_
     13: #define M5SETX   m5setx_
     14: #define M6RDEL   m6rdel_
     15: #elif !defined(PETSC_HAVE_FORTRAN_CAPS)
     16: #define LU1FAC   lu1fac
     17: #define LU6SOL   lu6sol
     18: #define M1PAGE   m1page
     19: #define M5SETX   m5setx
     20: #define M6RDEL   m6rdel
     21: #endif

     24: /*
     25:     Dummy symbols that the MINOS files mi25bfac.f and mi15blas.f may require
     26: */
     27: void PETSC_STDCALL M1PAGE() {
     28:   ;
     29: }
     30: void PETSC_STDCALL M5SETX() {
     31:   ;
     32: }

     34: void PETSC_STDCALL M6RDEL() {
     35:   ;
     36: }

     39:                         double *parmlu, double *data, int *indc, int *indr,
     40:                         int *rowperm, int *colperm, int *collen, int *rowlen,
     41:                         int *colstart, int *rowstart, int *rploc, int *cploc,
     42:                         int *rpinv, int *cpinv, double *w, int *inform);

     45:                         int *size, int *luparm, double *parmlu, double *data,
     46:                         int *indc, int *indr, int *rowperm, int *colperm,
     47:                         int *collen, int *rowlen, int *colstart, int *rowstart,
     48:                         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:   PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
     81:   PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
     82:   PetscErrorCode (*MatDestroy)(Mat);
     83:   PetscTruth CleanUpLUSOL;

     85: } Mat_LUSOL;

     87: /*  LUSOL input/Output Parameters (Description uses C-style indexes
     88:  *
     89:  *  Input parameters                                        Typical value
     90:  *
     91:  *  luparm(0) = nout     File number for printed messages.         6
     92:  *  luparm(1) = lprint   Print level.                              0
     93:  *                    < 0 suppresses output.
     94:  *                    = 0 gives error messages.
     95:  *                    = 1 gives debug output from some of the
     96:  *                        other routines in LUSOL.
     97:  *                   >= 2 gives the pivot row and column and the
     98:  *                        no. of rows and columns involved at
     99:  *                        each elimination step in lu1fac.
    100:  *  luparm(2) = maxcol   lu1fac: maximum number of columns         5
    101:  *                        searched allowed in a Markowitz-type
    102:  *                        search for the next pivot element.
    103:  *                        For some of the factorization, the
    104:  *                        number of rows searched is
    105:  *                        maxrow = maxcol - 1.
    106:  *
    107:  *
    108:  *  Output parameters
    109:  *
    110:  *  luparm(9) = inform   Return code from last call to any LU routine.
    111:  *  luparm(10) = nsing    No. of singularities marked in the
    112:  *                        output array w(*).
    113:  *  luparm(11) = jsing    Column index of last singularity.
    114:  *  luparm(12) = minlen   Minimum recommended value for  lena.
    115:  *  luparm(13) = maxlen   ?
    116:  *  luparm(14) = nupdat   No. of updates performed by the lu8 routines.
    117:  *  luparm(15) = nrank    No. of nonempty rows of U.
    118:  *  luparm(16) = ndens1   No. of columns remaining when the density of
    119:  *                        the matrix being factorized reached dens1.
    120:  *  luparm(17) = ndens2   No. of columns remaining when the density of
    121:  *                        the matrix being factorized reached dens2.
    122:  *  luparm(18) = jumin    The column index associated with dumin.
    123:  *  luparm(19) = numl0    No. of columns in initial  L.
    124:  *  luparm(20) = lenl0    Size of initial  L  (no. of nonzeros).
    125:  *  luparm(21) = lenu0    Size of initial  U.
    126:  *  luparm(22) = lenl     Size of current  L.
    127:  *  luparm(23) = lenu     Size of current  U.
    128:  *  luparm(24) = lrow     Length of row file.
    129:  *  luparm(25) = ncp      No. of compressions of LU data structures.
    130:  *  luparm(26) = mersum   lu1fac: sum of Markowitz merit counts.
    131:  *  luparm(27) = nutri    lu1fac: triangular rows in U.
    132:  *  luparm(28) = nltri    lu1fac: triangular rows in L.
    133:  *  luparm(29) =
    134:  *
    135:  *
    136:  *  Input parameters                                        Typical value
    137:  *
    138:  *  parmlu(0) = elmax1   Max multiplier allowed in  L           10.0
    139:  *                        during factor.
    140:  *  parmlu(1) = elmax2   Max multiplier allowed in  L           10.0
    141:  *                        during updates.
    142:  *  parmlu(2) = small    Absolute tolerance for             eps**0.8
    143:  *                        treating reals as zero.     IBM double: 3.0d-13
    144:  *  parmlu(3) = utol1    Absolute tol for flagging          eps**0.66667
    145:  *                        small diagonals of U.       IBM double: 3.7d-11
    146:  *  parmlu(4) = utol2    Relative tol for flagging          eps**0.66667
    147:  *                        small diagonals of U.       IBM double: 3.7d-11
    148:  *  parmlu(5) = uspace   Factor limiting waste space in  U.      3.0
    149:  *                        In lu1fac, the row or column lists
    150:  *                        are compressed if their length
    151:  *                        exceeds uspace times the length of
    152:  *                        either file after the last compression.
    153:  *  parmlu(6) = dens1    The density at which the Markowitz      0.3
    154:  *                        strategy should search maxcol columns
    155:  *                        and no rows.
    156:  *  parmlu(7) = dens2    the density at which the Markowitz      0.6
    157:  *                        strategy should search only 1 column
    158:  *                        or (preferably) use a dense LU for
    159:  *                        all the remaining rows and columns.
    160:  *
    161:  *
    162:  *  Output parameters
    163:  *
    164:  *  parmlu(9) = amax     Maximum element in  A.
    165:  *  parmlu(10) = elmax    Maximum multiplier in current  L.
    166:  *  parmlu(11) = umax     Maximum element in current  U.
    167:  *  parmlu(12) = dumax    Maximum diagonal in  U.
    168:  *  parmlu(13) = dumin    Minimum diagonal in  U.
    169:  *  parmlu(14) =
    170:  *  parmlu(15) =
    171:  *  parmlu(16) =
    172:  *  parmlu(17) =
    173:  *  parmlu(18) =
    174:  *  parmlu(19) = resid    lu6sol: residual after solve with U or U'.
    175:  *  ...
    176:  *  parmlu(29) =
    177:  */

    179: #define Factorization_Tolerance       1e-1
    180: #define Factorization_Pivot_Tolerance pow(2.2204460492503131E-16, 2.0 / 3.0)
    181: #define Factorization_Small_Tolerance 1e-15 /* pow(DBL_EPSILON, 0.8) */

    186: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_LUSOL_SeqAIJ(Mat A,const MatType type,MatReuse reuse,Mat *newmat)
    187: {
    188:   /* This routine is only called to convert an unfactored PETSc-LUSOL matrix */
    189:   /* to its base PETSc type, so we will ignore 'MatType type'. */
    191:   Mat            B=*newmat;
    192:   Mat_LUSOL      *lusol=(Mat_LUSOL *)A->spptr;

    195:   if (reuse == MAT_INITIAL_MATRIX) {
    196:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    197:   }
    198:   B->ops->duplicate        = lusol->MatDuplicate;
    199:   B->ops->lufactorsymbolic = lusol->MatLUFactorSymbolic;
    200:   B->ops->destroy          = lusol->MatDestroy;
    201:
    202:   PetscFree(lusol);

    204:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_lusol_C","",PETSC_NULL);
    205:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_lusol_seqaij_C","",PETSC_NULL);

    207:   PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
    208:   *newmat = B;
    209:   return(0);
    210: }

    215: PetscErrorCode MatDestroy_LUSOL(Mat A)
    216: {
    218:   Mat_LUSOL *lusol=(Mat_LUSOL *)A->spptr;

    221:   if (lusol->CleanUpLUSOL) {
    222:     PetscFree(lusol->ip);
    223:     PetscFree(lusol->iq);
    224:     PetscFree(lusol->lenc);
    225:     PetscFree(lusol->lenr);
    226:     PetscFree(lusol->locc);
    227:     PetscFree(lusol->locr);
    228:     PetscFree(lusol->iploc);
    229:     PetscFree(lusol->iqloc);
    230:     PetscFree(lusol->ipinv);
    231:     PetscFree(lusol->iqinv);
    232:     PetscFree(lusol->mnsw);
    233:     PetscFree(lusol->mnsv);
    234:
    235:     PetscFree(lusol->indc);
    236:   }

    238:   MatConvert_LUSOL_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
    239:   (*A->ops->destroy)(A);
    240:   return(0);
    241: }

    245: PetscErrorCode MatSolve_LUSOL(Mat A,Vec b,Vec x)
    246: {
    247:   Mat_LUSOL *lusol=(Mat_LUSOL*)A->spptr;
    248:   double    *bb,*xx;
    249:   int       mode=5;
    251:   int       i,m,n,nnz,status;

    254:   VecGetArray(x, &xx);
    255:   VecGetArray(b, &bb);

    257:   m = n = lusol->n;
    258:   nnz = lusol->nnz;

    260:   for (i = 0; i < m; i++)
    261:     {
    262:       lusol->mnsv[i] = bb[i];
    263:     }

    265:   LU6SOL(&mode, &m, &n, lusol->mnsv, xx, &nnz,
    266:          lusol->luparm, lusol->parmlu, lusol->data,
    267:          lusol->indc, lusol->indr, lusol->ip, lusol->iq,
    268:          lusol->lenc, lusol->lenr, lusol->locc, lusol->locr, &status);

    270:   if (status != 0)
    271:     {
    272:       SETERRQ(PETSC_ERR_ARG_SIZ,"solve failed");
    273:     }

    275:   VecRestoreArray(x, &xx);
    276:   VecRestoreArray(b, &bb);
    277:   return(0);
    278: }

    282: PetscErrorCode MatLUFactorNumeric_LUSOL(Mat A,MatFactorInfo *info,Mat *F)
    283: {
    284:   Mat_SeqAIJ     *a;
    285:   Mat_LUSOL      *lusol = (Mat_LUSOL*)(*F)->spptr;
    287:   int            m, n, nz, nnz, status;
    288:   int            i, rs, re;
    289:   int            factorizations;

    292:   MatGetSize(A,&m,&n);
    293:   a = (Mat_SeqAIJ *)A->data;

    295:   if (m != lusol->n) {
    296:     SETERRQ(PETSC_ERR_ARG_SIZ,"factorization struct inconsistent");
    297:   }

    299:   factorizations = 0;
    300:   do
    301:     {
    302:       /*******************************************************************/
    303:       /* Check the workspace allocation.                                 */
    304:       /*******************************************************************/

    306:       nz = a->nz;
    307:       nnz = PetscMax(lusol->nnz, (int)(lusol->elbowroom*nz));
    308:       nnz = PetscMax(nnz, 5*n);

    310:       if (nnz < lusol->luparm[12]){
    311:         nnz = (int)(lusol->luroom * lusol->luparm[12]);
    312:       } else if ((factorizations > 0) && (lusol->luroom < 6)){
    313:         lusol->luroom += 0.1;
    314:       }

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

    318:       if (nnz > lusol->nnz){
    319:         PetscFree(lusol->indc);
    320:         PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);
    321:         lusol->indr = lusol->indc + nnz;
    322:         lusol->data = (double *)(lusol->indr + nnz);
    323:         lusol->nnz  = nnz;
    324:       }

    326:       /*******************************************************************/
    327:       /* Fill in the data for the problem.      (1-based Fortran style)  */
    328:       /*******************************************************************/

    330:       nz = 0;
    331:       for (i = 0; i < n; i++)
    332:         {
    333:           rs = a->i[i];
    334:           re = a->i[i+1];

    336:           while (rs < re)
    337:             {
    338:               if (a->a[rs] != 0.0)
    339:                 {
    340:                   lusol->indc[nz] = i + 1;
    341:                   lusol->indr[nz] = a->j[rs] + 1;
    342:                   lusol->data[nz] = a->a[rs];
    343:                   nz++;
    344:                 }
    345:               rs++;
    346:             }
    347:         }

    349:       /*******************************************************************/
    350:       /* Do the factorization.                                           */
    351:       /*******************************************************************/

    353:       LU1FAC(&m, &n, &nz, &nnz,
    354:              lusol->luparm, lusol->parmlu, lusol->data,
    355:              lusol->indc, lusol->indr, lusol->ip, lusol->iq,
    356:              lusol->lenc, lusol->lenr, lusol->locc, lusol->locr,
    357:              lusol->iploc, lusol->iqloc, lusol->ipinv,
    358:              lusol->iqinv, lusol->mnsw, &status);
    359:
    360:       switch(status)
    361:         {
    362:         case 0:                /* factored */
    363:           break;

    365:         case 7:                /* insufficient memory */
    366:           break;

    368:         case 1:
    369:         case -1:                /* singular */
    370:           SETERRQ(PETSC_ERR_LIB,"Singular matrix");

    372:         case 3:
    373:         case 4:                /* error conditions */
    374:           SETERRQ(PETSC_ERR_LIB,"matrix error");

    376:         default:                /* unknown condition */
    377:           SETERRQ(PETSC_ERR_LIB,"matrix unknown return code");
    378:         }

    380:       factorizations++;
    381:     } while (status == 7);
    382:   (*F)->assembled = PETSC_TRUE;
    383:   return(0);
    384: }

    388: PetscErrorCode MatLUFactorSymbolic_LUSOL(Mat A, IS r, IS c,MatFactorInfo *info, Mat *F) {
    389:   /************************************************************************/
    390:   /* Input                                                                */
    391:   /*     A  - matrix to factor                                            */
    392:   /*     r  - row permutation (ignored)                                   */
    393:   /*     c  - column permutation (ignored)                                */
    394:   /*                                                                      */
    395:   /* Output                                                               */
    396:   /*     F  - matrix storing the factorization;                           */
    397:   /************************************************************************/
    398:   Mat       B;
    399:   Mat_LUSOL *lusol;
    401:   int        i, m, n, nz, nnz;

    404:
    405:   /************************************************************************/
    406:   /* Check the arguments.                                                 */
    407:   /************************************************************************/

    409:   MatGetSize(A, &m, &n);
    410:   nz = ((Mat_SeqAIJ *)A->data)->nz;

    412:   /************************************************************************/
    413:   /* Create the factorization.                                            */
    414:   /************************************************************************/

    416:   MatCreate(A->comm,&B);
    417:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);
    418:   MatSetType(B,A->type_name);
    419:   MatSeqAIJSetPreallocation(B,0,PETSC_NULL);

    421:   B->ops->lufactornumeric = MatLUFactorNumeric_LUSOL;
    422:   B->ops->solve           = MatSolve_LUSOL;
    423:   B->factor               = FACTOR_LU;
    424:   lusol                   = (Mat_LUSOL*)(B->spptr);

    426:   /************************************************************************/
    427:   /* Initialize parameters                                                */
    428:   /************************************************************************/

    430:   for (i = 0; i < 30; i++)
    431:     {
    432:       lusol->luparm[i] = 0;
    433:       lusol->parmlu[i] = 0;
    434:     }

    436:   lusol->luparm[1] = -1;
    437:   lusol->luparm[2] = 5;
    438:   lusol->luparm[7] = 1;

    440:   lusol->parmlu[0] = 1 / Factorization_Tolerance;
    441:   lusol->parmlu[1] = 1 / Factorization_Tolerance;
    442:   lusol->parmlu[2] = Factorization_Small_Tolerance;
    443:   lusol->parmlu[3] = Factorization_Pivot_Tolerance;
    444:   lusol->parmlu[4] = Factorization_Pivot_Tolerance;
    445:   lusol->parmlu[5] = 3.0;
    446:   lusol->parmlu[6] = 0.3;
    447:   lusol->parmlu[7] = 0.6;

    449:   /************************************************************************/
    450:   /* Allocate the workspace needed by LUSOL.                              */
    451:   /************************************************************************/

    453:   lusol->elbowroom = PetscMax(lusol->elbowroom, info->fill);
    454:   nnz = PetscMax((int)(lusol->elbowroom*nz), 5*n);
    455:
    456:   lusol->n = n;
    457:   lusol->nz = nz;
    458:   lusol->nnz = nnz;
    459:   lusol->luroom = 1.75;

    461:   PetscMalloc(sizeof(int)*n,&lusol->ip);
    462:   PetscMalloc(sizeof(int)*n,&lusol->iq);
    463:   PetscMalloc(sizeof(int)*n,&lusol->lenc);
    464:   PetscMalloc(sizeof(int)*n,&lusol->lenr);
    465:   PetscMalloc(sizeof(int)*n,&lusol->locc);
    466:   PetscMalloc(sizeof(int)*n,&lusol->locr);
    467:   PetscMalloc(sizeof(int)*n,&lusol->iploc);
    468:   PetscMalloc(sizeof(int)*n,&lusol->iqloc);
    469:   PetscMalloc(sizeof(int)*n,&lusol->ipinv);
    470:   PetscMalloc(sizeof(int)*n,&lusol->iqinv);
    471:   PetscMalloc(sizeof(double)*n,&lusol->mnsw);
    472:   PetscMalloc(sizeof(double)*n,&lusol->mnsv);

    474:   PetscMalloc((sizeof(double)+2*sizeof(int))*nnz,&lusol->indc);
    475:   lusol->indr = lusol->indc + nnz;
    476:   lusol->data = (double *)(lusol->indr + nnz);
    477:   lusol->CleanUpLUSOL = PETSC_TRUE;
    478:   *F = B;
    479:   return(0);
    480: }

    485: PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_LUSOL(Mat A,const MatType type,MatReuse reuse,Mat *newmat)
    486: {
    488:   PetscInt       m, n;
    489:   Mat_LUSOL      *lusol;
    490:   Mat            B=*newmat;

    493:   MatGetSize(A, &m, &n);
    494:   if (m != n) {
    495:     SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
    496:   }
    497:   if (reuse == MAT_INITIAL_MATRIX) {
    498:     MatDuplicate(A,MAT_COPY_VALUES,&B);
    499:   }
    500:
    501:   PetscNew(Mat_LUSOL,&lusol);
    502:   lusol->MatDuplicate        = A->ops->duplicate;
    503:   lusol->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
    504:   lusol->MatDestroy          = A->ops->destroy;
    505:   lusol->CleanUpLUSOL        = PETSC_FALSE;

    507:   B->spptr                   = (void*)lusol;
    508:   B->ops->duplicate          = MatDuplicate_LUSOL;
    509:   B->ops->lufactorsymbolic   = MatLUFactorSymbolic_LUSOL;
    510:   B->ops->destroy            = MatDestroy_LUSOL;

    512:   PetscLogInfo((0,"MatConvert_SeqAIJ_LUSOL:Using LUSOL for LU factorization and solves.\n"));
    513:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_lusol_C",
    514:                                            "MatConvert_SeqAIJ_LUSOL",MatConvert_SeqAIJ_LUSOL);
    515:   PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_lusol_seqaij_C",
    516:                                            "MatConvert_LUSOL_SeqAIJ",MatConvert_LUSOL_SeqAIJ);
    517:   PetscObjectChangeTypeName((PetscObject)B,type);
    518:   *newmat = B;
    519:   return(0);
    520: }

    525: PetscErrorCode MatDuplicate_LUSOL(Mat A, MatDuplicateOption op, Mat *M) {
    527:   Mat_LUSOL *lu=(Mat_LUSOL *)A->spptr;
    529:   (*lu->MatDuplicate)(A,op,M);
    530:   PetscMemcpy((*M)->spptr,lu,sizeof(Mat_LUSOL));
    531:   return(0);
    532: }

    534: /*MC
    535:   MATLUSOL - MATLUSOL = "lusol" - A matrix type providing direct solvers (LU) for sequential matrices
    536:   via the external package LUSOL.

    538:   If LUSOL is installed (see the manual for
    539:   instructions on how to declare the existence of external packages),
    540:   a matrix type can be constructed which invokes LUSOL solvers.
    541:   After calling MatCreate(...,A), simply call MatSetType(A,MATLUSOL).
    542:   This matrix type is only supported for double precision real.

    544:   This matrix inherits from MATSEQAIJ.  As a result, MatSeqAIJSetPreallocation is
    545:   supported for this matrix type.  MatConvert can be called for a fast inplace conversion
    546:   to and from the MATSEQAIJ matrix type.

    548:   Options Database Keys:
    549: . -mat_type lusol - sets the matrix type to "lusol" during a call to MatSetFromOptions()

    551:    Level: beginner

    553: .seealso: PCLU
    554: M*/

    559: PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_LUSOL(Mat A)
    560: {

    564:   /* Change type name before calling MatSetType to force proper construction of SeqAIJ and LUSOL types */
    565:   PetscObjectChangeTypeName((PetscObject)A,MATLUSOL);
    566:   MatSetType(A,MATSEQAIJ);
    567:   MatConvert_SeqAIJ_LUSOL(A,MATLUSOL,MAT_REUSE_MATRIX,&A);
    568:   return(0);
    569: }