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 F,Mat A,const MatFactorInfo *info) 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->ops->solve = MatSolve_Essl; 93 (F)->assembled = PETSC_TRUE; 94 (F)->preallocated = PETSC_TRUE; 95 PetscFunctionReturn(0); 96 } 97 98 99 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 103 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 104 { 105 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 106 PetscErrorCode ierr; 107 Mat_Essl *essl; 108 PetscReal f = 1.0; 109 110 PetscFunctionBegin; 111 essl = (Mat_Essl *)(B->spptr); 112 113 /* allocate the work arrays required by ESSL */ 114 f = info->fill; 115 essl->nz = a->nz; 116 essl->lna = (int)a->nz*f; 117 essl->naux = 100 + 10*A->rmap->n; 118 119 /* since malloc is slow on IBM we try a single malloc */ 120 ierr = PetscMalloc(essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar),&essl->a);CHKERRQ(ierr); 121 essl->aux = essl->a + essl->lna; 122 essl->ia = (int*)(essl->aux + essl->naux); 123 essl->ja = essl->ia + essl->lna; 124 essl->CleanUpESSL = PETSC_TRUE; 125 126 ierr = PetscLogObjectMemory(B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 127 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 128 PetscFunctionReturn(0); 129 } 130 131 EXTERN_C_BEGIN 132 #undef __FUNCT__ 133 #define __FUNCT__ "MatFactorGetSolverPackage_essl" 134 PetscErrorCode MatFactorGetSolverPackage_essl(Mat A,const MatSolverPackage *type) 135 { 136 PetscFunctionBegin; 137 *type = MAT_SOLVER_ESSL; 138 PetscFunctionReturn(0); 139 } 140 141 /*MC 142 MAT_SOLVER_ESSL - "essl" - Provides 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 148 Works with MATSEQAIJ matrices 149 150 Level: beginner 151 152 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 153 M*/ 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatGetFactor_seqaij_essl" 157 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 158 { 159 Mat B; 160 PetscErrorCode ierr; 161 Mat_Essl *essl; 162 163 PetscFunctionBegin; 164 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 165 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 166 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 167 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 168 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 169 170 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 171 B->spptr = essl; 172 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 173 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_essl",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 174 B->factor = MAT_FACTOR_LU; 175 *F = B; 176 PetscFunctionReturn(0); 177 } 178 EXTERN_C_END 179