1 2 /* 3 Provides an interface to the IBM RS6000 Essl sparse solver 4 5 */ 6 #include <../src/mat/impls/aij/seq/aij.h> 7 8 /* #include <essl.h> This doesn't work! */ 9 10 EXTERN_C_BEGIN 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*); 13 EXTERN_C_END 14 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; 26 27 PetscBool CleanUpESSL; 28 } Mat_Essl; 29 30 #undef __FUNCT__ 31 #define __FUNCT__ "MatDestroy_Essl" 32 PetscErrorCode MatDestroy_Essl(Mat A) 33 { 34 PetscErrorCode ierr; 35 Mat_Essl *essl=(Mat_Essl*)A->spptr; 36 37 PetscFunctionBegin; 38 if (essl && essl->CleanUpESSL) { 39 ierr = PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);CHKERRQ(ierr); 40 } 41 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 42 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 #undef __FUNCT__ 47 #define __FUNCT__ "MatSolve_Essl" 48 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 49 { 50 Mat_Essl *essl = (Mat_Essl*)A->spptr; 51 PetscScalar *xx; 52 PetscErrorCode ierr; 53 int nessl,zero = 0; 54 55 PetscFunctionBegin; 56 ierr = PetscBLASIntCast(A->cmap->n,&nessl);CHKERRQ(ierr); 57 ierr = VecCopy(b,x);CHKERRQ(ierr); 58 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 59 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 60 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "MatLUFactorNumeric_Essl" 66 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 67 { 68 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data; 69 Mat_Essl *essl=(Mat_Essl*)(F)->spptr; 70 PetscErrorCode ierr; 71 int nessl,i,one = 1; 72 73 PetscFunctionBegin; 74 ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr); 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 89 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 90 91 dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 92 93 F->ops->solve = MatSolve_Essl; 94 (F)->assembled = PETSC_TRUE; 95 (F)->preallocated = PETSC_TRUE; 96 PetscFunctionReturn(0); 97 } 98 99 100 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 104 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 105 { 106 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 107 PetscErrorCode ierr; 108 Mat_Essl *essl; 109 PetscReal f = 1.0; 110 111 PetscFunctionBegin; 112 essl = (Mat_Essl*)(B->spptr); 113 114 /* allocate the work arrays required by ESSL */ 115 f = info->fill; 116 ierr = PetscBLASIntCast(a->nz,&essl->nz);CHKERRQ(ierr); 117 ierr = PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);CHKERRQ(ierr); 118 ierr = PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);CHKERRQ(ierr); 119 120 /* since malloc is slow on IBM we try a single malloc */ 121 ierr = PetscMalloc4(essl->lna,PetscScalar,&essl->a,essl->naux,PetscScalar,&essl->aux,essl->lna,int,&essl->ia,essl->lna,int,&essl->ja);CHKERRQ(ierr); 122 123 essl->CleanUpESSL = PETSC_TRUE; 124 125 ierr = PetscLogObjectMemory(B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 126 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 = MATSOLVERESSL; 138 PetscFunctionReturn(0); 139 } 140 141 /*MC 142 MATSOLVERESSL - "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_COMM_SELF,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 172 B->spptr = essl; 173 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 174 175 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_essl",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 176 177 B->factortype = MAT_FACTOR_LU; 178 *F = B; 179 PetscFunctionReturn(0); 180 } 181 EXTERN_C_END 182