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