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 8 /* #include <essl.h> This doesn't work! */ 9 10 PETSC_EXTERN void dgss(int*,int*,double*,int*,int*,int*,double*,double*,int*); 11 PETSC_EXTERN void dgsf(int*,int*,int*,double*,int*,int*,int*,int*,double*,double*,double*,int*); 12 13 typedef struct { 14 int n,nz; 15 PetscScalar *a; 16 int *ia; 17 int *ja; 18 int lna; 19 int iparm[5]; 20 PetscReal rparm[5]; 21 PetscReal oparm[5]; 22 PetscScalar *aux; 23 int naux; 24 25 PetscBool CleanUpESSL; 26 } Mat_Essl; 27 28 PetscErrorCode MatDestroy_Essl(Mat A) 29 { 30 PetscErrorCode ierr; 31 Mat_Essl *essl=(Mat_Essl*)A->data; 32 33 PetscFunctionBegin; 34 if (essl->CleanUpESSL) { 35 ierr = PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);CHKERRQ(ierr); 36 } 37 ierr = PetscFree(A->data);CHKERRQ(ierr); 38 PetscFunctionReturn(0); 39 } 40 41 PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x) 42 { 43 Mat_Essl *essl = (Mat_Essl*)A->data; 44 PetscScalar *xx; 45 PetscErrorCode ierr; 46 int nessl,zero = 0; 47 48 PetscFunctionBegin; 49 ierr = PetscBLASIntCast(A->cmap->n,&nessl);CHKERRQ(ierr); 50 ierr = VecCopy(b,x);CHKERRQ(ierr); 51 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 52 dgss(&zero,&nessl,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 PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info) 58 { 59 Mat_SeqAIJ *aa =(Mat_SeqAIJ*)(A)->data; 60 Mat_Essl *essl=(Mat_Essl*)(F)->data; 61 PetscErrorCode ierr; 62 int nessl,i,one = 1; 63 64 PetscFunctionBegin; 65 ierr = PetscBLASIntCast(A->rmap->n,&nessl);CHKERRQ(ierr); 66 /* copy matrix data into silly ESSL data structure (1-based Frotran style) */ 67 for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1; 68 for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1; 69 70 ierr = PetscArraycpy(essl->a,aa->a,aa->nz);CHKERRQ(ierr); 71 72 /* set Essl options */ 73 essl->iparm[0] = 1; 74 essl->iparm[1] = 5; 75 essl->iparm[2] = 1; 76 essl->iparm[3] = 0; 77 essl->rparm[0] = 1.e-12; 78 essl->rparm[1] = 1.0; 79 80 ierr = PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);CHKERRQ(ierr); 81 82 dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux); 83 84 F->ops->solve = MatSolve_Essl; 85 (F)->assembled = PETSC_TRUE; 86 (F)->preallocated = PETSC_TRUE; 87 PetscFunctionReturn(0); 88 } 89 90 PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info) 91 { 92 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 93 PetscErrorCode ierr; 94 Mat_Essl *essl; 95 PetscReal f = 1.0; 96 97 PetscFunctionBegin; 98 essl = (Mat_Essl*)(B->data); 99 100 /* allocate the work arrays required by ESSL */ 101 f = info->fill; 102 ierr = PetscBLASIntCast(a->nz,&essl->nz);CHKERRQ(ierr); 103 ierr = PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);CHKERRQ(ierr); 104 ierr = PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);CHKERRQ(ierr); 105 106 /* since malloc is slow on IBM we try a single malloc */ 107 ierr = PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);CHKERRQ(ierr); 108 109 essl->CleanUpESSL = PETSC_TRUE; 110 111 ierr = PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));CHKERRQ(ierr); 112 113 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 114 PetscFunctionReturn(0); 115 } 116 117 PetscErrorCode MatFactorGetSolverType_essl(Mat A,MatSolverType *type) 118 { 119 PetscFunctionBegin; 120 *type = MATSOLVERESSL; 121 PetscFunctionReturn(0); 122 } 123 124 /*MC 125 MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices 126 via the external package ESSL. 127 128 If ESSL is installed (see the manual for 129 instructions on how to declare the existence of external packages), 130 131 Works with MATSEQAIJ matrices 132 133 Level: beginner 134 135 .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType 136 M*/ 137 138 PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F) 139 { 140 Mat B; 141 PetscErrorCode ierr; 142 Mat_Essl *essl; 143 144 PetscFunctionBegin; 145 if (A->cmap->N != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix must be square"); 146 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 147 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 148 ierr = PetscStrallocpy("essl",&((PetscObject)B)->type_name);CHKERRQ(ierr); 149 ierr = MatSetUp(B);CHKERRQ(ierr); 150 151 ierr = PetscNewLog(B,&essl);CHKERRQ(ierr); 152 153 B->data = essl; 154 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 155 B->ops->destroy = MatDestroy_Essl; 156 B->ops->getinfo = MatGetInfo_External; 157 158 ierr = PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl);CHKERRQ(ierr); 159 160 B->factortype = MAT_FACTOR_LU; 161 ierr = PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU]);CHKERRQ(ierr); 162 ierr = PetscFree(B->solvertype);CHKERRQ(ierr); 163 ierr = PetscStrallocpy(MATSOLVERESSL,&B->solvertype);CHKERRQ(ierr); 164 165 *F = B; 166 PetscFunctionReturn(0); 167 } 168 169 PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void) 170 { 171 PetscErrorCode ierr; 172 PetscFunctionBegin; 173 ierr = MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ, MAT_FACTOR_LU,MatGetFactor_seqaij_essl);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176