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

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

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

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

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

 25:   PetscBool CleanUpESSL;
 26: } Mat_Essl;

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

 32:   if (essl->CleanUpESSL) PetscFree4(essl->a, essl->aux, essl->ia, essl->ja);
 33:   PetscFree(A->data);
 34:   return 0;
 35: }

 37: 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:   PetscBLASIntCast(A->cmap->n, &nessl);
 44:   VecCopy(b, x);
 45:   VecGetArray(x, &xx);
 46:   dgss(&zero, &nessl, essl->a, essl->ia, essl->ja, &essl->lna, xx, essl->aux, &essl->naux);
 47:   VecRestoreArray(x, &xx);
 48:   return 0;
 49: }

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

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

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

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

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

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

 76:   F->ops->solve     = MatSolve_Essl;
 77:   (F)->assembled    = PETSC_TRUE;
 78:   (F)->preallocated = PETSC_TRUE;
 79:   return 0;
 80: }

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

 88:   essl = (Mat_Essl *)(B->data);

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

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

 99:   essl->CleanUpESSL = PETSC_TRUE;

101:   B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
102:   return 0;
103: }

105: PetscErrorCode MatFactorGetSolverType_essl(Mat A, MatSolverType *type)
106: {
107:   *type = MATSOLVERESSL;
108:   return 0;
109: }

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

115:   If ESSL is installed (see the manual for
116:   instructions on how to declare the existence of external packages),

118:   Works with `MATSEQAIJ` matrices

120:    Level: beginner

122: .seealso: `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;

131:   MatCreate(PetscObjectComm((PetscObject)A), &B);
132:   MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, A->rmap->n, A->cmap->n);
133:   PetscStrallocpy("essl", &((PetscObject)B)->type_name);
134:   MatSetUp(B);

136:   PetscNew(&essl);

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

143:   PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_essl);

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

150:   *F = B;
151:   return 0;
152: }

154: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void)
155: {
156:   MatSolverTypeRegister(MATSOLVERESSL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_essl);
157:   return 0;
158: }