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 CHKERRQ(PetscFree4(essl->a,essl->aux,essl->ia,essl->ja)); 36 } 37 CHKERRQ(PetscFree(A->data)); 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 CHKERRQ(PetscBLASIntCast(A->cmap->n,&nessl)); 50 CHKERRQ(VecCopy(b,x)); 51 CHKERRQ(VecGetArray(x,&xx)); 52 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 53 CHKERRQ(VecRestoreArray(x,&xx)); 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 CHKERRQ(PetscBLASIntCast(A->rmap->n,&nessl)); 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 CHKERRQ(PetscArraycpy(essl->a,aa->a,aa->nz)); 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 CHKERRQ(PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL)); 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 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 91 { 92 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 93 PetscErrorCode ierr; 94 Mat_Essl *essl; 95 PetscReal f = 1.0; 96 97 PetscFunctionBegin; 98 essl = (Mat_Essl*)(B->data); 99 100 /* allocate the work arrays required by ESSL */ 101 f = info->fill; 102 CHKERRQ(PetscBLASIntCast(a->nz,&essl->nz)); 103 CHKERRQ(PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna)); 104 CHKERRQ(PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux)); 105 106 /* since malloc is slow on IBM we try a single malloc */ 107 CHKERRQ(PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja)); 108 109 essl->CleanUpESSL = PETSC_TRUE; 110 111 CHKERRQ(PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar))); 112 113 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFactorGetSolverType_essl(Mat A,MatSolverType *type) 118 { 119 PetscFunctionBegin; 120 *type = MATSOLVERESSL; 121 PetscFunctionReturn(0); 122 } 123 124 /*MC 125 MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices 126 via the external package ESSL. 127 128 If ESSL is installed (see the manual for 129 instructions on how to declare the existence of external packages), 130 131 Works with MATSEQAIJ matrices 132 133 Level: beginner 134 135 .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType 136 M*/ 137 138 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 139 { 140 Mat B; 141 PetscErrorCode ierr; 142 Mat_Essl *essl; 143 144 PetscFunctionBegin; 145 PetscCheckFalse(A->cmap->N != A->rmap->N,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 146 CHKERRQ(MatCreate(PetscObjectComm((PetscObject)A),&B)); 147 CHKERRQ(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n)); 148 CHKERRQ(PetscStrallocpy("essl",&((PetscObject)B)->type_name)); 149 CHKERRQ(MatSetUp(B)); 150 151 CHKERRQ(PetscNewLog(B,&essl)); 152 153 B->data = essl; 154 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 155 B->ops->destroy = MatDestroy_Essl; 156 B->ops->getinfo = MatGetInfo_External; 157 158 CHKERRQ(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl)); 159 160 B->factortype = MAT_FACTOR_LU; 161 CHKERRQ(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU])); 162 CHKERRQ(PetscFree(B->solvertype)); 163 CHKERRQ(PetscStrallocpy(MATSOLVERESSL,&B->solvertype)); 164 165 *F = B; 166 PetscFunctionReturn(0); 167 } 168 169 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void) 170 { 171 PetscErrorCode ierr; 172 PetscFunctionBegin; 173 CHKERRQ(MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl)); 174 PetscFunctionReturn(0); 175 } 176