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