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 PetscErrorCode MatDestroy_Essl(Mat A) 29 { 30 PetscErrorCode ierr; 31 Mat_Essl *essl=(Mat_Essl*)A->data; 32 33 PetscFunctionBegin; 34 if (essl->CleanUpESSL) { 35 ierr = PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);CHKERRQ(ierr); 36 } 37 ierr = PetscFree(A->data);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 42 { 43 Mat_Essl *essl = (Mat_Essl*)A->data; 44 PetscScalar *xx; 45 PetscErrorCode ierr; 46 int nessl,zero = 0; 47 48 PetscFunctionBegin; 49 ierr = PetscBLASIntCast(A->cmap->n,&nessl);CHKERRQ(ierr); 50 ierr = VecCopy(b,x);CHKERRQ(ierr); 51 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 52 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 53 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 58 { 59 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data; 60 Mat_Essl *essl=(Mat_Essl*)(F)->data; 61 PetscErrorCode ierr; 62 int nessl,i,one = 1; 63 64 PetscFunctionBegin; 65 ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr); 66 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 67 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 68 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 69 70 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 71 72 /* set Essl options */ 73 essl->iparm[0] = 1; 74 essl->iparm[1] = 5; 75 essl->iparm[2] = 1; 76 essl->iparm[3] = 0; 77 essl->rparm[0] = 1.e-12; 78 essl->rparm[1] = 1.0; 79 80 ierr = PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);CHKERRQ(ierr); 81 82 dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 83 84 F->ops->solve = MatSolve_Essl; 85 (F)->assembled = PETSC_TRUE; 86 (F)->preallocated = PETSC_TRUE; 87 PetscFunctionReturn(0); 88 } 89 90 91 92 93 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 94 { 95 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 96 PetscErrorCode ierr; 97 Mat_Essl *essl; 98 PetscReal f = 1.0; 99 100 PetscFunctionBegin; 101 essl = (Mat_Essl*)(B->data); 102 103 /* allocate the work arrays required by ESSL */ 104 f = info->fill; 105 ierr = PetscBLASIntCast(a->nz,&essl->nz);CHKERRQ(ierr); 106 ierr = PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);CHKERRQ(ierr); 107 ierr = PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);CHKERRQ(ierr); 108 109 /* since malloc is slow on IBM we try a single malloc */ 110 ierr = PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);CHKERRQ(ierr); 111 112 essl->CleanUpESSL = PETSC_TRUE; 113 114 ierr = PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 115 116 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 117 PetscFunctionReturn(0); 118 } 119 120 PetscErrorCode MatFactorGetSolverPackage_essl(Mat A,const MatSolverPackage *type) 121 { 122 PetscFunctionBegin; 123 *type = MATSOLVERESSL; 124 PetscFunctionReturn(0); 125 } 126 127 /*MC 128 MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices 129 via the external package ESSL. 130 131 If ESSL is installed (see the manual for 132 instructions on how to declare the existence of external packages), 133 134 Works with MATSEQAIJ matrices 135 136 Level: beginner 137 138 .seealso: PCLU, PCFactorSetMatSolverPackage(), MatSolverPackage 139 M*/ 140 141 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 142 { 143 Mat B; 144 PetscErrorCode ierr; 145 Mat_Essl *essl; 146 147 PetscFunctionBegin; 148 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 149 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 150 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 151 ierr = PetscStrallocpy("essl",&((PetscObject)B)->type_name);CHKERRQ(ierr); 152 ierr = MatSetUp(B);CHKERRQ(ierr); 153 154 ierr = PetscNewLog(B,&essl);CHKERRQ(ierr); 155 156 B->data = essl; 157 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 158 B->ops->destroy = MatDestroy_Essl; 159 B->ops->getinfo = MatGetInfo_External; 160 161 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverPackage_C",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 162 163 B->factortype = MAT_FACTOR_LU; 164 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 165 ierr = PetscStrallocpy(MATSOLVERESSL,&B->solvertype);CHKERRQ(ierr); 166 167 *F = B; 168 PetscFunctionReturn(0); 169 } 170 171 PETSC_EXTERN PetscErrorCode MatSolverPackageRegister_Essl(void) 172 { 173 PetscErrorCode ierr; 174 PetscFunctionBegin; 175 ierr = MatSolverPackageRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);CHKERRQ(ierr); 176 PetscFunctionReturn(0); 177 } 178