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,Mat *newmat) { 34 PetscErrorCode ierr; 35 Mat B=*newmat; 36 Mat_Essl *essl=(Mat_Essl*)A->spptr; 37 38 PetscFunctionBegin; 39 if (B != A) { 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,&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,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] = A->lupivotthreshold; 114 115 dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 116 essl->rparm,essl->oparm,essl->aux,&essl->naux); 117 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatLUFactorSymbolic_Essl" 123 PetscErrorCode MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) { 124 Mat B; 125 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 126 PetscErrorCode ierr; 127 int len; 128 Mat_Essl *essl; 129 PetscReal f = 1.0; 130 131 PetscFunctionBegin; 132 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 133 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 134 ierr = MatSetType(B,A->type_name);CHKERRQ(ierr); 135 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 136 137 B->ops->solve = MatSolve_Essl; 138 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 139 B->factor = FACTOR_LU; 140 141 essl = (Mat_Essl *)(B->spptr); 142 143 /* allocate the work arrays required by ESSL */ 144 f = info->fill; 145 essl->nz = a->nz; 146 essl->lna = (int)a->nz*f; 147 essl->naux = 100 + 10*A->m; 148 149 /* since malloc is slow on IBM we try a single malloc */ 150 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 151 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 152 essl->aux = essl->a + essl->lna; 153 essl->ia = (int*)(essl->aux + essl->naux); 154 essl->ja = essl->ia + essl->lna; 155 essl->CleanUpESSL = PETSC_TRUE; 156 157 PetscLogObjectMemory(B,len+sizeof(Mat_Essl)); 158 *F = B; 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatAssemblyEnd_Essl" 164 PetscErrorCode MatAssemblyEnd_Essl(Mat A,MatAssemblyType mode) { 165 PetscErrorCode ierr; 166 Mat_Essl *essl=(Mat_Essl*)(A->spptr); 167 168 PetscFunctionBegin; 169 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 170 171 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 172 A->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 173 PetscLogInfo(0,"Using ESSL for LU factorization and solves"); 174 PetscFunctionReturn(0); 175 } 176 177 EXTERN_C_BEGIN 178 #undef __FUNCT__ 179 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 180 PetscErrorCode MatConvert_SeqAIJ_Essl(Mat A,const MatType type,Mat *newmat) 181 { 182 Mat B=*newmat; 183 PetscErrorCode ierr; 184 Mat_Essl *essl; 185 186 PetscFunctionBegin; 187 if (B != A) { 188 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 189 } 190 191 ierr = PetscNew(Mat_Essl,&essl);CHKERRQ(ierr); 192 essl->MatDuplicate = A->ops->duplicate; 193 essl->MatAssemblyEnd = A->ops->assemblyend; 194 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 195 essl->MatDestroy = A->ops->destroy; 196 essl->CleanUpESSL = PETSC_FALSE; 197 198 B->spptr = (void*)essl; 199 B->ops->duplicate = MatDuplicate_Essl; 200 B->ops->assemblyend = MatAssemblyEnd_Essl; 201 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 202 B->ops->destroy = MatDestroy_Essl; 203 204 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 205 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 206 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 207 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 208 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 209 *newmat = B; 210 PetscFunctionReturn(0); 211 } 212 EXTERN_C_END 213 214 #undef __FUNCT__ 215 #define __FUNCT__ "MatDuplicate_Essl" 216 PetscErrorCode MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M) { 217 PetscErrorCode ierr; 218 Mat_Essl *lu = (Mat_Essl *)A->spptr; 219 220 PetscFunctionBegin; 221 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 222 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 226 /*MC 227 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 228 via the external package ESSL. 229 230 If ESSL is installed (see the manual for 231 instructions on how to declare the existence of external packages), 232 a matrix type can be constructed which invokes ESSL solvers. 233 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 234 This matrix type is only supported for double precision real. 235 236 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 237 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 238 the MATSEQAIJ type without data copy. 239 240 Options Database Keys: 241 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 242 243 Level: beginner 244 245 .seealso: PCLU 246 M*/ 247 248 EXTERN_C_BEGIN 249 #undef __FUNCT__ 250 #define __FUNCT__ "MatCreate_Essl" 251 PetscErrorCode MatCreate_Essl(Mat A) 252 { 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */ 257 ierr = PetscObjectChangeTypeName((PetscObject)A,MATESSL);CHKERRQ(ierr); 258 ierr = MatSetType(A,MATSEQAIJ); 259 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 EXTERN_C_END 263