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 A,MatFactorInfo *info,Mat *F) 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)->assembled = PETSC_TRUE; 93 PetscFunctionReturn(0); 94 } 95 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 99 PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 100 { 101 Mat B; 102 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 103 PetscErrorCode ierr; 104 int len; 105 Mat_Essl *essl; 106 PetscReal f = 1.0; 107 108 PetscFunctionBegin; 109 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 110 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 111 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 112 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 113 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 114 115 B->ops->solve = MatSolve_Essl; 116 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 117 B->factor = FACTOR_LU; 118 119 essl = (Mat_Essl *)(B->spptr); 120 121 /* allocate the work arrays required by ESSL */ 122 f = info->fill; 123 essl->nz = a->nz; 124 essl->lna = (int)a->nz*f; 125 essl->naux = 100 + 10*A->rmap.n; 126 127 /* since malloc is slow on IBM we try a single malloc */ 128 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 129 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 130 essl->aux = essl->a + essl->lna; 131 essl->ia = (int*)(essl->aux + essl->naux); 132 essl->ja = essl->ia + essl->lna; 133 essl->CleanUpESSL = PETSC_TRUE; 134 135 ierr = PetscLogObjectMemory(B,len);CHKERRQ(ierr); 136 *F = B; 137 PetscFunctionReturn(0); 138 } 139 140 /*MC 141 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 142 via the external package ESSL. 143 144 If ESSL is installed (see the manual for 145 instructions on how to declare the existence of external packages), 146 a matrix type can be constructed which invokes ESSL solvers. 147 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 148 This matrix type is only supported for double precision real. 149 150 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 151 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 152 the MATSEQAIJ type without data copy. 153 154 Options Database Keys: 155 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 156 157 Level: beginner 158 159 .seealso: PCLU 160 M*/ 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatGetFactor_seqaij_essl" 164 PetscErrorCode MatGetFactor_seqaij_essl(Mat A,Mat *F) 165 { 166 Mat B; 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 168 PetscErrorCode ierr; 169 Mat_Essl *essl; 170 171 PetscFunctionBegin; 172 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 173 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 174 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 175 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 176 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 177 178 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 179 B->spptr = essl; 180 B->ops->solve = MatSolve_Essl; 181 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 182 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 183 B->factor = FACTOR_LU; 184 *F = B; 185 PetscFunctionReturn(0); 186 } 187