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: {
31: Mat_Essl *essl=(Mat_Essl*)A->data;
34: if (essl->CleanUpESSL) {
35: PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);
36: }
37: PetscFree(A->data);
38: return(0);
39: }
41: PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x)
42: {
43: Mat_Essl *essl = (Mat_Essl*)A->data;
44: PetscScalar *xx;
46: int nessl,zero = 0;
49: PetscBLASIntCast(A->cmap->n,&nessl);
50: VecCopy(b,x);
51: VecGetArray(x,&xx);
52: dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux);
53: VecRestoreArray(x,&xx);
54: return(0);
55: }
57: PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info)
58: {
59: Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data;
60: Mat_Essl *essl=(Mat_Essl*)(F)->data;
62: int nessl,i,one = 1;
65: PetscBLASIntCast(A->rmap->n,&nessl);
66: /* copy matrix data into silly ESSL data structure (1-based Frotran style) */
67: for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1;
68: for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1;
70: PetscArraycpy(essl->a,aa->a,aa->nz);
72: /* set Essl options */
73: essl->iparm[0] = 1;
74: essl->iparm[1] = 5;
75: essl->iparm[2] = 1;
76: essl->iparm[3] = 0;
77: essl->rparm[0] = 1.e-12;
78: essl->rparm[1] = 1.0;
80: PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);
82: dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux);
84: F->ops->solve = MatSolve_Essl;
85: (F)->assembled = PETSC_TRUE;
86: (F)->preallocated = PETSC_TRUE;
87: return(0);
88: }
90: PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info)
91: {
92: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
94: Mat_Essl *essl;
95: PetscReal f = 1.0;
98: essl = (Mat_Essl*)(B->data);
100: /* allocate the work arrays required by ESSL */
101: f = info->fill;
102: PetscBLASIntCast(a->nz,&essl->nz);
103: PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);
104: PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);
106: /* since malloc is slow on IBM we try a single malloc */
107: PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);
109: essl->CleanUpESSL = PETSC_TRUE;
111: PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));
113: B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
114: return(0);
115: }
117: PetscErrorCode MatFactorGetSolverType_essl(Mat A,MatSolverType *type)
118: {
120: *type = MATSOLVERESSL;
121: return(0);
122: }
124: /*MC
125: MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices
126: via the external package ESSL.
128: If ESSL is installed (see the manual for
129: instructions on how to declare the existence of external packages),
131: Works with MATSEQAIJ matrices
133: Level: beginner
135: .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType
136: M*/
138: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F)
139: {
140: Mat B;
142: Mat_Essl *essl;
145: if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square");
146: MatCreate(PetscObjectComm((PetscObject)A),&B);
147: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);
148: PetscStrallocpy("essl",&((PetscObject)B)->type_name);
149: MatSetUp(B);
151: PetscNewLog(B,&essl);
153: B->data = essl;
154: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
155: B->ops->destroy = MatDestroy_Essl;
156: B->ops->getinfo = MatGetInfo_External;
158: PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl);
160: B->factortype = MAT_FACTOR_LU;
161: PetscFree(B->solvertype);
162: PetscStrallocpy(MATSOLVERESSL,&B->solvertype);
164: *F = B;
165: return(0);
166: }
168: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void)
169: {
172: MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);
173: return(0);
174: }