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