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 PetscCall(PetscLogObjectMemory((PetscObject)B, essl->lna * (2 * sizeof(int) + sizeof(PetscScalar)) + essl->naux * sizeof(PetscScalar))); 102 103 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 104 PetscFunctionReturn(0); 105 } 106 107 PetscErrorCode MatFactorGetSolverType_essl(Mat A, MatSolverType *type) { 108 PetscFunctionBegin; 109 *type = MATSOLVERESSL; 110 PetscFunctionReturn(0); 111 } 112 113 /*MC 114 MATSOLVERESSL - "essl" - Provides direct solvers, LU, for sequential matrices 115 via the external package ESSL. 116 117 If ESSL is installed (see the manual for 118 instructions on how to declare the existence of external packages), 119 120 Works with `MATSEQAIJ` matrices 121 122 Level: beginner 123 124 .seealso: `PCLU`, `PCFactorSetMatSolverType()`, `MatSolverType` 125 M*/ 126 127 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A, MatFactorType ftype, Mat *F) { 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(PetscNewLog(B, &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(0); 154 } 155 156 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void) { 157 PetscFunctionBegin; 158 PetscCall(MatSolverTypeRegister(MATSOLVERESSL, MATSEQAIJ, MAT_FACTOR_LU, MatGetFactor_seqaij_essl)); 159 PetscFunctionReturn(0); 160 } 161