1 /*$Id: essl.c,v 1.49 2001/08/07 03:02:47 balay Exp $*/ 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 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 23 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 24 int (*MatDestroy)(Mat); 25 PetscTruth CleanUpESSL; 26 } Mat_SeqAIJ_Essl; 27 28 EXTERN_C_BEGIN 29 #undef __FUNCT__ 30 #define __FUNCT__ "MatConvert_Essl_SeqAIJ" 31 int MatConvert_Essl_SeqAIJ(Mat A,MatType type,Mat *newmat) { 32 int ierr; 33 Mat B=*newmat; 34 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 35 36 PetscFunctionBegin; 37 if (B != A) { 38 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); 39 } else { 40 B->ops->matassemblyend = essl->MatAssemblyEnd; 41 B->ops->lufactorsymbolic = essl->MatLUFactorSymbolic; 42 B->ops->destroy = essl->MatDestroy; 43 44 /* free the Essl datastructures */ 45 ierr = PetscFree(essl);CHKERRQ(ierr); 46 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 47 } 48 *newmat = B; 49 PetscFunctionReturn(0); 50 } 51 EXTERN_C_END 52 53 #undef __FUNCT__ 54 #define __FUNCT__ "MatDestroy_SeqAIJ_Essl" 55 int MatDestroy_SeqAIJ_Essl(Mat A) 56 { 57 int ierr; 58 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 59 60 PetscFunctionBegin; 61 if (essl->CleanUpESSL) { 62 ierr = PetscFree(essl->a);CHKERRQ(ierr); 63 } 64 ierr = MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,&A); 65 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "MatSolve_SeqAIJ_Essl" 71 int MatSolve_SeqAIJ_Essl(Mat A,Vec b,Vec x) 72 { 73 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_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_SeqAIJ_Essl" 88 int MatLUFactorNumeric_SeqAIJ_Essl(Mat A,Mat *F) 89 { 90 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(A)->data; 91 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)(*F)->spptr; 92 int i,ierr,one = 1; 93 94 PetscFunctionBegin; 95 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 96 for (i=0; i<A->m+1; i++) essl->ia[i] = aa->i[i] + 1; 97 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 98 99 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 100 101 /* set Essl options */ 102 essl->iparm[0] = 1; 103 essl->iparm[1] = 5; 104 essl->iparm[2] = 1; 105 essl->iparm[3] = 0; 106 essl->rparm[0] = 1.e-12; 107 essl->rparm[1] = A->lupivotthreshold; 108 109 dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 110 essl->rparm,essl->oparm,essl->aux,&essl->naux); 111 112 PetscFunctionReturn(0); 113 } 114 115 #undef __FUNCT__ 116 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Essl" 117 int MatLUFactorSymbolic_SeqAIJ_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 118 { 119 Mat B; 120 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 121 int ierr,len; 122 Mat_SeqAIJ_Essl *essl; 123 PetscReal f = 1.0; 124 125 PetscFunctionBegin; 126 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 127 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 128 ierr = MatSetType(B,MATESSL);CHKERRQ(ierr); 129 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 130 131 B->ops->solve = MatSolve_SeqAIJ_Essl; 132 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Essl; 133 B->factor = FACTOR_LU; 134 135 essl = (Mat_SeqAIJ_Essl *)(B->spptr); 136 137 /* allocate the work arrays required by ESSL */ 138 f = info->fill; 139 essl->nz = a->nz; 140 essl->lna = (int)a->nz*f; 141 essl->naux = 100 + 10*A->m; 142 143 /* since malloc is slow on IBM we try a single malloc */ 144 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 145 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 146 essl->aux = essl->a + essl->lna; 147 essl->ia = (int*)(essl->aux + essl->naux); 148 essl->ja = essl->ia + essl->lna; 149 essl->CleanUpESSL = PETSC_TRUE; 150 151 PetscLogObjectMemory(B,len+sizeof(Mat_SeqAIJ_Essl)); 152 *F = B; 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ_Essl" 158 int MatAssemblyEnd_SeqAIJ_Essl(Mat A,MatAssemblyType mode) { 159 int ierr; 160 Mat_SeqAIJ_Essl *essl=(Mat_SeqAIJ_Essl*)(A->spptr); 161 162 PetscFunctionBegin; 163 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 164 165 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 166 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 167 PetscLogInfo(0,"Using ESSL for SeqAIJ LU factorization and solves"); 168 PetscFunctionReturn(0); 169 } 170 171 EXTERN_C_BEGIN 172 #undef __FUNCT__ 173 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 174 int MatConvert_SeqAIJ_Essl(Mat A,MatType type,Mat *newmat) { 175 Mat B=*newmat; 176 int ierr; 177 Mat_SeqAIJ_Essl *essl; 178 179 PetscFunctionBegin; 180 181 if (B != A) { 182 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 183 } 184 185 ierr = PetscNew(Mat_SeqAIJ_Essl,&essl);CHKERRQ(ierr); 186 essl->MatAssemblyEnd = A->ops->assemblyend; 187 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 188 essl->MatDestroy = A->ops->destroy; 189 essl->CleanUpESSL = PETSC_FALSE; 190 191 B->spptr = (void *)essl; 192 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ_Essl; 193 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 194 B->ops->destroy = MatDestroy_SeqAIJ_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 /*MC 207 MATESSL - a matrix type providing direct solvers (LU) for sequential matrices 208 via the external package ESSL. 209 210 If ESSL is installed (see the manual for 211 instructions on how to declare the existence of external packages), 212 a matrix type can be constructed which invokes ESSL solvers. 213 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 214 This matrix type is only supported for double precision real. 215 216 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 217 supported for this matrix type. 218 219 Options Database Keys: 220 . -mat_type essl - sets the matrix type to essl during a call to MatSetFromOptions() 221 222 Level: beginner 223 224 .seealso: PCLU 225 M*/ 226 227 EXTERN_C_BEGIN 228 #undef __FUNCT__ 229 #define __FUNCT__ "MatCreate_SeqAIJ_Essl" 230 int MatCreate_SeqAIJ_Essl(Mat A) { 231 int ierr; 232 233 PetscFunctionBegin; 234 ierr = MatSetType(A,MATSEQAIJ); 235 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 236 PetscFunctionReturn(0); 237 } 238 EXTERN_C_END 239