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