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 *aa = (Mat_SeqAIJ*)(A)->data; 102 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)(*F)->spptr; 103 int i,ierr,one = 1; 104 105 PetscFunctionBegin; 106 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 107 for (i=0; i<A->m+1; i++) essl->ia[i] = aa->i[i] + 1; 108 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 109 110 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 111 112 /* set Essl options */ 113 essl->iparm[0] = 1; 114 essl->iparm[1] = 5; 115 essl->iparm[2] = 1; 116 essl->iparm[3] = 0; 117 essl->rparm[0] = 1.e-12; 118 essl->rparm[1] = A->lupivotthreshold; 119 120 dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 121 essl->rparm,essl->oparm,essl->aux,&essl->naux); 122 123 PetscFunctionReturn(0); 124 } 125 126 #undef __FUNCT__ 127 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Essl" 128 int MatLUFactorSymbolic_SeqAIJ_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 129 { 130 Mat B; 131 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 132 int ierr,len; 133 Mat_SeqAIJ_Essl *essl; 134 PetscReal f = 1.0; 135 136 PetscFunctionBegin; 137 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 138 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 139 ierr = MatSetType(B,MATESSL);CHKERRQ(ierr); 140 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 141 142 B->ops->solve = MatSolve_SeqAIJ_Essl; 143 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Essl; 144 B->factor = FACTOR_LU; 145 146 essl = (Mat_SeqAIJ_Essl *)(B->spptr); 147 148 /* allocate the work arrays required by ESSL */ 149 f = info->fill; 150 essl->nz = a->nz; 151 essl->lna = (int)a->nz*f; 152 essl->naux = 100 + 10*A->m; 153 154 /* since malloc is slow on IBM we try a single malloc */ 155 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 156 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 157 essl->aux = essl->a + essl->lna; 158 essl->ia = (int*)(essl->aux + essl->naux); 159 essl->ja = essl->ia + essl->lna; 160 161 PetscLogObjectMemory(B,len+sizeof(Mat_SeqAIJ_Essl)); 162 *F = B; 163 PetscFunctionReturn(0); 164 } 165 166 EXTERN_C_BEGIN 167 #undef __FUNCT__ 168 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 169 int MatConvert_SeqAIJ_Essl(Mat A,MatType type,Mat *newmat) { 170 Mat B=*newmat; 171 int ierr; 172 Mat_SeqAIJ_Essl *essl; 173 174 PetscFunctionBegin; 175 176 if (B != A) { 177 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 178 } 179 180 ierr = PetscNew(Mat_SeqAIJ_Essl,&essl);CHKERRQ(ierr); 181 essl->MatAssemblyEnd = A->ops->assemblyend; 182 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 183 essl->MatDestroy = A->ops->destroy; 184 B->spptr = (void *)essl; 185 186 B->ops->assemblyend = MatAssemblyEnd_SeqAIJ_Essl; 187 B->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 188 B->ops->destroy = MatDestroy_SeqAIJ_Essl; 189 190 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 191 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 192 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 193 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 194 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 195 *newmat = B; 196 PetscFunctionReturn(0); 197 } 198 EXTERN_C_END 199 200 EXTERN_C_BEGIN 201 #undef __FUNCT__ 202 #define __FUNCT__ "MatCreate_SeqAIJ_Essl" 203 int MatCreate_SeqAIJ_Essl(Mat A) { 204 int ierr; 205 206 PetscFunctionBegin; 207 ierr = MatSetType(A,MATSEQAIJ); 208 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 EXTERN_C_END 212