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 /* #include <essl.h> This doesn't work! */ 9 10 EXTERN_C_BEGIN 11 void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*); 12 void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*); 13 EXTERN_C_END 14 15 typedef struct { 16 int n,nz; 17 PetscScalar *a; 18 int *ia; 19 int *ja; 20 int lna; 21 int iparm[5]; 22 PetscReal rparm[5]; 23 PetscReal oparm[5]; 24 PetscScalar *aux; 25 int naux; 26 27 PetscErrorCode (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 28 PetscErrorCode (*MatAssemblyEnd)(Mat,MatAssemblyType); 29 PetscErrorCode (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 30 PetscErrorCode (*MatDestroy)(Mat); 31 PetscTruth CleanUpESSL; 32 } Mat_Essl; 33 34 EXTERN PetscErrorCode MatDuplicate_Essl(Mat,MatDuplicateOption,Mat*); 35 36 EXTERN_C_BEGIN 37 #undef __FUNCT__ 38 #define __FUNCT__ "MatConvert_Essl_SeqAIJ" 39 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_Essl_SeqAIJ(Mat A,MatType type,MatReuse reuse,Mat *newmat) { 40 PetscErrorCode ierr; 41 Mat B=*newmat; 42 Mat_Essl *essl=(Mat_Essl*)A->spptr; 43 44 PetscFunctionBegin; 45 if (reuse == MAT_INITIAL_MATRIX) { 46 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); 47 } 48 B->ops->duplicate = essl->MatDuplicate; 49 B->ops->assemblyend = essl->MatAssemblyEnd; 50 B->ops->lufactorsymbolic = essl->MatLUFactorSymbolic; 51 B->ops->destroy = essl->MatDestroy; 52 53 /* free the Essl datastructures */ 54 ierr = PetscFree(essl);CHKERRQ(ierr); 55 56 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_essl_C","",PETSC_NULL);CHKERRQ(ierr); 57 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_essl_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 58 59 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 60 *newmat = B; 61 PetscFunctionReturn(0); 62 } 63 EXTERN_C_END 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "MatDestroy_Essl" 67 PetscErrorCode MatDestroy_Essl(Mat A) 68 { 69 PetscErrorCode ierr; 70 Mat_Essl *essl=(Mat_Essl*)A->spptr; 71 72 PetscFunctionBegin; 73 if (essl->CleanUpESSL) { 74 ierr = PetscFree(essl->a);CHKERRQ(ierr); 75 } 76 ierr = MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,MAT_REUSE_MATRIX,&A); 77 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNCT__ 82 #define __FUNCT__ "MatSolve_Essl" 83 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) { 84 Mat_Essl *essl = (Mat_Essl*)A->spptr; 85 PetscScalar *xx; 86 PetscErrorCode ierr; 87 int m,zero = 0; 88 89 PetscFunctionBegin; 90 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 91 ierr = VecCopy(b,x);CHKERRQ(ierr); 92 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 93 dgss(&zero,&A->cmap.n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 94 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 #undef __FUNCT__ 99 #define __FUNCT__ "MatLUFactorNumeric_Essl" 100 PetscErrorCode MatLUFactorNumeric_Essl(Mat A,MatFactorInfo *info,Mat *F) { 101 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(A)->data; 102 Mat_Essl *essl=(Mat_Essl*)(*F)->spptr; 103 PetscErrorCode ierr; 104 int i,one = 1; 105 106 PetscFunctionBegin; 107 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 108 for (i=0; i<A->rmap.n+1; i++) essl->ia[i] = aa->i[i] + 1; 109 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 110 111 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 112 113 /* set Essl options */ 114 essl->iparm[0] = 1; 115 essl->iparm[1] = 5; 116 essl->iparm[2] = 1; 117 essl->iparm[3] = 0; 118 essl->rparm[0] = 1.e-12; 119 essl->rparm[1] = 1.0; 120 ierr = PetscOptionsGetReal(((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 121 122 dgsf(&one,&A->rmap.n,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 123 essl->rparm,essl->oparm,essl->aux,&essl->naux); 124 125 (*F)->assembled = PETSC_TRUE; 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 131 PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) { 132 Mat B; 133 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 134 PetscErrorCode ierr; 135 int len; 136 Mat_Essl *essl; 137 PetscReal f = 1.0; 138 139 PetscFunctionBegin; 140 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 141 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 142 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 143 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 144 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 145 146 B->ops->solve = MatSolve_Essl; 147 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 148 B->factor = FACTOR_LU; 149 150 essl = (Mat_Essl *)(B->spptr); 151 152 /* allocate the work arrays required by ESSL */ 153 f = info->fill; 154 essl->nz = a->nz; 155 essl->lna = (int)a->nz*f; 156 essl->naux = 100 + 10*A->rmap.n; 157 158 /* since malloc is slow on IBM we try a single malloc */ 159 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 160 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 161 essl->aux = essl->a + essl->lna; 162 essl->ia = (int*)(essl->aux + essl->naux); 163 essl->ja = essl->ia + essl->lna; 164 essl->CleanUpESSL = PETSC_TRUE; 165 166 ierr = PetscLogObjectMemory(B,len);CHKERRQ(ierr); 167 *F = B; 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "MatAssemblyEnd_Essl" 173 PetscErrorCode MatAssemblyEnd_Essl(Mat A,MatAssemblyType mode) 174 { 175 PetscErrorCode ierr; 176 Mat_Essl *essl=(Mat_Essl*)(A->spptr); 177 178 PetscFunctionBegin; 179 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 180 181 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 182 A->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 183 ierr = PetscInfo(A,"Using ESSL for LU factorization and solves\n");CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 EXTERN_C_BEGIN 188 #undef __FUNCT__ 189 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 190 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_Essl(Mat A,MatType type,MatReuse reuse,Mat *newmat) 191 { 192 Mat B=*newmat; 193 PetscErrorCode ierr; 194 Mat_Essl *essl; 195 196 PetscFunctionBegin; 197 if (reuse == MAT_INITIAL_MATRIX) { 198 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 199 } 200 B->ops->matsolve = 0; 201 202 ierr = PetscNewLog(B,Mat_Essl,&essl);CHKERRQ(ierr); 203 essl->MatDuplicate = A->ops->duplicate; 204 essl->MatAssemblyEnd = A->ops->assemblyend; 205 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 206 essl->MatDestroy = A->ops->destroy; 207 essl->CleanUpESSL = PETSC_FALSE; 208 209 B->spptr = (void*)essl; 210 B->ops->duplicate = MatDuplicate_Essl; 211 B->ops->assemblyend = MatAssemblyEnd_Essl; 212 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 213 B->ops->destroy = MatDestroy_Essl; 214 215 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 216 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 217 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 218 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 219 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 220 221 *newmat = B; 222 PetscFunctionReturn(0); 223 } 224 EXTERN_C_END 225 226 #undef __FUNCT__ 227 #define __FUNCT__ "MatDuplicate_Essl" 228 PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M) 229 { 230 PetscErrorCode ierr; 231 Mat_Essl *lu = (Mat_Essl *)A->spptr; 232 233 PetscFunctionBegin; 234 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 235 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 239 /*MC 240 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 241 via the external package ESSL. 242 243 If ESSL is installed (see the manual for 244 instructions on how to declare the existence of external packages), 245 a matrix type can be constructed which invokes ESSL solvers. 246 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 247 This matrix type is only supported for double precision real. 248 249 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 250 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 251 the MATSEQAIJ type without data copy. 252 253 Options Database Keys: 254 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 255 256 Level: beginner 257 258 .seealso: PCLU 259 M*/ 260 261 EXTERN_C_BEGIN 262 #undef __FUNCT__ 263 #define __FUNCT__ "MatCreate_Essl" 264 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Essl(Mat A) 265 { 266 PetscErrorCode ierr; 267 268 PetscFunctionBegin; 269 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 270 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 EXTERN_C_END 274