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