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 (*MatDestroy)(Mat); 23 } Mat_SeqAIJ_Essl; 24 25 #undef __FUNCT__ 26 #define __FUNCT__ "MatDestroy_SeqAIJ_Essl" 27 int MatDestroy_SeqAIJ_Essl(Mat A) 28 { 29 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 30 int ierr,(*destroy)(Mat); 31 32 PetscFunctionBegin; 33 /* free the Essl datastructures */ 34 destroy = essl->MatDestroy; 35 ierr = PetscFree(essl->a);CHKERRQ(ierr); 36 ierr = (*destroy)(A);CHKERRQ(ierr); 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "MatSolve_SeqAIJ_Essl" 42 int MatSolve_SeqAIJ_Essl(Mat A,Vec b,Vec x) 43 { 44 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)A->spptr; 45 PetscScalar *xx; 46 int ierr,m,zero = 0; 47 48 PetscFunctionBegin; 49 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 50 ierr = VecCopy(b,x);CHKERRQ(ierr); 51 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 52 dgss(&zero,&A->n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 53 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 54 PetscFunctionReturn(0); 55 } 56 57 #undef __FUNCT__ 58 #define __FUNCT__ "MatLUFactorNumeric_SeqAIJ_Essl" 59 int MatLUFactorNumeric_SeqAIJ_Essl(Mat A,Mat *F) 60 { 61 Mat_SeqAIJ *a = (Mat_SeqAIJ*)(*F)->data; 62 Mat_SeqAIJ *aa = (Mat_SeqAIJ*)(A)->data; 63 Mat_SeqAIJ_Essl *essl = (Mat_SeqAIJ_Essl*)(*F)->spptr; 64 int i,ierr,one = 1; 65 66 PetscFunctionBegin; 67 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 68 for (i=0; i<A->m+1; i++) essl->ia[i] = aa->i[i] + 1; 69 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 70 71 ierr = PetscMemcpy(essl->a,aa->a,(aa->nz)*sizeof(PetscScalar));CHKERRQ(ierr); 72 73 /* set Essl options */ 74 essl->iparm[0] = 1; 75 essl->iparm[1] = 5; 76 essl->iparm[2] = 1; 77 essl->iparm[3] = 0; 78 essl->rparm[0] = 1.e-12; 79 essl->rparm[1] = A->lupivotthreshold; 80 81 dgsf(&one,&A->m,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm, 82 essl->rparm,essl->oparm,essl->aux,&essl->naux); 83 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatLUFactorSymbolic_SeqAIJ_Essl" 89 int MatLUFactorSymbolic_SeqAIJ_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) 90 { 91 Mat B; 92 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b; 93 int ierr,*ridx,*cidx,i,len; 94 Mat_SeqAIJ_Essl *essl; 95 PetscReal f = 1.0; 96 97 PetscFunctionBegin; 98 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 99 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 100 ierr = MatSetType(B,MATESSL);CHKERRQ(ierr); 101 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 102 103 B->ops->solve = MatSolve_SeqAIJ_Essl; 104 B->ops->lufactornumeric = MatLUFactorNumeric_SeqAIJ_Essl; 105 B->factor = FACTOR_LU; 106 107 essl = (Mat_SeqAIJ_Essl *)(B->spptr); 108 109 /* allocate the work arrays required by ESSL */ 110 f = info->fill; 111 essl->nz = a->nz; 112 essl->lna = (int)a->nz*f; 113 essl->naux = 100 + 10*A->m; 114 115 /* since malloc is slow on IBM we try a single malloc */ 116 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 117 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 118 essl->aux = essl->a + essl->lna; 119 essl->ia = (int*)(essl->aux + essl->naux); 120 essl->ja = essl->ia + essl->lna; 121 122 PetscLogObjectMemory(B,len+sizeof(Mat_SeqAIJ_Essl)); 123 *F = B; 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "MatUseEssl_SeqAIJ" 129 int MatUseEssl_SeqAIJ(Mat A) 130 { 131 PetscFunctionBegin; 132 A->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqAIJ_Essl; 133 PetscLogInfo(0,"Using ESSL for SeqAIJ LU factorization and solves"); 134 PetscFunctionReturn(0); 135 } 136 137 EXTERN_C_BEGIN 138 #undef __FUNCT__ 139 #define __FUNCT__ "MatCreate_SeqAIJ_Essl" 140 int MatCreate_SeqAIJ_Essl(Mat A) { 141 int ierr; 142 Mat_SeqAIJ_Essl *essl; 143 144 PetscFunctionBegin; 145 ierr = MatSetType(A,MATSEQAIJ); 146 ierr = MatUseEssl_SeqAIJ(A); 147 148 ierr = PetscNew(Mat_SeqAIJ_Essl,&essl);CHKERRQ(ierr); 149 essl->MatDestroy = A->ops->destroy; 150 A->spptr = (void *)essl; 151 152 A->ops->destroy = MatDestroy_SeqAIJ_Essl; 153 PetscFunctionReturn(0); 154 } 155