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