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) {
33: PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);
34: }
35: PetscFree(A->data);
36: return 0;
37: }
39: PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x)
40: {
41: Mat_Essl *essl = (Mat_Essl*)A->data;
42: PetscScalar *xx;
43: int nessl,zero = 0;
45: PetscBLASIntCast(A->cmap->n,&nessl);
46: VecCopy(b,x);
47: VecGetArray(x,&xx);
48: dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux);
49: VecRestoreArray(x,&xx);
50: return 0;
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: PetscBLASIntCast(A->rmap->n,&nessl);
60: /* copy matrix data into silly ESSL data structure (1-based Frotran 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;
64: PetscArraycpy(essl->a,aa->a,aa->nz);
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;
74: PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);
76: dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux);
78: F->ops->solve = MatSolve_Essl;
79: (F)->assembled = PETSC_TRUE;
80: (F)->preallocated = PETSC_TRUE;
81: return 0;
82: }
84: 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;
90: essl = (Mat_Essl*)(B->data);
92: /* allocate the work arrays required by ESSL */
93: f = info->fill;
94: PetscBLASIntCast(a->nz,&essl->nz);
95: PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);
96: PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);
98: /* since malloc is slow on IBM we try a single malloc */
99: PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);
101: essl->CleanUpESSL = PETSC_TRUE;
103: PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));
105: B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
106: return 0;
107: }
109: PetscErrorCode MatFactorGetSolverType_essl(Mat A,MatSolverType *type)
110: {
111: *type = MATSOLVERESSL;
112: return 0;
113: }
115: /*MC
116: MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices
117: via the external package ESSL.
119: If ESSL is installed (see the manual for
120: instructions on how to declare the existence of external packages),
122: Works with MATSEQAIJ matrices
124: Level: beginner
126: .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType
127: M*/
129: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F)
130: {
131: Mat B;
132: Mat_Essl *essl;
135: MatCreate(PetscObjectComm((PetscObject)A),&B);
136: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);
137: PetscStrallocpy("essl",&((PetscObject)B)->type_name);
138: MatSetUp(B);
140: PetscNewLog(B,&essl);
142: B->data = essl;
143: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
144: B->ops->destroy = MatDestroy_Essl;
145: B->ops->getinfo = MatGetInfo_External;
147: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl);
149: B->factortype = MAT_FACTOR_LU;
150: PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU]);
151: PetscFree(B->solvertype);
152: PetscStrallocpy(MATSOLVERESSL,&B->solvertype);
154: *F = B;
155: return 0;
156: }
158: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void)
159: {
160: MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);
161: return 0;
162: }