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 } Mat_SeqAIJ_Essl; 26 27 EXTERN int MatLUFactorSymbolic_SeqAIJ_Essl(Mat,IS,IS,MatFactorInfo*,Mat*); 28 29 EXTERN_C_BEGIN 30 #undef __FUNCT__ 31 #define __FUNCT__ "MatConvert_Essl_SeqAIJ" 32 int MatConvert_Essl_SeqAIJ(Mat A,MatType type,Mat *newmat) { 33 int ierr; 34 Mat B=*newmat; 35 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 36 37 PetscFunctionBegin; 38 if (B != A) { 39 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); 40 } else { 41 /* free the Essl datastructures */ 42 ierr = PetscFree(essl->a);CHKERRQ(ierr); 43 ierr = PetscFree(essl);CHKERRQ(ierr); 44 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 45 } 46 *newmat = B; 47 PetscFunctionReturn(0); 48 } 49 EXTERN_C_END 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "MatDestroy_SeqAIJ_Essl" 53 int MatDestroy_SeqAIJ_Essl(Mat A) 54 { 55 int ierr; 56 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 57 int (*destroy)(Mat)=essl->MatDestroy; 58 59 PetscFunctionBegin; 60 ierr = MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,&A); 61 ierr = (*destroy)(A);CHKERRQ(ierr); 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ_Essl" 67 int MatAssemblyEnd_SeqAIJ_Essl(Mat A,MatAssemblyType mode) { 68 int ierr; 69 Mat_SeqAIJ_Essl *essl=(Mat_SeqAIJ_Essl*)(A->spptr); 70 71 PetscFunctionBegin; 72 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 73 74 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 75 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 76 PetscLogInfo(0,"Using ESSL for SeqAIJ LU factorization and solves"); 77 PetscFunctionReturn(0); 78 } 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatSolve_SeqAIJ_Essl" 82 int MatSolve_SeqAIJ_Essl(Mat A,Vec b,Vec x) 83 { 84 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 85 PetscScalar *xx; 86 int ierr,m,zero = 0; 87 88 PetscFunctionBegin; 89 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 90 ierr = VecCopy(b,x);CHKERRQ(ierr); 91 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 92 dgss(&zero,&A->n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 93 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_Essl" 99 int MatLUFactorNumeric_SeqAIJ_Essl(Mat A,Mat *F) 100 { 101 Mat_SeqAIJ *a = (Mat_SeqAIJ*)(*F)->data; 102 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(A)->data; 103 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)(*F)->spptr; 104 int i,ierr,one = 1; 105 106 PetscFunctionBegin; 107 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 108 for (i=0; i<A->m+1; i++) essl->ia[i] = aa->i[i] + 1; 109 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 110 111 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 112 113 /* set Essl options */ 114 essl->iparm[0] = 1; 115 essl->iparm[1] = 5; 116 essl->iparm[2] = 1; 117 essl->iparm[3] = 0; 118 essl->rparm[0] = 1.e-12; 119 essl->rparm[1] = A->lupivotthreshold; 120 121 dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 122 essl->rparm,essl->oparm,essl->aux,&essl->naux); 123 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Essl" 129 int MatLUFactorSymbolic_SeqAIJ_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 130 { 131 Mat B; 132 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b; 133 int ierr,*ridx,*cidx,i,len; 134 Mat_SeqAIJ_Essl *essl; 135 PetscReal f = 1.0; 136 137 PetscFunctionBegin; 138 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 139 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 140 ierr = MatSetType(B,MATESSL);CHKERRQ(ierr); 141 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 142 143 B->ops->solve = MatSolve_SeqAIJ_Essl; 144 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Essl; 145 B->factor = FACTOR_LU; 146 147 essl = (Mat_SeqAIJ_Essl *)(B->spptr); 148 149 /* allocate the work arrays required by ESSL */ 150 f = info->fill; 151 essl->nz = a->nz; 152 essl->lna = (int)a->nz*f; 153 essl->naux = 100 + 10*A->m; 154 155 /* since malloc is slow on IBM we try a single malloc */ 156 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 157 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 158 essl->aux = essl->a + essl->lna; 159 essl->ia = (int*)(essl->aux + essl->naux); 160 essl->ja = essl->ia + essl->lna; 161 162 PetscLogObjectMemory(B,len+sizeof(Mat_SeqAIJ_Essl)); 163 *F = B; 164 PetscFunctionReturn(0); 165 } 166 167 EXTERN_C_BEGIN 168 #undef __FUNCT__ 169 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 170 int MatConvert_SeqAIJ_Essl(Mat A,MatType type,Mat *newmat) { 171 Mat B=*newmat; 172 int ierr; 173 Mat_SeqAIJ_Essl *essl; 174 175 PetscFunctionBegin; 176 177 if (B != A) { 178 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 179 } 180 181 ierr = PetscNew(Mat_SeqAIJ_Essl,&essl);CHKERRQ(ierr); 182 essl->MatAssemblyEnd = A->ops->assemblyend; 183 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 184 essl->MatDestroy = A->ops->destroy; 185 B->spptr = (void *)essl; 186 187 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ_Essl; 188 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 189 B->ops->destroy = MatDestroy_SeqAIJ_Essl; 190 191 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 192 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 193 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 194 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 195 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 196 *newmat = B; 197 PetscFunctionReturn(0); 198 } 199 EXTERN_C_END 200 201 EXTERN_C_BEGIN 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatCreate_SeqAIJ_Essl" 204 int MatCreate_SeqAIJ_Essl(Mat A) { 205 int ierr; 206 Mat_SeqAIJ_Essl *essl; 207 208 PetscFunctionBegin; 209 ierr = MatSetType(A,MATSEQAIJ); 210 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 EXTERN_C_END 214