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 nessl = PetscBLASIntCast(A->cmap->n); 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 nessl = PetscBLASIntCast(A->rmap->n); 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,&nessl,&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 = PetscBLASIntCast(a->nz); 116 essl->lna = PetscBLASIntCast((PetscInt)(a->nz*f)); 117 essl->naux = PetscBLASIntCast(100 + 10*A->rmap->n); 118 119 /* since malloc is slow on IBM we try a single malloc */ 120 ierr = PetscMalloc4(essl->lna,PetscScalar,&essl->a,essl->naux,PetscScalar,&essl->aux,essl->lna,int,&essl->ia,essl->lna,int,&essl->ja);CHKERRQ(ierr); 121 essl->CleanUpESSL = PETSC_TRUE; 122 123 ierr = PetscLogObjectMemory(B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 124 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 125 PetscFunctionReturn(0); 126 } 127 128 EXTERN_C_BEGIN 129 #undef __FUNCT__ 130 #define __FUNCT__ "MatFactorGetSolverPackage_essl" 131 PetscErrorCode MatFactorGetSolverPackage_essl(Mat A,const MatSolverPackage *type) 132 { 133 PetscFunctionBegin; 134 *type = MATSOLVERESSL; 135 PetscFunctionReturn(0); 136 } 137 138 /*MC 139 MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices 140 via the external package ESSL. 141 142 If ESSL is installed (see the manual for 143 instructions on how to declare the existence of external packages), 144 145 Works with MATSEQAIJ matrices 146 147 Level: beginner 148 149 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 150 M*/ 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "MatGetFactor_seqaij_essl" 154 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 155 { 156 Mat B; 157 PetscErrorCode ierr; 158 Mat_Essl *essl; 159 160 PetscFunctionBegin; 161 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 162 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 163 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 164 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 165 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 166 167 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 168 B->spptr = essl; 169 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 170 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_essl",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 171 B->factortype = MAT_FACTOR_LU; 172 *F = B; 173 PetscFunctionReturn(0); 174 } 175 EXTERN_C_END 176