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->data; 34 35 PetscFunctionBegin; 36 if (essl->CleanUpESSL) { 37 ierr = PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);CHKERRQ(ierr); 38 } 39 ierr = PetscFree(A->data);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 43 #undef __FUNCT__ 44 #define __FUNCT__ "MatSolve_Essl" 45 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 46 { 47 Mat_Essl *essl = (Mat_Essl*)A->data; 48 PetscScalar *xx; 49 PetscErrorCode ierr; 50 int nessl,zero = 0; 51 52 PetscFunctionBegin; 53 ierr = PetscBLASIntCast(A->cmap->n,&nessl);CHKERRQ(ierr); 54 ierr = VecCopy(b,x);CHKERRQ(ierr); 55 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 56 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 57 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "MatLUFactorNumeric_Essl" 63 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 64 { 65 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data; 66 Mat_Essl *essl=(Mat_Essl*)(F)->data; 67 PetscErrorCode ierr; 68 int nessl,i,one = 1; 69 70 PetscFunctionBegin; 71 ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr); 72 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 73 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 74 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 75 76 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 77 78 /* set Essl options */ 79 essl->iparm[0] = 1; 80 essl->iparm[1] = 5; 81 essl->iparm[2] = 1; 82 essl->iparm[3] = 0; 83 essl->rparm[0] = 1.e-12; 84 essl->rparm[1] = 1.0; 85 86 ierr = PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);CHKERRQ(ierr); 87 88 dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 89 90 F->ops->solve = MatSolve_Essl; 91 (F)->assembled = PETSC_TRUE; 92 (F)->preallocated = PETSC_TRUE; 93 PetscFunctionReturn(0); 94 } 95 96 97 98 99 #undef __FUNCT__ 100 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 101 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 102 { 103 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 104 PetscErrorCode ierr; 105 Mat_Essl *essl; 106 PetscReal f = 1.0; 107 108 PetscFunctionBegin; 109 essl = (Mat_Essl*)(B->data); 110 111 /* allocate the work arrays required by ESSL */ 112 f = info->fill; 113 ierr = PetscBLASIntCast(a->nz,&essl->nz);CHKERRQ(ierr); 114 ierr = PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);CHKERRQ(ierr); 115 ierr = PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);CHKERRQ(ierr); 116 117 /* since malloc is slow on IBM we try a single malloc */ 118 ierr = PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);CHKERRQ(ierr); 119 120 essl->CleanUpESSL = PETSC_TRUE; 121 122 ierr = PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 123 124 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 125 PetscFunctionReturn(0); 126 } 127 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 PETSC_EXTERN 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(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 162 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 163 ierr = PetscStrallocpy("essl",&((PetscObject)B)->type_name);CHKERRQ(ierr); 164 ierr = MatSetUp(B);CHKERRQ(ierr); 165 166 ierr = PetscNewLog(B,&essl);CHKERRQ(ierr); 167 168 B->data = essl; 169 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 170 B->ops->destroy = MatDestroy_Essl; 171 B->ops->getinfo = MatGetInfo_External; 172 173 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 174 175 B->factortype = MAT_FACTOR_LU; 176 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 177 ierr = PetscStrallocpy(MATSOLVERESSL,&B->solvertype);CHKERRQ(ierr); 178 179 *F = B; 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "MatSolverPackageRegister_Essl" 185 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Essl(void) 186 { 187 PetscErrorCode ierr; 188 PetscFunctionBegin; 189 ierr = MatSolverPackageRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192