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 Fortran 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 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 106 PetscFunctionReturn(0); 107 } 108 109 PetscErrorCode MatFactorGetSolverType_essl(Mat A, MatSolverType *type) 110 { 111 PetscFunctionBegin; 112 *type = MATSOLVERESSL; 113 PetscFunctionReturn(0); 114 } 115 116 /*MC 117 MATSOLVERESSL - "essl" - Provides direct solvers, LU, for sequential matrices 118 via the external package ESSL. 119 120 If ESSL is installed (see the manual for 121 instructions on how to declare the existence of external packages), 122 123 Works with `MATSEQAIJ` matrices 124 125 Level: beginner 126 127 .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType` 128 M*/ 129 130 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A, MatFactorType ftype, Mat *F) 131 { 132 Mat B; 133 Mat_Essl *essl; 134 135 PetscFunctionBegin; 136 PetscCheck(A->cmap->N == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix must be square"); 137 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 138 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, A->rmap->n, A->cmap->n)); 139 PetscCall(PetscStrallocpy("essl", &((PetscObject)B)->type_name)); 140 PetscCall(MatSetUp(B)); 141 142 PetscCall(PetscNew(&essl)); 143 144 B->data = essl; 145 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 146 B->ops->destroy = MatDestroy_Essl; 147 B->ops->getinfo = MatGetInfo_External; 148 149 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_essl)); 150 151 B->factortype = MAT_FACTOR_LU; 152 PetscCall(PetscStrallocpy(MATORDERINGEXTERNAL, (char **)&B->preferredordering[MAT_FACTOR_LU])); 153 PetscCall(PetscFree(B->solvertype)); 154 PetscCall(PetscStrallocpy(MATSOLVERESSL, &B->solvertype)); 155 156 *F = B; 157 PetscFunctionReturn(0); 158 } 159 160 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void) 161 { 162 PetscFunctionBegin; 163 PetscCall(MatSolverTypeRegister(MATSOLVERESSL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_essl)); 164 PetscFunctionReturn(0); 165 } 166