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 { 169 PetscErrorCode ierr; 170 Mat_Essl *essl=(Mat_Essl*)(A->spptr); 171 172 PetscFunctionBegin; 173 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 174 175 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 176 A->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 177 ierr = PetscLogInfo((0,"MatAssemblyEnd_Essl:Using ESSL for LU factorization and solves\n"));CHKERRQ(ierr); 178 PetscFunctionReturn(0); 179 } 180 181 EXTERN_C_BEGIN 182 #undef __FUNCT__ 183 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 184 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_Essl(Mat A,const MatType type,MatReuse reuse,Mat *newmat) 185 { 186 Mat B=*newmat; 187 PetscErrorCode ierr; 188 Mat_Essl *essl; 189 190 PetscFunctionBegin; 191 if (reuse == MAT_INITIAL_MATRIX) { 192 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 193 } 194 195 ierr = PetscNew(Mat_Essl,&essl);CHKERRQ(ierr); 196 essl->MatDuplicate = A->ops->duplicate; 197 essl->MatAssemblyEnd = A->ops->assemblyend; 198 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 199 essl->MatDestroy = A->ops->destroy; 200 essl->CleanUpESSL = PETSC_FALSE; 201 202 B->spptr = (void*)essl; 203 B->ops->duplicate = MatDuplicate_Essl; 204 B->ops->assemblyend = MatAssemblyEnd_Essl; 205 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 206 B->ops->destroy = MatDestroy_Essl; 207 208 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 209 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 210 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 211 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 212 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 213 214 *newmat = B; 215 PetscFunctionReturn(0); 216 } 217 EXTERN_C_END 218 219 #undef __FUNCT__ 220 #define __FUNCT__ "MatDuplicate_Essl" 221 PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M) 222 { 223 PetscErrorCode ierr; 224 Mat_Essl *lu = (Mat_Essl *)A->spptr; 225 226 PetscFunctionBegin; 227 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 228 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 /*MC 233 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 234 via the external package ESSL. 235 236 If ESSL is installed (see the manual for 237 instructions on how to declare the existence of external packages), 238 a matrix type can be constructed which invokes ESSL solvers. 239 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 240 This matrix type is only supported for double precision real. 241 242 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 243 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 244 the MATSEQAIJ type without data copy. 245 246 Options Database Keys: 247 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 248 249 Level: beginner 250 251 .seealso: PCLU 252 M*/ 253 254 EXTERN_C_BEGIN 255 #undef __FUNCT__ 256 #define __FUNCT__ "MatCreate_Essl" 257 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Essl(Mat A) 258 { 259 PetscErrorCode ierr; 260 261 PetscFunctionBegin; 262 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */ 263 ierr = PetscObjectChangeTypeName((PetscObject)A,MATESSL);CHKERRQ(ierr); 264 ierr = MatSetType(A,MATSEQAIJ); 265 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); 266 PetscFunctionReturn(0); 267 } 268 EXTERN_C_END 269