1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the IBM RS6000 Essl sparse solver 5 6 */ 7 #include "src/mat/impls/aij/seq/aij.h" 8 9 /* #include <essl.h> This doesn't work! */ 10 11 EXTERN_C_BEGIN 12 void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*); 13 void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*); 14 EXTERN_C_END 15 16 typedef struct { 17 int n,nz; 18 PetscScalar *a; 19 int *ia; 20 int *ja; 21 int lna; 22 int iparm[5]; 23 PetscReal rparm[5]; 24 PetscReal oparm[5]; 25 PetscScalar *aux; 26 int naux; 27 28 PetscTruth CleanUpESSL; 29 } Mat_Essl; 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "MatDestroy_Essl" 33 PetscErrorCode MatDestroy_Essl(Mat A) 34 { 35 PetscErrorCode ierr; 36 Mat_Essl *essl=(Mat_Essl*)A->spptr; 37 38 PetscFunctionBegin; 39 if (essl->CleanUpESSL) { 40 ierr = PetscFree(essl->a);CHKERRQ(ierr); 41 } 42 ierr = PetscFree(essl);CHKERRQ(ierr); 43 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatSolve_Essl" 49 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 50 { 51 Mat_Essl *essl = (Mat_Essl*)A->spptr; 52 PetscScalar *xx; 53 PetscErrorCode ierr; 54 int m,zero = 0; 55 56 PetscFunctionBegin; 57 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 58 ierr = VecCopy(b,x);CHKERRQ(ierr); 59 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 60 dgss(&zero,&A->cmap.n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 61 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "MatLUFactorNumeric_Essl" 67 PetscErrorCode MatLUFactorNumeric_Essl(Mat A,MatFactorInfo *info,Mat *F) 68 { 69 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(A)->data; 70 Mat_Essl *essl=(Mat_Essl*)(*F)->spptr; 71 PetscErrorCode ierr; 72 int i,one = 1; 73 74 PetscFunctionBegin; 75 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 76 for (i=0; i<A->rmap.n+1; i++) essl->ia[i] = aa->i[i] + 1; 77 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 78 79 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 80 81 /* set Essl options */ 82 essl->iparm[0] = 1; 83 essl->iparm[1] = 5; 84 essl->iparm[2] = 1; 85 essl->iparm[3] = 0; 86 essl->rparm[0] = 1.e-12; 87 essl->rparm[1] = 1.0; 88 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 89 90 dgsf(&one,&A->rmap.n,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 91 92 (*F)->assembled = PETSC_TRUE; 93 (*F)->preallocated = PETSC_TRUE; 94 PetscFunctionReturn(0); 95 } 96 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 100 PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 101 { 102 Mat B; 103 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 104 PetscErrorCode ierr; 105 int len; 106 Mat_Essl *essl; 107 PetscReal f = 1.0; 108 109 PetscFunctionBegin; 110 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 111 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 112 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 113 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 114 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 115 116 B->ops->solve = MatSolve_Essl; 117 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 118 B->factor = MAT_FACTOR_LU; 119 120 essl = (Mat_Essl *)(B->spptr); 121 122 /* allocate the work arrays required by ESSL */ 123 f = info->fill; 124 essl->nz = a->nz; 125 essl->lna = (int)a->nz*f; 126 essl->naux = 100 + 10*A->rmap.n; 127 128 /* since malloc is slow on IBM we try a single malloc */ 129 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 130 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 131 essl->aux = essl->a + essl->lna; 132 essl->ia = (int*)(essl->aux + essl->naux); 133 essl->ja = essl->ia + essl->lna; 134 essl->CleanUpESSL = PETSC_TRUE; 135 136 ierr = PetscLogObjectMemory(B,len);CHKERRQ(ierr); 137 *F = B; 138 PetscFunctionReturn(0); 139 } 140 141 /*MC 142 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 143 via the external package ESSL. 144 145 If ESSL is installed (see the manual for 146 instructions on how to declare the existence of external packages), 147 a matrix type can be constructed which invokes ESSL solvers. 148 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 149 This matrix type is only supported for double precision real. 150 151 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 152 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 153 the MATSEQAIJ type without data copy. 154 155 Options Database Keys: 156 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 157 158 Level: beginner 159 160 .seealso: PCLU 161 M*/ 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatGetFactor_seqaij_essl" 165 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 166 { 167 Mat B; 168 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 169 PetscErrorCode ierr; 170 Mat_Essl *essl; 171 172 PetscFunctionBegin; 173 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 174 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 175 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 176 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 177 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 178 179 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 180 B->spptr = essl; 181 B->ops->solve = MatSolve_Essl; 182 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 183 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 184 B->factor = MAT_FACTOR_LU; 185 *F = B; 186 PetscFunctionReturn(0); 187 } 188