1 /*
2 Provides an interface to the IBM RS6000 Essl sparse solver
3
4 */
5 #include <../src/mat/impls/aij/seq/aij.h>
6
7 /* #include <essl.h> This doesn't work! */
8
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 *);
11
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;
23
24 PetscBool CleanUpESSL;
25 } Mat_Essl;
26
MatDestroy_Essl(Mat A)27 static PetscErrorCode MatDestroy_Essl(Mat A)
28 {
29 Mat_Essl *essl = (Mat_Essl *)A->data;
30
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 }
36
MatSolve_Essl(Mat A,Vec b,Vec x)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;
42
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 }
51
MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo * info)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;
57
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;
63
64 PetscCall(PetscArraycpy(essl->a, aa->a, aa->nz));
65
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;
73
74 PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)A)->prefix, "-matessl_lu_threshold", &essl->rparm[1], NULL));
75
76 dgsf(&one, &nessl, &essl->nz, essl->a, essl->ia, essl->ja, &essl->lna, essl->iparm, essl->rparm, essl->oparm, essl->aux, &essl->naux);
77
78 F->ops->solve = MatSolve_Essl;
79 F->assembled = PETSC_TRUE;
80 F->preallocated = PETSC_TRUE;
81 PetscFunctionReturn(PETSC_SUCCESS);
82 }
83
MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo * info)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;
89
90 PetscFunctionBegin;
91 essl = (Mat_Essl *)B->data;
92
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));
98
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));
101
102 essl->CleanUpESSL = PETSC_TRUE;
103
104 B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
MatFactorGetSolverType_essl(Mat A,MatSolverType * type)108 static PetscErrorCode MatFactorGetSolverType_essl(Mat A, MatSolverType *type)
109 {
110 PetscFunctionBegin;
111 *type = MATSOLVERESSL;
112 PetscFunctionReturn(PETSC_SUCCESS);
113 }
114
115 /*MC
116 MATSOLVERESSL - "essl" - Provides direct solvers, LU, for sequential matrices via the external package ESSL.
117
118 Works with `MATSEQAIJ` matrices
119
120 Level: beginner
121
122 .seealso: [](ch_matrices), `Mat`, `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType`
123 M*/
124
MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat * F)125 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A, MatFactorType ftype, Mat *F)
126 {
127 Mat B;
128 Mat_Essl *essl;
129
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));
136
137 PetscCall(PetscNew(&essl));
138
139 B->data = essl;
140 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
141 B->ops->destroy = MatDestroy_Essl;
142 B->ops->getinfo = MatGetInfo_External;
143
144 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_essl));
145
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));
150
151 *F = B;
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
MatSolverTypeRegister_Essl(void)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 }
161