Subversion Repositories slepc-dev

/*                      

       This file implements a wrapper to the PLANSO package

*/

#include "src/eps/impls/planso/plansop.h"

/* Nasty global variable to access EPS data from PLANop_ and PLANopm_ */

static EPS globaleps;

#undef __FUNCT__  

#define __FUNCT__ "EPSSetUp_PLANSO"

static int EPSSetUp_PLANSO(EPS eps)

{

  int        ierr, n;

  EPS_PLANSO *pl = (EPS_PLANSO *)eps->data;

  PetscFunctionBegin;

#if defined(PETSC_USE_COMPLEX)

  SETERRQ(PETSC_ERR_SUP,"Requested method is not available for complex problems");

#endif

  if (!eps->ishermitian)

    SETERRQ(PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");

  ierr = VecGetLocalSize(eps->vec_initial,&n); CHKERRQ(ierr);

  pl->lwork = 5*n+1+4*eps->ncv+PetscMax(n,eps->ncv+1);

  ierr = PetscMalloc(pl->lwork*sizeof(PetscReal),&pl->work);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSSetDefaults_PLANSO"

static int EPSSetDefaults_PLANSO(EPS eps)

{

  int         ierr, N;

  PetscFunctionBegin;

  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);

  if (eps->ncv) {

    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");

  }

  else eps->ncv = eps->nev;

  if (!eps->max_it) eps->max_it = PetscMax(100,N);

  if (!eps->tol) eps->tol = 1.e-7;

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PLANop_"

int  PLANop_(int *n,PetscReal *s, PetscReal *q, PetscReal *p)

{

  Vec    x,y;

  int    ierr;

  PetscFunctionBegin;

  ierr = VecCreateMPIWithArray(globaleps->comm,*n,PETSC_DECIDE,(PetscScalar*)q,&x);CHKERRQ(ierr);

  ierr = VecCreateMPIWithArray(globaleps->comm,*n,PETSC_DECIDE,(PetscScalar*)p,&y);CHKERRQ(ierr);

  ierr = STApply(globaleps->OP,x,y);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "PLANopm_"

int  PLANopm_(int *n,PetscReal *q, PetscReal *s)

{

  Vec    x,y;

  int    ierr;

  PetscFunctionBegin;

  ierr = VecCreateMPIWithArray(globaleps->comm,*n,PETSC_DECIDE,(PetscScalar*)q,&x);CHKERRQ(ierr);

  ierr = VecCreateMPIWithArray(globaleps->comm,*n,PETSC_DECIDE,(PetscScalar*)s,&y);CHKERRQ(ierr);

  ierr = STApplyB(globaleps->OP,x,y);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSSolve_PLANSO"

static int  EPSSolve_PLANSO(EPS eps,int *its)

{

  int        i, n, msglvl, lohi, ierr;

  PetscReal  condm;

  EPS_PLANSO *pl = (EPS_PLANSO *)eps->data;

  MPI_Fint    fcomm;

  PetscFunctionBegin;

  ierr = VecGetLocalSize(eps->vec_initial,&n); CHKERRQ(ierr);

  if (eps->which==EPS_LARGEST_MAGNITUDE) lohi = 1;

  else if (eps->which==EPS_SMALLEST_MAGNITUDE) lohi = -1;

  condm = 1.0;         /* estimated condition number: we have no information */

  msglvl = 0;

  globaleps = eps;

  fcomm = MPI_Comm_c2f(eps->comm);

  PLANdr2_( &n, &eps->ncv, &eps->nev, &lohi, &condm, &eps->tol, &eps->its, &eps->nconv,

            eps->eigr, eps->eigi, pl->work, &pl->lwork, &ierr, &msglvl, &fcomm );

  for (i=0;i<eps->nconv;i++) eps->eigi[i]=0.0;

  *its       = eps->its;

  eps->reason = EPS_CONVERGED_TOL;

  if (ierr!=0) { SETERRQ1(PETSC_ERR_LIB,"Error in PLANSO (code=%d)",ierr);}

  PetscFunctionReturn(0);

}

#undef __FUNCT__  

#define __FUNCT__ "EPSDestroy_PLANSO"

/*

  EPSDestroy_PLANSO - Destroys the context variable for PLANSO.

  Input Parameter:

. eps - the iterative context

*/

int EPSDestroy_PLANSO(EPS eps)

{

  EPS_PLANSO *pl = (EPS_PLANSO *)eps->data;

  int         ierr;

  PetscFunctionBegin;

  PetscValidHeaderSpecific(eps,EPS_COOKIE);

  if (pl->work)  { ierr = PetscFree(pl->work);CHKERRQ(ierr); }

  if (eps->data) { ierr = PetscFree(eps->data);CHKERRQ(ierr); }

  PetscFunctionReturn(0);

}

EXTERN_C_BEGIN

#undef __FUNCT__  

#define __FUNCT__ "EPSCreate_PLANSO"

int EPSCreate_PLANSO(EPS eps)

{

  EPS_PLANSO *planso;

  int        ierr;

  PetscFunctionBegin;

  ierr = PetscNew(EPS_PLANSO,&planso);CHKERRQ(ierr);

  PetscMemzero(planso,sizeof(EPS_PLANSO));

  PetscLogObjectMemory(eps,sizeof(EPS_PLANSO));

  eps->data                      = (void *) planso;

  eps->ops->setup                = EPSSetUp_PLANSO;

  eps->ops->setdefaults          = EPSSetDefaults_PLANSO;

  eps->ops->solve                = EPSSolve_PLANSO;

  eps->ops->destroy              = EPSDestroy_PLANSO;

  eps->ops->view                 = 0;

  PetscFunctionReturn(0);

}

EXTERN_C_END