Actual source code: essl.c
1: #define PETSCMAT_DLL
3: /*
4: Provides an interface to the IBM RS6000 Essl sparse solver
6: */
7: #include src/mat/impls/aij/seq/aij.h
8: /* #include <essl.h> This doesn't work! */
11: void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*);
12: void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*);
15: typedef struct {
16: int n,nz;
17: PetscScalar *a;
18: int *ia;
19: int *ja;
20: int lna;
21: int iparm[5];
22: PetscReal rparm[5];
23: PetscReal oparm[5];
24: PetscScalar *aux;
25: int naux;
27: PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*);
28: PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType);
29: PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*);
30: PetscErrorCode (*MatDestroy)(Mat);
31: PetscTruth CleanUpESSL;
32: } Mat_Essl;
34: EXTERN PetscErrorCode MatDuplicate_Essl(Mat,MatDuplicateOption,Mat*);
39: PetscErrorCode MatConvert_Essl_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) {
41: Mat B=*newmat;
42: Mat_Essl *essl=(Mat_Essl*)A->spptr;
43:
45: if (reuse == MAT_INITIAL_MATRIX) {
46: MatDuplicate(A,MAT_COPY_VALUES,&B);
47: }
48: B->ops->duplicate = essl->MatDuplicate;
49: B->ops->assemblyend = essl->MatAssemblyEnd;
50: B->ops->lufactorsymbolic = essl->MatLUFactorSymbolic;
51: B->ops->destroy = essl->MatDestroy;
53: /* free the Essl datastructures */
54: PetscFree(essl);
56: PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_essl_C","",PETSC_NULL);
57: PetscObjectComposeFunction((PetscObject)B,"MatConvert_essl_seqaij_C","",PETSC_NULL);
59: PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);
60: *newmat = B;
61: return(0);
62: }
67: PetscErrorCode MatDestroy_Essl(Mat A)
68: {
70: Mat_Essl *essl=(Mat_Essl*)A->spptr;
73: if (essl->CleanUpESSL) {
74: PetscFree(essl->a);
75: }
76: MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A);
77: (*A->ops->destroy)(A);
78: return(0);
79: }
83: PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) {
84: Mat_Essl *essl = (Mat_Essl*)A->spptr;
85: PetscScalar *xx;
87: int m,zero = 0;
90: VecGetLocalSize(b,&m);
91: VecCopy(b,x);
92: VecGetArray(x,&xx);
93: dgss(&zero,&A->cmap.n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux);
94: VecRestoreArray(x,&xx);
95: return(0);
96: }
100: PetscErrorCode MatLUFactorNumeric_Essl(Mat A,MatFactorInfo *info,Mat *F) {
101: Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(A)->data;
102: Mat_Essl *essl=(Mat_Essl*)(*F)->spptr;
104: int i,one = 1;
107: /* copy matrix data into silly ESSL data structure (1-based Frotran style) */
108: for (i=0; i<A->rmap.n+1; i++) essl->ia[i] = aa->i[i] + 1;
109: for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1;
110:
111: PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));
112:
113: /* set Essl options */
114: essl->iparm[0] = 1;
115: essl->iparm[1] = 5;
116: essl->iparm[2] = 1;
117: essl->iparm[3] = 0;
118: essl->rparm[0] = 1.e-12;
119: essl->rparm[1] = 1.0;
120: PetscOptionsGetReal(A->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);
122: dgsf(&one,&A->rmap.n,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,
123: essl->rparm,essl->oparm,essl->aux,&essl->naux);
125: (*F)->assembled = PETSC_TRUE;
126: return(0);
127: }
131: PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) {
132: Mat B;
133: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
135: int len;
136: Mat_Essl *essl;
137: PetscReal f = 1.0;
140: if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square");
141: MatCreate(A->comm,&B);
142: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);
143: MatSetType(B,A->type_name);
144: MatSeqAIJSetPreallocation(B,0,PETSC_NULL);
146: B->ops->solve = MatSolve_Essl;
147: B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
148: B->factor = FACTOR_LU;
149:
150: essl = (Mat_Essl *)(B->spptr);
152: /* allocate the work arrays required by ESSL */
153: f = info->fill;
154: essl->nz = a->nz;
155: essl->lna = (int)a->nz*f;
156: essl->naux = 100 + 10*A->rmap.n;
158: /* since malloc is slow on IBM we try a single malloc */
159: len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar);
160: PetscMalloc(len,&essl->a);
161: essl->aux = essl->a + essl->lna;
162: essl->ia = (int*)(essl->aux + essl->naux);
163: essl->ja = essl->ia + essl->lna;
164: essl->CleanUpESSL = PETSC_TRUE;
166: PetscLogObjectMemory(B,len+sizeof(Mat_Essl));
167: *F = B;
168: return(0);
169: }
173: PetscErrorCode MatAssemblyEnd_Essl(Mat A,MatAssemblyType mode)
174: {
176: Mat_Essl *essl=(Mat_Essl*)(A->spptr);
179: (*essl->MatAssemblyEnd)(A,mode);
181: essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
182: A->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
183: PetscInfo(0,"Using ESSL for LU factorization and solves\n");
184: return(0);
185: }
190: PetscErrorCode MatConvert_SeqAIJ_Essl(Mat A,MatType type,MatReuse reuse,Mat *newmat)
191: {
192: Mat B=*newmat;
194: Mat_Essl *essl;
197: if (reuse == MAT_INITIAL_MATRIX) {
198: MatDuplicate(A,MAT_COPY_VALUES,&B);
199: }
201: PetscNew(Mat_Essl,&essl);
202: essl->MatDuplicate = A->ops->duplicate;
203: essl->MatAssemblyEnd = A->ops->assemblyend;
204: essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic;
205: essl->MatDestroy = A->ops->destroy;
206: essl->CleanUpESSL = PETSC_FALSE;
208: B->spptr = (void*)essl;
209: B->ops->duplicate = MatDuplicate_Essl;
210: B->ops->assemblyend = MatAssemblyEnd_Essl;
211: B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
212: B->ops->destroy = MatDestroy_Essl;
214: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C",
215: "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);
216: PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C",
217: "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);
218: PetscObjectChangeTypeName((PetscObject)B,type);
220: *newmat = B;
221: return(0);
222: }
227: PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M)
228: {
230: Mat_Essl *lu = (Mat_Essl *)A->spptr;
233: (*lu->MatDuplicate)(A,op,M);
234: PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));
235: return(0);
236: }
238: /*MC
239: MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices
240: via the external package ESSL.
242: If ESSL is installed (see the manual for
243: instructions on how to declare the existence of external packages),
244: a matrix type can be constructed which invokes ESSL solvers.
245: After calling MatCreate(...,A), simply call MatSetType(A,MATESSL).
246: This matrix type is only supported for double precision real.
248: This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is
249: supported for this matrix type. One can also call MatConvert for an inplace conversion to or from
250: the MATSEQAIJ type without data copy.
252: Options Database Keys:
253: . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions()
255: Level: beginner
257: .seealso: PCLU
258: M*/
263: PetscErrorCode MatCreate_Essl(Mat A)
264: {
268: MatSetType(A,MATSEQAIJ);
269: MatConvert_SeqAIJ_Essl(A,MATESSL,MAT_REUSE_MATRIX,&A);
270: return(0);
271: }