Actual source code: essl.c

      1: /*
      2:         Provides an interface to the IBM RS6000 Essl sparse solver

      4: */
      5: #include <../src/mat/impls/aij/seq/aij.h>

      7: /* #include <essl.h> This doesn't work!  */

      9: PETSC_EXTERN void dgss(int *, int *, double *, int *, int *, int *, double *, double *, int *);
     10: PETSC_EXTERN void dgsf(int *, int *, int *, double *, int *, int *, int *, int *, double *, double *, double *, int *);

     12: typedef struct {
     13:   int          n, nz;
     14:   PetscScalar *a;
     15:   int         *ia;
     16:   int         *ja;
     17:   int          lna;
     18:   int          iparm[5];
     19:   PetscReal    rparm[5];
     20:   PetscReal    oparm[5];
     21:   PetscScalar *aux;
     22:   int          naux;

     24:   PetscBool CleanUpESSL;
     25: } Mat_Essl;

     27: static PetscErrorCode MatDestroy_Essl(Mat A)
     28: {
     29:   Mat_Essl *essl = (Mat_Essl *)A->data;

     31:   PetscFunctionBegin;
     32:   if (essl->CleanUpESSL) PetscCall(PetscFree4(essl->a, essl->aux, essl->ia, essl->ja));
     33:   PetscCall(PetscFree(A->data));
     34:   PetscFunctionReturn(PETSC_SUCCESS);
     35: }

     37: static PetscErrorCode MatSolve_Essl(Mat A, Vec b, Vec x)
     38: {
     39:   Mat_Essl    *essl = (Mat_Essl *)A->data;
     40:   PetscScalar *xx;
     41:   int          nessl, zero = 0;

     43:   PetscFunctionBegin;
     44:   PetscCall(PetscBLASIntCast(A->cmap->n, &nessl));
     45:   PetscCall(VecCopy(b, x));
     46:   PetscCall(VecGetArray(x, &xx));
     47:   dgss(&zero, &nessl, essl->a, essl->ia, essl->ja, &essl->lna, xx, essl->aux, &essl->naux);
     48:   PetscCall(VecRestoreArray(x, &xx));
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: static PetscErrorCode MatLUFactorNumeric_Essl(Mat F, Mat A, const MatFactorInfo *info)
     53: {
     54:   Mat_SeqAIJ *aa   = (Mat_SeqAIJ *)A->data;
     55:   Mat_Essl   *essl = (Mat_Essl *)F->data;
     56:   int         nessl, i, one = 1;

     58:   PetscFunctionBegin;
     59:   PetscCall(PetscBLASIntCast(A->rmap->n, &nessl));
     60:   /* copy matrix data into silly ESSL data structure (1-based Fortran style) */
     61:   for (i = 0; i < A->rmap->n + 1; i++) essl->ia[i] = aa->i[i] + 1;
     62:   for (i = 0; i < aa->nz; i++) essl->ja[i] = aa->j[i] + 1;

     64:   PetscCall(PetscArraycpy(essl->a, aa->a, aa->nz));

     66:   /* set Essl options */
     67:   essl->iparm[0] = 1;
     68:   essl->iparm[1] = 5;
     69:   essl->iparm[2] = 1;
     70:   essl->iparm[3] = 0;
     71:   essl->rparm[0] = 1.e-12;
     72:   essl->rparm[1] = 1.0;

     74:   PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)A)->prefix, "-matessl_lu_threshold", &essl->rparm[1], NULL));

     76:   dgsf(&one, &nessl, &essl->nz, essl->a, essl->ia, essl->ja, &essl->lna, essl->iparm, essl->rparm, essl->oparm, essl->aux, &essl->naux);

     78:   F->ops->solve   = MatSolve_Essl;
     79:   F->assembled    = PETSC_TRUE;
     80:   F->preallocated = PETSC_TRUE;
     81:   PetscFunctionReturn(PETSC_SUCCESS);
     82: }

     84: static PetscErrorCode MatLUFactorSymbolic_Essl(Mat B, Mat A, IS r, IS c, const MatFactorInfo *info)
     85: {
     86:   Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
     87:   Mat_Essl   *essl;
     88:   PetscReal   f = 1.0;

     90:   PetscFunctionBegin;
     91:   essl = (Mat_Essl *)B->data;

     93:   /* allocate the work arrays required by ESSL */
     94:   f = info->fill;
     95:   PetscCall(PetscBLASIntCast(a->nz, &essl->nz));
     96:   PetscCall(PetscBLASIntCast((PetscInt)(a->nz * f), &essl->lna));
     97:   PetscCall(PetscBLASIntCast(100 + 10 * A->rmap->n, &essl->naux));

     99:   /* since malloc is slow on IBM we try a single malloc */
    100:   PetscCall(PetscMalloc4(essl->lna, &essl->a, essl->naux, &essl->aux, essl->lna, &essl->ia, essl->lna, &essl->ja));

    102:   essl->CleanUpESSL = PETSC_TRUE;

    104:   B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: static PetscErrorCode MatFactorGetSolverType_essl(Mat A, MatSolverType *type)
    109: {
    110:   PetscFunctionBegin;
    111:   *type = MATSOLVERESSL;
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: /*MC
    116:   MATSOLVERESSL - "essl" - Provides direct solvers, LU, for sequential matrices via the external package ESSL.

    118:   Works with `MATSEQAIJ` matrices

    120:    Level: beginner

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

    125: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A, MatFactorType ftype, Mat *F)
    126: {
    127:   Mat       B;
    128:   Mat_Essl *essl;

    130:   PetscFunctionBegin;
    131:   PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square");
    132:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    133:   PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, A->rmap->n, A->cmap->n));
    134:   PetscCall(PetscStrallocpy("essl", &((PetscObject)B)->type_name));
    135:   PetscCall(MatSetUp(B));

    137:   PetscCall(PetscNew(&essl));

    139:   B->data                  = essl;
    140:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
    141:   B->ops->destroy          = MatDestroy_Essl;
    142:   B->ops->getinfo          = MatGetInfo_External;

    144:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_essl));

    146:   B->factortype = MAT_FACTOR_LU;
    147:   PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU]));
    148:   PetscCall(PetscFree(B->solvertype));
    149:   PetscCall(PetscStrallocpy(MATSOLVERESSL, &B->solvertype));

    151:   *F = B;
    152:   PetscFunctionReturn(PETSC_SUCCESS);
    153: }

    155: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_Essl(void)
    156: {
    157:   PetscFunctionBegin;
    158:   PetscCall(MatSolverTypeRegister(MATSOLVERESSL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_essl));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }