Subversion Repositories slepc-dev

/*

  SLEPc eigensolver: "davidson"

  Method: General Davidson Method

  References:

    - Ernest R. Davidson. Super-matrix methods. Computer Physics Communications,

      53:49–60, May 1989.

  TODO:

    - If the problem is symmetric, H=V*A*V is symmetric too, we can do some ops:

        *) Use a LAPACK eigensolver that only modify the upper triangular part

           of H, so we can save H in the lower one, and in a vector for the

           diagonal. That expends only 2*n^2 + n of memory.

    - In the interface, support static preconditioners of PETSc (PC) and others

      with shift parameter.

    - (DONE) Start with a krylov subspace, the H matrix will be used as the starting

      projected problem matrix, and A*V <- V*H + f*e_m^T.

    - Implements with BLAS3 prods Z <- V*U

*/

/*

   Dashboard struct: contains the methods that will be employed and the tunning

   options.

*/

#include "petsc.h"

#include "private/epsimpl.h"

typedef struct _dvdFunctionList {

  void *f;

  void *d;

  struct _dvdFunctionList *next;

} dvdFunctionList;

typedef struct _dvdDashboard {

  /**** Function steps ****/

  /* Initialize V */

  PetscInt (*initV)(struct _dvdDashboard*);

  void *initV_data;

  /* Find the approximate eigenpairs from V */

  PetscInt (*calcPairs)(struct _dvdDashboard*);

  void *calcPairs_data;

  /* Test for convergence */

  PetscTruth (*testConv)(struct _dvdDashboard*, PetscScalar eigvr,

                         PetscScalar eigvi, PetscScalar res, PetscReal *error);

  void *testConv_data;

  /* Number of converged eigenpairs */

  PetscInt nconv;

  /* Improve the selected eigenpairs */

  PetscInt (*improveX)(struct _dvdDashboard*);

  void *improveX_data;

  /* Check for restarting */

  PetscTruth (*isRestarting)(struct _dvdDashboard*);

  void *isRestarting_data;

  /* Perform restarting */

  PetscInt (*restartV)(struct _dvdDashboard*);

  void *restartV_data;

  /* Update V */

  PetscInt (*updateV)(struct _dvdDashboard*);

  void *updateV_data;

  /**** Problem specification ****/

  Mat A, B;         /* Problem matrices */

  PetscInt nev;     /* number of eigenpairs */

  EPSWhich which;   /* spectrum selection */

  PetscTruth

    withTarget;     /* if there is a target */

  PetscScalar

    target;         /* target value */

  PetscReal tol;    /* tolerance */

  /**** Subspaces specification ****/

  Vec *V,           /* searching subspace */

    *cX,            /* converged eigenvectors */

    *X,             /* best eigenvectors in V */

    *oldX,          /* eigenvectors from the previous iteration */

    *D,             /* improved X vectors */

    *W;             /* A*V */

  PetscInt size_V,  /* size of V */

    size_X,         /* size of X */

    size_oldX,      /* size of oldX */

    size_D,         /* size of D */

    size_cX,        /* size of cX */

    size_W,         /* size of W */

    max_size_V,     /* max size of V */

    max_size_X,     /* max size of X */

    max_size_oldX,  /* max size of oldX */

    max_size_D,     /* max size of D */

    max_size_W,     /* max size of W */

    size_uW;        /* size of updated W from V */

  /**** Auxiliary space ****/

  Vec *auxV;        /* auxiliary vectors */

  PetscScalar

    *auxS;          /* auxiliary scalars */

  PetscInt

    size_auxV,      /* max size of auxV */

    size_auxS;      /* max size of auxS */

  /**** Eigenvalues and errors ****/

  PetscScalar

    *ceigr, *ceigi, /* converged eigenvalues */

    *eigr, *eigi;   /* current eigenvalues */

  PetscReal

    *errest;        /* error eigenpairs */

  /**** Shared function and variables ****/

  PetscInt (*e_Vchanged)(struct _dvdDashboard*, PetscInt s_imm,

                         PetscInt e_imm, PetscInt s_new, PetscInt e_new);

  void *e_Vchanged_data;

  PetscInt (*calcpairs_residual)(struct _dvdDashboard*, PetscInt s, PetscInt e,

                                 Vec *R);

  void *calcpairs_residual_data;

  PetscInt (*e_newIteration)(struct _dvdDashboard*);

  void *e_newIteration_data;

  dvdFunctionList

    *destroyList;   /* destructor list */

  PetscInt

    n_calcPairs,    /* says to calcPairs how many ones */

    n_oldX,         /* says to calcParis how many old eigenvalues computes */

    n_convPairs;    /* says to updateV how many eigenvectors in X converged */

  PetscScalar *H;   /* Projected problem */

  PetscInt ldH,     /* leading dimension of H */

    size_H;         /* columns in H */

} dvdDashboard;

#define DVD_FL_ADD(list, fun) { \

  dvdFunctionList *fl=(list); \

  PetscErrorCode ierr; \

  ierr = PetscMalloc(sizeof(dvdFunctionList), &(list)); CHKERRQ(ierr); \

  (list)->f = (fun); \

  (list)->next = fl; }

#define DVD_FL_CALL1(list, t, arg0) { \

  dvdFunctionList *fl; \

  for(fl=(list); fl; fl=fl->next) (*(t)fl->f)((arg0)); }

#define DVD_FL_DEL(list) { \

  dvdFunctionList *fl=(list), *oldfl; \

  PetscErrorCode ierr; \

  while(fl) { \

    oldfl = fl; fl = fl->next; ierr = PetscFree(oldfl); CHKERRQ(ierr); }}

/*

  The blackboard configuration structure: saves information about the memory

  and other requirements

*/

typedef struct {

  PetscInt max_size_V,  /* max size of V */

    max_size_X,         /* max size of X */

    max_size_oldX,      /* max size of oldX */

    max_size_auxV,      /* max size of auxiliary vecs */

    max_size_auxS,      /* max size of auxiliary scalars */

    own_vecs,           /* number of global vecs */

    own_scalars;        /* number of local scalars */

  Vec *free_vecs;       /* free global vectors */

  PetscScalar

    *free_scalars;      /* free scalars */

  PetscInt state;       /* method states:

                            0: preconfiguring

                            1: configuring

                            2: running

                        */

} dvdBlackboard;

#define DVD_STATE_PRECONF 0

#define DVD_STATE_CONF 1

#define DVD_STATE_RUN 2

/* Shared types */

typedef void* dvdPrecondData;

typedef PetscErrorCode (*dvdPrecond)(dvdPrecondData, PetscInt i, Vec x, Vec Px);

typedef PetscInt (*dvdDestructor)(dvdDashboard*);

typedef PetscInt (*e_Vchanged_type)(dvdDashboard*, PetscInt s_imm,

                         PetscInt e_imm, PetscInt s_new, PetscInt e_new);

typedef PetscTruth (*isRestarting_type)(dvdDashboard*);

typedef PetscInt (*e_newIteration_type)(dvdDashboard*);

/* Routines for initV step */

PetscInt dvd_initV_classic(dvdDashboard *d, dvdBlackboard *b, PetscInt k);

PetscInt dvd_initV_krylov(dvdDashboard *d, dvdBlackboard *b, PetscInt k,

                          Mat OP, IP ip);

/* Routines for calcPairs step */

PetscInt dvd_calcpairs_rr(dvdDashboard *d, dvdBlackboard *b);

/* Routines for improveX step */

PetscInt dvd_improvex_gdbasic(dvdDashboard *d, dvdBlackboard *b, dvdPrecond P,

                              dvdPrecondData Pd);

/* Routines for testConv step */

PetscInt dvd_testconv_basic(dvdDashboard *d, dvdBlackboard *b);

/* Routines for management of V */

PetscInt dvd_managementV_basic(dvdDashboard *d, dvdBlackboard *b,

                               PetscInt bs, PetscInt max_size_V,

                               PetscInt restart_size_X, PetscInt plusk);

/* Some utilities */

PetscInt dvd_orthV(dvdDashboard *d, dvdBlackboard *b, IP ip);

PetscErrorCode dvd_static_precond_PC(dvdDashboard *d, PC pc, dvdPrecond *pf,

                                     dvdPrecondData *pd);

PetscErrorCode dvd_jacobi_precond(dvdDashboard *d, Vec diagA, dvdPrecond *pf,

                                  dvdPrecondData *pd);

PetscInt dvd_isrestarting_whenfinish(dvdDashboard *d, dvdBlackboard *b);

/* Methods */

PetscErrorCode dvd_schm_basic_preconf(dvdDashboard *d, dvdBlackboard *b,

  PetscInt max_size_V, PetscInt min_size_V, PetscInt bs, PetscInt size_initV,

  PetscInt plusk);

PetscErrorCode dvd_schm_basic_conf(dvdDashboard *d, dvdBlackboard *b,

  PetscInt max_size_V, PetscInt min_size_V, PetscInt bs, PetscInt size_initV,

  PetscInt plusk, PC pc, Vec diagA, IP ip);

/* BLAS routines */

PetscErrorCode SlepcDenseMatProd(PetscScalar *C, PetscInt _ldC, PetscScalar b,

  PetscScalar a,

  const PetscScalar *A, PetscInt _ldA, PetscInt rA, PetscInt cA, PetscTruth At,

  const PetscScalar *B, PetscInt _ldB, PetscInt rB, PetscInt cB, PetscTruth Bt);

PetscErrorCode SlepcUpdateVectorsZ(Vec *Y, Vec *X, PetscInt cX,

  const PetscScalar *M, PetscInt ldM, PetscInt rM, PetscInt cM);