text: /* File: PCG.h */
  #ifndef SOLVERS_H
  #define SOLVERS_H
  
  /* Append underscore to fortran modual name */
  #ifdef _UF
    #define resprint   resprint_
  #endif
  #ifndef _UF
    #define resprint   RESPRINT
  #endif
  
  #include "r_vector.h"
  
  /* 
   *  Generic PCG:
   *
   *    Iterative method for approximating the solution to Ax=b  
   *    where x is an r_vector structure representing the initial guess
   *    and b is an r_vector structure representing the right-hand side.
   *    Upon completion, the x vector contains the approximate solution
   *    computed by the PCG method.
   *
   *    The matrix A and the preconditioner B are passed in 
   *    as pointers to generic operators.
   *
   *    Specifications:
   *
   *     IITER      Maximum number of iterations allowed.
   *
   *     RCLOSE     Stopping criterion for PCG.
   *
   *     pf         Print flag: If pf is non-zero, then the 
   *                l2-norm of the preconditioned residual, 
   *                the residual and the max-norm of the residual
   *                are printed.  The average reduction is also
   *                printed at the end of the iterations.
   *
   *     IOUT       FORTRAN unit number.
   *
   *     ow         If ow is non-zero, then the right-hand side
   *                is overwritten with the residual.  This
   *                saves memory.  For non-zero ow a zero 
   *                initial guess is used.
   *
   *  The iteration parameters are set using PCG_set.  For example:
   *
   *    PCG_set(&PCG,50,1.0e-10,1,IOUT);
   *
   *  will set the maximum number of iteration to 50,
   *  the stopping criterion will 1.0e-10 and the reduction
   *  history is printed to IOUT.
   *
   *  The PCG_eval call returns the number of iterations:
   *
   *    iter=PCG_eval(&x,&b,&PCG);
   *
   *  The final l2-norm of the residual squared 
   *  is stored in the PCG structure, PCG.BIGR
  */
  
  typedef struct PCG_operator
  {
    /* Work vectors. */
    r_vector r,p,z;
  
    /* Iteration parameters. */
    int IITER;
    double RCLOSE;
    int pf;
    int IOUT;
  
    /* 
      Pointers to generic operators.
    */
    GEN_operator *A_ptr,*B_ptr;
  
    int ow; /* Overwrite flag. */
  
    double BIGR; /* Final Residuals */
  }PCG_operator; 
  
  int PCG_allocate(GEN_operator* PCG_GEN_ptr, int ow, r_data* rdp);
  void PCG_free(void* A_ptr);
  int PCG_assemble(GEN_operator* PCG_GEN_ptr,
                   GEN_operator* A_GEN_ptr,
                   GEN_operator* B_GEN_ptr);
  int PCG_set(GEN_operator* PCG_GEN_ptr, int IITER,
              double RCLOSE, int pf, int IOUT);
  
  int PCG_eval(r_vector* d2_ptr, r_vector* d1_ptr, void* A_ptr);
  
  
  /*
   *  Generic multigrid algorithm:
   *
   *  MG_r_data contains a list of r_data structures and number of levels.
   *
   *  An mg_vector contains a list of r_vector's.  
   *  The r_vectors are allocated for a range of levels
   *  i0,..,i1.
   *
   *  An MG_GEN_operator contains a list of GEN_operator's.
   *  A list of GEN_operator's is used to define the multi-level
   *  smoothers, stiffness matrices, or prolongation/restriction
   *  pairs.
   *
   *  An MG_matrix contains mg_vector's, pointers to MG_GEN_operator's,
   *  and iteration parameters.  The MG_GEN_operators are defined
   *  outside the contex of the generic MG operator and depend on the
   *  multigrid method being used shuch as cell-centered finite-difference.
   *
   *  The generic MG structures allow for recursive definitions
   *  such as alternating plane smoothers which are 2-D multigrid
   *  solvers.
   */
  
  typedef struct MG_r_data
  {
    r_data* rd_list; /* Array of r_data objects */
    int levels;      /* Number of levels */
  }MG_r_data;
  
  int MG_r_data_allocate(MG_r_data* mgdp, int levels);
  void MG_r_data_free(MG_r_data* mgdp);
  
  /* Multilevel r_vector */
  typedef struct mg_vector
  {
    MG_r_data* mgrdp; /* Pointer to MG_r_data object */
    int i0,i1;        /* Levels i0 through i1 */
    r_vector* r_list; /* Array of r_vector objects */
  }mg_vector;
  
  int mg_vector_allocate(mg_vector* mgp, int i0, int i1, MG_r_data* mgrdp);
  void mg_vector_free(mg_vector* mgp);
  
  /* Multilevel GEN_operator */
  typedef struct MG_GEN_operator
  {
    MG_r_data* mgrdp;       /* Pointer to MG_r_data object */
    GEN_operator* GEN_list; /* Array of GEN_operator's */
  }MG_GEN_operator;
  
  int MG_GEN_allocate(MG_GEN_operator* MG_GEN_ptr, MG_r_data* mgrdp);
  void MG_GEN_free(MG_GEN_operator* MG_GEN_ptr);
  
  typedef struct MG_operator
  {
    MG_r_data* mgrdp;    /* Pointer to MG_r_data object */
  
    /* Overwrite Flag */
    int ow; /* ow=0-> allocate multilevel work vector */
  
    /* Lists of mg vectors */
    mg_vector mgp; /* Multilevel solution vector */
    mg_vector mgb; /* Multilevel right-hand side */
    mg_vector mgr; /* Multilevel residual vector */
    mg_vector mgw; /* Multilevel work vector */
  
    MG_GEN_operator* A_ptr;   /* Multilevel coefficient matrix */
    MG_GEN_operator* B_ptr;   /* Multilevel Smoother */
    MG_GEN_operator* P_ptr;   /* Multilevel prolongation operator */
    MG_GEN_operator* R_ptr;   /* Multilevel restriction operator */
  
    int mu0;        /* Number of smoothing iterations. */
    int mu1;        /* Number of multigrid cycles. */
    int nu;         /* V-Cycle(nu=1), W-Cycle(nu=2) */
  
  }MG_operator;
  
  int MG_allocate(GEN_operator* MG_ptr, int ow, MG_r_data* mgrdp);
  
  void MG_free(void* MG_free);
  
  int MG_assemble(GEN_operator* MG_GEN_ptr,
                  MG_GEN_operator* A_ptr,
                  MG_GEN_operator* B_ptr,
                  MG_GEN_operator* P_ptr,
                  MG_GEN_operator* R_ptr);
  
  void MG_set(GEN_operator* MG_ptr, int mu0, int mu1, int nu);
             
  int MG_eval(r_vector* p2_ptr, r_vector* p1_ptr, void* A_ptr);
  
  #endif
  