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,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,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 PetscFunctionBegin; 212 ierr = (*A->ops->duplicate)(A,op,M);CHKERRQ(ierr); 213 ierr = MatConvert_SeqAIJ_Essl(*M,MATESSL,M);CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 217 /*MC 218 MATESSL = "essl" - A matrix type providing direct solvers (LU) for sequential matrices 219 via the external package ESSL. 220 221 If ESSL is installed (see the manual for 222 instructions on how to declare the existence of external packages), 223 a matrix type can be constructed which invokes ESSL solvers. 224 After calling MatCreate(...,A), simply call MatSetType(A,MATESSL). 225 This matrix type is only supported for double precision real. 226 227 This matrix inherits from MATSEQAIJ. As a result, MatSeqAIJSetPreallocation is 228 supported for this matrix type. One can also call MatConvert for an inplace conversion to or from 229 the MATSEQAIJ type without data copy. 230 231 Options Database Keys: 232 . -mat_type essl - sets the matrix type to "essl" during a call to MatSetFromOptions() 233 234 Level: beginner 235 236 .seealso: PCLU 237 M*/ 238 239 EXTERN_C_BEGIN 240 #undef __FUNCT__ 241 #define __FUNCT__ "MatCreate_Essl" 242 int MatCreate_Essl(Mat A) { 243 int ierr; 244 245 PetscFunctionBegin; 246 /* Change type name before calling MatSetType to force proper construction of SeqAIJ and Essl types */ 247 ierr = PetscObjectChangeTypeName((PetscObject)A,MATESSL);CHKERRQ(ierr); 248 ierr = MatSetType(A,MATSEQAIJ); 249 ierr = MatConvert_SeqAIJ_Essl(A,MATESSL,&A);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 EXTERN_C_END 253