2: /*
  3:         Provides an interface to the IBM RS6000 Essl sparse solver

  5: */
  6: #include <../src/mat/impls/aij/seq/aij.h>

  8: /* #include <essl.h> This doesn't work!  */

 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*);

 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;

 25:   PetscBool   CleanUpESSL;
 26: } Mat_Essl;

 28: PetscErrorCode MatDestroy_Essl(Mat A)
 29: {
 30:   Mat_Essl       *essl=(Mat_Essl*)A->data;

 32:   if (essl->CleanUpESSL) {
 33:     PetscFree4(essl->a,essl->aux,essl->ia,essl->ja);
 34:   }
 35:   PetscFree(A->data);
 36:   return 0;
 37: }

 39: PetscErrorCode MatSolve_Essl(Mat A,Vec b,Vec x)
 40: {
 41:   Mat_Essl       *essl = (Mat_Essl*)A->data;
 42:   PetscScalar    *xx;
 43:   int            nessl,zero = 0;

 45:   PetscBLASIntCast(A->cmap->n,&nessl);
 46:   VecCopy(b,x);
 47:   VecGetArray(x,&xx);
 48:   dgss(&zero,&nessl,essl->a,essl->ia,essl->ja,&essl->lna,xx,essl->aux,&essl->naux);
 49:   VecRestoreArray(x,&xx);
 50:   return 0;
 51: }

 53: PetscErrorCode MatLUFactorNumeric_Essl(Mat F,Mat A,const MatFactorInfo *info)
 54: {
 55:   Mat_SeqAIJ     *aa  =(Mat_SeqAIJ*)(A)->data;
 56:   Mat_Essl       *essl=(Mat_Essl*)(F)->data;
 57:   int            nessl,i,one = 1;

 59:   PetscBLASIntCast(A->rmap->n,&nessl);
 60:   /* copy matrix data into silly ESSL data structure (1-based Frotran style) */
 61:   for (i=0; i<A->rmap->n+1; i++) essl->ia[i] = aa->i[i] + 1;
 62:   for (i=0; i<aa->nz; i++) essl->ja[i] = aa->j[i] + 1;

 64:   PetscArraycpy(essl->a,aa->a,aa->nz);

 66:   /* set Essl options */
 67:   essl->iparm[0] = 1;
 68:   essl->iparm[1] = 5;
 69:   essl->iparm[2] = 1;
 70:   essl->iparm[3] = 0;
 71:   essl->rparm[0] = 1.e-12;
 72:   essl->rparm[1] = 1.0;

 74:   PetscOptionsGetReal(NULL,((PetscObject)A)->prefix,"-matessl_lu_threshold",&essl->rparm[1],NULL);

 76:   dgsf(&one,&nessl,&essl->nz,essl->a,essl->ia,essl->ja,&essl->lna,essl->iparm,essl->rparm,essl->oparm,essl->aux,&essl->naux);

 78:   F->ops->solve     = MatSolve_Essl;
 79:   (F)->assembled    = PETSC_TRUE;
 80:   (F)->preallocated = PETSC_TRUE;
 81:   return 0;
 82: }

 84: PetscErrorCode MatLUFactorSymbolic_Essl(Mat B,Mat A,IS r,IS c,const MatFactorInfo *info)
 85: {
 86:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
 87:   Mat_Essl       *essl;
 88:   PetscReal      f = 1.0;

 90:   essl = (Mat_Essl*)(B->data);

 92:   /* allocate the work arrays required by ESSL */
 93:   f    = info->fill;
 94:   PetscBLASIntCast(a->nz,&essl->nz);
 95:   PetscBLASIntCast((PetscInt)(a->nz*f),&essl->lna);
 96:   PetscBLASIntCast(100 + 10*A->rmap->n,&essl->naux);

 98:   /* since malloc is slow on IBM we try a single malloc */
 99:   PetscMalloc4(essl->lna,&essl->a,essl->naux,&essl->aux,essl->lna,&essl->ia,essl->lna,&essl->ja);

101:   essl->CleanUpESSL = PETSC_TRUE;

103:   PetscLogObjectMemory((PetscObject)B,essl->lna*(2*sizeof(int)+sizeof(PetscScalar)) + essl->naux*sizeof(PetscScalar));

105:   B->ops->lufactornumeric = MatLUFactorNumeric_Essl;
106:   return 0;
107: }

109: PetscErrorCode MatFactorGetSolverType_essl(Mat A,MatSolverType *type)
110: {
111:   *type = MATSOLVERESSL;
112:   return 0;
113: }

115: /*MC
116:   MATSOLVERESSL - "essl" - Provides direct solvers (LU) for sequential matrices
117:                               via the external package ESSL.

119:   If ESSL is installed (see the manual for
120:   instructions on how to declare the existence of external packages),

122:   Works with MATSEQAIJ matrices

124:    Level: beginner

126: .seealso: PCLU, PCFactorSetMatSolverType(), MatSolverType
127: M*/

129: PETSC_EXTERN PetscErrorCode MatGetFactor_seqaij_essl(Mat A,MatFactorType ftype,Mat *F)
130: {
131:   Mat            B;
132:   Mat_Essl       *essl;

135:   MatCreate(PetscObjectComm((PetscObject)A),&B);
136:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,A->rmap->n,A->cmap->n);
137:   PetscStrallocpy("essl",&((PetscObject)B)->type_name);
138:   MatSetUp(B);

140:   PetscNewLog(B,&essl);

142:   B->data                  = essl;
143:   B->ops->lufactorsymbolic = MatLUFactorSymbolic_Essl;
144:   B->ops->destroy          = MatDestroy_Essl;
145:   B->ops->getinfo          = MatGetInfo_External;

147:   PetscObjectComposeFunction((PetscObject)B,"MatFactorGetSolverType_C",MatFactorGetSolverType_essl);

149:   B->factortype = MAT_FACTOR_LU;
150:   PetscStrallocpy(MATORDERINGEXTERNAL,(char**)&B->preferredordering[MAT_FACTOR_LU]);
151:   PetscFree(B->solvertype);
152:   PetscStrallocpy(MATSOLVERESSL,&B->solvertype);

154:   *F            = B;
155:   return 0;
156: }

158: PETSC_EXTERN PetscErrorCode MatSolverTypeRegister_Essl(void)
159: {
160:   MatSolverTypeRegister(MATSOLVERESSL,MATSEQAIJ,          MAT_FACTOR_LU,MatGetFactor_seqaij_essl);
161:   return 0;
162: }