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