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 PetscCall(PetscFree4(essl->a,essl->aux,essl->ia,essl->ja)); 36 } 37 PetscCall(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 PetscCall(PetscBLASIntCast(A->cmap->n,&nessl)); 50 PetscCall(VecCopy(b,x)); 51 PetscCall(VecGetArray(x,&xx)); 52 dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 53 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(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 PetscCall(PetscBLASIntCast(a->nz,&essl->nz)); 103 PetscCall(PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna)); 104 PetscCall(PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux)); 105 106 /* since malloc is slow on IBM we try a single malloc */ 107 PetscCall(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 PetscCall(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 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&B)); 147 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n)); 148 PetscCall(PetscStrallocpy("essl",&((PetscObject)B)->type_name)); 149 PetscCall(MatSetUp(B)); 150 151 PetscCall(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 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl)); 159 160 B->factortype = MAT_FACTOR_LU; 161 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU])); 162 PetscCall(PetscFree(B->solvertype)); 163 PetscCall(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 PetscCall(MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl)); 174 PetscFunctionReturn(0); 175 } 176