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->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->m+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(A->prefix,"-matessl_lu_threshold",&essl->rparm[1],PETSC_NULL);CHKERRQ(ierr); 121 122 dgsf(&one,&A->m,&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->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 141 ierr = MatCreate(A->comm,&B);CHKERRQ(ierr); 142 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n);CHKERRQ(ierr); 143 ierr = MatSetType(B,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->m; 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+sizeof(Mat_Essl));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 = PetscVerboseInfo((0,"MatAssemblyEnd_Essl: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 201 ierr = PetscNew(Mat_Essl,&essl);CHKERRQ(ierr); 202 essl->MatDuplicate = A->ops->duplicate; 203 essl->MatAssemblyEnd = A->ops->assemblyend; 204 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 205 essl->MatDestroy = A->ops->destroy; 206 essl->CleanUpESSL = PETSC_FALSE; 207 208 B->spptr = (void*)essl; 209 B->ops->duplicate = MatDuplicate_Essl; 210 B->ops->assemblyend = MatAssemblyEnd_Essl; 211 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 212 B->ops->destroy = MatDestroy_Essl; 213 214 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 215 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 216 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 217 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 218 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 219 220 *newmat = B; 221 PetscFunctionReturn(0); 222 } 223 EXTERN_C_END 224 225 #undef __FUNCT__ 226 #define __FUNCT__ "MatDuplicate_Essl" 227 PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M) 228 { 229 PetscErrorCode ierr; 230 Mat_Essl *lu = (Mat_Essl *)A->spptr; 231 232 PetscFunctionBegin; 233 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 234 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));CHKERRQ(ierr); 235 PetscFunctionReturn(0); 236 } 237 238 /*MC 239 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 240 via the external package ESSL. 241 242 If ESSL is installed (see the manual for 243 instructions on how to declare the existence of external packages), 244 a matrix type can be constructed which invokes ESSL solvers. 245 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 246 This matrix type is only supported for double precision real. 247 248 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 249 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 250 the MATSEQAIJ type without data copy. 251 252 Options Database Keys: 253 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 254 255 Level: beginner 256 257 .seealso: PCLU 258 M*/ 259 260 EXTERN_C_BEGIN 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatCreate_Essl" 263 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Essl(Mat A) 264 { 265 PetscErrorCode ierr; 266 267 PetscFunctionBegin; 268 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */ 269 ierr = PetscObjectChangeTypeName((PetscObject)A,MATESSL);CHKERRQ(ierr); 270 ierr = MatSetType(A,MATSEQAIJ); 271 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 EXTERN_C_END 275