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 PETSC_EXTERN void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*); 11 PETSC_EXTERN void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*); 12 13 typedef struct { 14 int n,nz; 15 PetscScalar *a; 16 int *ia; 17 int *ja; 18 int lna; 19 int iparm[5]; 20 PetscReal rparm[5]; 21 PetscReal oparm[5]; 22 PetscScalar *aux; 23 int naux; 24 25 PetscBool CleanUpESSL; 26 } Mat_Essl; 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "MatDestroy_Essl" 30 PetscErrorCode MatDestroy_Essl(Mat A) 31 { 32 PetscErrorCode ierr; 33 Mat_Essl *essl=(Mat_Essl*)A->spptr; 34 35 PetscFunctionBegin; 36 if (essl && essl->CleanUpESSL) { 37 ierr = PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);CHKERRQ(ierr); 38 } 39 ierr = PetscFree(A->spptr);CHKERRQ(ierr); 40 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "MatSolve_Essl" 46 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 47 { 48 Mat_Essl *essl = (Mat_Essl*)A->spptr; 49 PetscScalar *xx; 50 PetscErrorCode ierr; 51 int nessl,zero = 0; 52 53 PetscFunctionBegin; 54 ierr = PetscBLASIntCast(A->cmap->n,&nessl);CHKERRQ(ierr); 55 ierr = VecCopy(b,x);CHKERRQ(ierr); 56 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 57 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 58 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 59 PetscFunctionReturn(0); 60 } 61 62 #undef __FUNCT__ 63 #define __FUNCT__ "MatLUFactorNumeric_Essl" 64 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 65 { 66 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data; 67 Mat_Essl *essl=(Mat_Essl*)(F)->spptr; 68 PetscErrorCode ierr; 69 int nessl,i,one = 1; 70 71 PetscFunctionBegin; 72 ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr); 73 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 74 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 75 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 76 77 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 78 79 /* set Essl options */ 80 essl->iparm[0] = 1; 81 essl->iparm[1] = 5; 82 essl->iparm[2] = 1; 83 essl->iparm[3] = 0; 84 essl->rparm[0] = 1.e-12; 85 essl->rparm[1] = 1.0; 86 87 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);CHKERRQ(ierr); 88 89 dgsf(&one,&nessl,&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 ierr = PetscBLASIntCast(a->nz,&essl->nz);CHKERRQ(ierr); 115 ierr = PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);CHKERRQ(ierr); 116 ierr = PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);CHKERRQ(ierr); 117 118 /* since malloc is slow on IBM we try a single malloc */ 119 ierr = PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);CHKERRQ(ierr); 120 121 essl->CleanUpESSL = PETSC_TRUE; 122 123 ierr = PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 124 125 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 126 PetscFunctionReturn(0); 127 } 128 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 PETSC_EXTERN 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(PetscObjectComm((PetscObject)A),&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,NULL);CHKERRQ(ierr); 166 167 ierr = PetscNewLog(B,&essl);CHKERRQ(ierr); 168 169 B->spptr = essl; 170 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 171 172 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 173 174 B->factortype = MAT_FACTOR_LU; 175 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 176 ierr = PetscStrallocpy(MATSOLVERESSL,&B->solvertype);CHKERRQ(ierr); 177 178 *F = B; 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "MatSolverPackageRegister_Essl" 184 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Essl(void) 185 { 186 PetscErrorCode ierr; 187 PetscFunctionBegin; 188 ierr = MatSolverPackageRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);CHKERRQ(ierr); 189 PetscFunctionReturn(0); 190 } 191