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: }