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