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 m,zero = 0; 54 55 PetscFunctionBegin; 56 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 57 ierr = VecCopy(b,x);CHKERRQ(ierr); 58 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 59 dgss(&zero,&A->cmap->n,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 i,one = 1; 72 73 PetscFunctionBegin; 74 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 75 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 76 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 77 78 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 79 80 /* set Essl options */ 81 essl->iparm[0] = 1; 82 essl->iparm[1] = 5; 83 essl->iparm[2] = 1; 84 essl->iparm[3] = 0; 85 essl->rparm[0] = 1.e-12; 86 essl->rparm[1] = 1.0; 87 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 88 89 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); 90 91 F->ops->solve = MatSolve_Essl; 92 (F)->assembled = PETSC_TRUE; 93 (F)->preallocated = PETSC_TRUE; 94 PetscFunctionReturn(0); 95 } 96 97 98 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 102 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 103 { 104 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 105 PetscErrorCode ierr; 106 Mat_Essl *essl; 107 PetscReal f = 1.0; 108 109 PetscFunctionBegin; 110 essl = (Mat_Essl *)(B->spptr); 111 112 /* allocate the work arrays required by ESSL */ 113 f = info->fill; 114 essl->nz = a->nz; 115 essl->lna = (int)a->nz*f; 116 essl->naux = 100 + 10*A->rmap->n; 117 118 /* since malloc is slow on IBM we try a single malloc */ 119 ierr = PetscMalloc4(essl->lna,PetscScalar,&essl->a,essl->naux,PetscScalar,&essl->aux,essl->lna,int,&essl->ia,essl->lna,int,&essl->ja);CHKERRQ(ierr); 120 essl->CleanUpESSL = PETSC_TRUE; 121 122 ierr = PetscLogObjectMemory(B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 123 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 124 PetscFunctionReturn(0); 125 } 126 127 EXTERN_C_BEGIN 128 #undef __FUNCT__ 129 #define __FUNCT__ "MatFactorGetSolverPackage_essl" 130 PetscErrorCode MatFactorGetSolverPackage_essl(Mat A,const MatSolverPackage *type) 131 { 132 PetscFunctionBegin; 133 *type = MATSOLVERESSL; 134 PetscFunctionReturn(0); 135 } 136 137 /*MC 138 MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices 139 via the external package ESSL. 140 141 If ESSL is installed (see the manual for 142 instructions on how to declare the existence of external packages), 143 144 Works with MATSEQAIJ matrices 145 146 Level: beginner 147 148 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 149 M*/ 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "MatGetFactor_seqaij_essl" 153 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 154 { 155 Mat B; 156 PetscErrorCode ierr; 157 Mat_Essl *essl; 158 159 PetscFunctionBegin; 160 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 161 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 162 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 163 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 164 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 165 166 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 167 B->spptr = essl; 168 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 169 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_essl",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 170 B->factortype = MAT_FACTOR_LU; 171 *F = B; 172 PetscFunctionReturn(0); 173 } 174 EXTERN_C_END 175