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