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 (*MatDuplicate)(Mat,MatDuplicateOption,Mat*); 23 int (*MatAssemblyEnd)(Mat,MatAssemblyType); 24 int (*MatLUFactorSymbolic)(Mat,IS,IS,MatFactorInfo*,Mat*); 25 int (*MatDestroy)(Mat); 26 PetscTruth CleanUpESSL; 27 } Mat_Essl; 28 29 EXTERN int MatDuplicate_Essl(Mat,MatDuplicateOption,Mat*); 30 31 EXTERN_C_BEGIN 32 #undef __FUNCT__ 33 #define __FUNCT__ "MatConvert_Essl_SeqAIJ" 34 int MatConvert_Essl_SeqAIJ(Mat A,const MatType type,Mat *newmat) { 35 int ierr; 36 Mat B=*newmat; 37 Mat_Essl *essl=(Mat_Essl*)A->spptr; 38 39 PetscFunctionBegin; 40 if (B != A) { 41 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); 42 } 43 B->ops->duplicate = essl->MatDuplicate; 44 B->ops->assemblyend = essl->MatAssemblyEnd; 45 B->ops->lufactorsymbolic = essl->MatLUFactorSymbolic; 46 B->ops->destroy = essl->MatDestroy; 47 48 /* free the Essl datastructures */ 49 ierr = PetscFree(essl);CHKERRQ(ierr); 50 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 51 *newmat = B; 52 PetscFunctionReturn(0); 53 } 54 EXTERN_C_END 55 56 #undef __FUNCT__ 57 #define __FUNCT__ "MatDestroy_Essl" 58 int MatDestroy_Essl(Mat A) { 59 int ierr; 60 Mat_Essl *essl=(Mat_Essl*)A->spptr; 61 62 PetscFunctionBegin; 63 if (essl->CleanUpESSL) { 64 ierr = PetscFree(essl->a);CHKERRQ(ierr); 65 } 66 ierr = MatConvert_Essl_SeqAIJ(A,MATSEQAIJ,&A); 67 ierr = (*A->ops->destroy)(A);CHKERRQ(ierr); 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "MatSolve_Essl" 73 int MatSolve_Essl(Mat A,Vec b,Vec x) { 74 Mat_Essl *essl = (Mat_Essl*)A->spptr; 75 PetscScalar *xx; 76 int ierr,m,zero = 0; 77 78 PetscFunctionBegin; 79 ierr = VecGetLocalSize(b,&m);CHKERRQ(ierr); 80 ierr = VecCopy(b,x);CHKERRQ(ierr); 81 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 82 dgss(&zero,&A->n,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux); 83 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatLUFactorNumeric_Essl" 89 int MatLUFactorNumeric_Essl(Mat A,Mat *F) { 90 Mat_SeqAIJ *aa=(Mat_SeqAIJ*)(A)->data; 91 Mat_Essl *essl=(Mat_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_Essl" 117 int MatLUFactorSymbolic_Essl(Mat A,IS r,IS c,MatFactorInfo *info,Mat *F) { 118 Mat B; 119 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 120 int ierr,len; 121 Mat_Essl *essl; 122 PetscReal f = 1.0; 123 124 PetscFunctionBegin; 125 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 126 ierr = MatCreate(A->comm,PETSC_DECIDE,PETSC_DECIDE,A->m,A->n,&B);CHKERRQ(ierr); 127 ierr = MatSetType(B,MATESSL);CHKERRQ(ierr); 128 ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr); 129 130 B->ops->solve = MatSolve_Essl; 131 B->ops->lufactornumeric = MatLUFactorNumeric_Essl; 132 B->factor = FACTOR_LU; 133 134 essl = (Mat_Essl *)(B->spptr); 135 136 /* allocate the work arrays required by ESSL */ 137 f = info->fill; 138 essl->nz = a->nz; 139 essl->lna = (int)a->nz*f; 140 essl->naux = 100 + 10*A->m; 141 142 /* since malloc is slow on IBM we try a single malloc */ 143 len = essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar); 144 ierr = PetscMalloc(len,&essl->a);CHKERRQ(ierr); 145 essl->aux = essl->a + essl->lna; 146 essl->ia = (int*)(essl->aux + essl->naux); 147 essl->ja = essl->ia + essl->lna; 148 essl->CleanUpESSL = PETSC_TRUE; 149 150 PetscLogObjectMemory(B,len+sizeof(Mat_Essl)); 151 *F = B; 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatAssemblyEnd_Essl" 157 int MatAssemblyEnd_Essl(Mat A,MatAssemblyType mode) { 158 int ierr; 159 Mat_Essl *essl=(Mat_Essl*)(A->spptr); 160 161 PetscFunctionBegin; 162 ierr = (*essl->MatAssemblyEnd)(A,mode);CHKERRQ(ierr); 163 164 essl->MatLUFactorSymbolic = A->ops->lufactorsymbolic; 165 A->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 166 PetscLogInfo(0,"Using ESSL for LU factorization and solves"); 167 PetscFunctionReturn(0); 168 } 169 170 EXTERN_C_BEGIN 171 #undef __FUNCT__ 172 #define __FUNCT__ "MatConvert_SeqAIJ_Essl" 173 int MatConvert_SeqAIJ_Essl(Mat A,const MatType type,Mat *newmat) { 174 Mat B=*newmat; 175 int ierr; 176 Mat_Essl *essl; 177 178 PetscFunctionBegin; 179 180 if (B != A) { 181 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 182 } 183 184 ierr = PetscNew(Mat_Essl,&essl);CHKERRQ(ierr); 185 essl->MatDuplicate = A->ops->duplicate; 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->duplicate = MatDuplicate_Essl; 193 B->ops->assemblyend = MatAssemblyEnd_Essl; 194 B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl; 195 B->ops->destroy = MatDestroy_Essl; 196 197 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_essl_C", 198 "MatConvert_SeqAIJ_Essl",MatConvert_SeqAIJ_Essl);CHKERRQ(ierr); 199 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_essl_seqaij_C", 200 "MatConvert_Essl_SeqAIJ",MatConvert_Essl_SeqAIJ);CHKERRQ(ierr); 201 ierr = PetscObjectChangeTypeName((PetscObject)B,type);CHKERRQ(ierr); 202 *newmat = B; 203 PetscFunctionReturn(0); 204 } 205 EXTERN_C_END 206 207 #undef __FUNCT__ 208 #define __FUNCT__ "MatDuplicate_Essl" 209 int MatDuplicate_Essl(Mat A, MatDuplicateOption op, Mat *M) { 210 int ierr; 211 Mat_Essl *lu = (Mat_Essl *)A->spptr; 212 213 PetscFunctionBegin; 214 ierr = (*lu->MatDuplicate)(A,op,M);CHKERRQ(ierr); 215 ierr = MatConvert_SeqAIJ_Essl(*M,MATESSL,M);CHKERRQ(ierr); 216 ierr = PetscMemcpy((*M)->spptr,lu,sizeof(Mat_Essl));CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 /*MC 221 MATESSL - MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 222 via the external package ESSL. 223 224 If ESSL is installed (see the manual for 225 instructions on how to declare the existence of external packages), 226 a matrix type can be constructed which invokes ESSL solvers. 227 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 228 This matrix type is only supported for double precision real. 229 230 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 231 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 232 the MATSEQAIJ type without data copy. 233 234 Options Database Keys: 235 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 236 237 Level: beginner 238 239 .seealso: PCLU 240 M*/ 241 242 EXTERN_C_BEGIN 243 #undef __FUNCT__ 244 #define __FUNCT__ "MatCreate_Essl" 245 int MatCreate_Essl(Mat A) { 246 int ierr; 247 248 PetscFunctionBegin; 249 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */ 250 ierr = PetscObjectChangeTypeName((PetscObject)A,MATESSL);CHKERRQ(ierr); 251 ierr = MatSetType(A,MATSEQAIJ); 252 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 EXTERN_C_END 256