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