Actual source code: essl.c


  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:   PetscFunctionBegin;
 33:   if (essl->CleanUpESSL) PetscCall(PetscFree4(essl->a, essl->aux, essl->ia, essl->ja));
 34:   PetscCall(PetscFree(A->data));
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

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

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

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

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

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

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

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

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

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

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

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

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

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

103:   essl->CleanUpESSL = PETSC_TRUE;

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

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

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

119:   Works with `MATSEQAIJ` matrices

121:    Level: beginner

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

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

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

138:   PetscCall(PetscNew(&essl));

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

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

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

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

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