1 #define PETSCMAT_DLL 2 3 /* 4 Provides an interface to the IBM RS6000 Essl sparse solver 5 6 */ 7 #include "../src/mat/impls/aij/seq/aij.h" 8 9 /* #include <essl.h> This doesn't work! */ 10 11 EXTERN_C_BEGIN 12 void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*); 13 void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*); 14 EXTERN_C_END 15 16 typedef struct { 17 int n,nz; 18 PetscScalar *a; 19 int *ia; 20 int *ja; 21 int lna; 22 int iparm[5]; 23 PetscReal rparm[5]; 24 PetscReal oparm[5]; 25 PetscScalar *aux; 26 int naux; 27 28 PetscTruth CleanUpESSL; 29 } Mat_Essl; 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "MatDestroy_Essl" 33 PetscErrorCode MatDestroy_Essl(Mat A) 34 { 35 PetscErrorCode ierr; 36 Mat_Essl *essl=(Mat_Essl*)A->spptr; 37 38 PetscFunctionBegin; 39 if (essl->CleanUpESSL) { 40 ierr = PetscFree(essl->a);CHKERRQ(ierr); 41 } 42 ierr = PetscFree(essl);CHKERRQ(ierr); 43 ierr = MatDestroy_SeqAIJ(A);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatSolve_Essl" 49 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 50 { 51 Mat_Essl *essl = (Mat_Essl*)A->spptr; 52 PetscScalar *xx; 53 PetscErrorCode ierr; 54 int m,zero = 0; 55 56 PetscFunctionBegin; 57 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 58 ierr = VecCopy(b,x);CHKERRQ(ierr); 59 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 60 dgss(&zero,&A->cmap->n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 61 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "MatLUFactorNumeric_Essl" 67 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 68 { 69 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(A)->data; 70 Mat_Essl *essl=(Mat_Essl*)(F)->spptr; 71 PetscErrorCode ierr; 72 int i,one = 1; 73 74 PetscFunctionBegin; 75 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 76 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 77 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 78 79 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 80 81 /* set Essl options */ 82 essl->iparm[0] = 1; 83 essl->iparm[1] = 5; 84 essl->iparm[2] = 1; 85 essl->iparm[3] = 0; 86 essl->rparm[0] = 1.e-12; 87 essl->rparm[1] = 1.0; 88 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 89 90 dgsf(&one,&A->rmap->n,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 91 92 F->ops->solve = MatSolve_Essl; 93 (F)->assembled = PETSC_TRUE; 94 (F)->preallocated = PETSC_TRUE; 95 PetscFunctionReturn(0); 96 } 97 98 99 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 103 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 104 { 105 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 106 PetscErrorCode ierr; 107 int len; 108 Mat_Essl *essl; 109 PetscReal f = 1.0; 110 111 PetscFunctionBegin; 112 essl = (Mat_Essl *)(B->spptr); 113 114 /* allocate the work arrays required by ESSL */ 115 f = info->fill; 116 essl->nz = a->nz; 117 essl->lna = (int)a->nz*f; 118 essl->naux = 100 + 10*A->rmap->n; 119 120 /* since malloc is slow on IBM we try a single malloc */ 121 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 122 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 123 essl->aux = essl->a + essl->lna; 124 essl->ia = (int*)(essl->aux + essl->naux); 125 essl->ja = essl->ia + essl->lna; 126 essl->CleanUpESSL = PETSC_TRUE; 127 128 ierr = PetscLogObjectMemory(B,len);CHKERRQ(ierr); 129 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 130 PetscFunctionReturn(0); 131 } 132 133 EXTERN_C_BEGIN 134 #undef __FUNCT__ 135 #define __FUNCT__ "MatFactorGetSolverPackage_essl" 136 PetscErrorCode MatFactorGetSolverPackage_essl(Mat A,const MatSolverPackage *type) 137 { 138 PetscFunctionBegin; 139 *type = MAT_SOLVER_ESSL; 140 PetscFunctionReturn(0); 141 } 142 EXTERN_C_END 143 144 /*MC 145 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 146 via the external package ESSL. 147 148 If ESSL is installed (see the manual for 149 instructions on how to declare the existence of external packages), 150 a matrix type can be constructed which invokes ESSL solvers. 151 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 152 This matrix type is only supported for double precision real. 153 154 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 155 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 156 the MATSEQAIJ type without data copy. 157 158 Options Database Keys: 159 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 160 161 Level: beginner 162 163 .seealso: PCLU 164 M*/ 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "MatGetFactor_seqaij_essl" 168 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 169 { 170 Mat B; 171 PetscErrorCode ierr; 172 Mat_Essl *essl; 173 174 PetscFunctionBegin; 175 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 176 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 177 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 178 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 179 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 180 181 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 182 B->spptr = essl; 183 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 184 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_essl",MatFactorGetSolverPackage_essl);CHKERRQ(ierr); 185 B->factor = MAT_FACTOR_LU; 186 *F = B; 187 PetscFunctionReturn(0); 188 } 189