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