#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: ml.h,v 1.6 2001/01/27 21:32:36 knepley Exp $";
#endif
/* 
   Private data structure for Multilevel preconditioner.
*/
#if !defined(__ML_H)
#define __ML_H

#include "gvec.h"
#ifdef PETSC_HAVE_PLAPACK
  #include "PLA.h"
#endif
#include "petscblaslapack.h"

enum {PC_ML_SetUpInit, PC_ML_SetUpConstrained, PC_ML_SetUpConstrainedBd, PC_ML_SetUpParallel,
      PC_ML_ReducePartitionMesh, PC_ML_ReducePartitionRowCol, PC_ML_ReduceFactor, PC_ML_ReduceBdGrad, PC_ML_ReduceBdGradExtract,
      PC_ML_ReduceBdGradRowPartLocalToGlobal, PC_ML_ReduceShrinkMesh, PC_ML_ApplySymmetricLeftParallel,
      PC_ML_ApplySymmetricRightParallel, PC_ML_QRFactorization, PC_ML_ApplyQR, PC_ML_CreateBdGrad, PC_ML_MAX_EVENTS};
extern int PC_MLEvents[PC_ML_MAX_EVENTS];
#define PC_MLLogEventBegin(e,o1,o2,o3,o4) PetscLogEventBegin(PC_MLEvents[e],o1,o2,o3,o4)
#define PC_MLLogEventEnd(e,o1,o2,o3,o4)   PetscLogEventEnd(PC_MLEvents[e],o1,o2,o3,o4)

enum {MESH_OFFSETS, MESH_ADJ, MESH_ELEM, NUM_MESH_DIV};
enum {PART_ROW_INT, PART_ROW_BD, PART_ROW_RES, NUM_PART_ROW_DIV};
enum {FACT_U, FACT_DINV, FACT_V, FACT_QR, FACT_TAU, NUM_FACT_DIV};

typedef struct {
  Grid             grid;                /* The problem grid */
  PetscTruth       useMath;             /* Use Mathematica functions instead */
  PetscViewer      mathViewer;          /* The link to Mathematica */
  PetscViewer      factorizationViewer; /* The visualization of the factorization process */

  /* Gradient variabels */
  int              gradOp;              /* The gradient operator */
  int              divOp;               /* The divergence operator */
  int              sField;              /* The shape function field */
  int              tField;              /* The test function field */
  VarOrdering      sOrder;              /* The global variable ordering for the shape function in B */
  LocalVarOrdering sLocOrder;           /* The local  variable ordering for the shape function in B */
  VarOrdering      tOrder;              /* The global variable ordering for the test  function in B */
  LocalVarOrdering tLocOrder;           /* The local  variable ordering for the test  function in B */
  PetscScalar      gradAlpha;           /* The scalar multiplier for the gradient operator */
  GMat             B;                   /* The gradient matrix */
  Mat              locB;                /* The local gradient matrix */
  GVec             reduceSol;           /* The solution in the range of B */
  GVec             nonlinearIterate;    /* The current iterate in a Newton solve */
  GVec             diag;                /* The diagonal of the jacobian (approximately) */

  /* Mesh variables */
  int              numMeshes;           /* The number of meshes corresponding to the factorization hierarchy */
  int              partSize;            /* The target size of column partitions */
  int             *numElements;         /* The number of elements in each mesh */
  int             *numVertices;         /* The number of vertices in each mesh */
  double          *vertices;            /* The coordinates for each vertex */
  int              numEdges;            /* The number of edges in the original mesh */
  int           ***meshes;              /* The mesh hierarchy in CSR format */
  
  /* Numbering variables */
  int             *numPartitions;       /* The number of partitions for each level */
  int            **numPartitionCols;    /* The number of columns in each partition and the final QR */
  int           ***colPartition;        /* The column or node partition */
  int            **numPartitionRows;    /* The number of rows in each partition followed by the boundary and residual */
  int          ****rowPartition;        /* The row or edge partition */
  int              rank;                /* The rank of B^p which is the length of range */
  int             *range;               /* The rows of P in the range of B ordered so that D is diagonal */
  VecScatter       rangeScatter;        /* The map from a solution vector to a column vector */
  int              numLocNullCols;      /* The number of columns corresponding to zero singular values */
  int             *nullCols;            /* The columns corresponding to zero singular values */

  /* Factorization variables */
  int              numRows;             /* The number of rows    in B */
  int              numCols;             /* The number of columns in B */
  int              numLevels;           /* The number of levels in the factorization */
  int              maxLevels;           /* The threshold for new allocation */
  int              QRthresh;            /* The node threshold for performing a final QR in this domain */
  PetscReal        DoQRFactor;          /* The multiplier m for row > cols*m which determines when to QR U */
  double           zeroTol;             /* The numeric tolerance used for determining zero singular values */
  PetscReal    ****factors;             /* The factor list */
  Mat             *grads;               /* The boundary gradients for each level */
  Vec             *bdReduceVecs;        /* The work vectors for the boundary+residual rows at each level */
  Vec             *colReduceVecs;       /* The work vectors for the remaining columns at each level */
  Vec             *colReduceVecs2;      /* A duplicate of colReduceVecs[] */
  PetscReal       *interiorWork;        /* A work vector the size of the maximum rows in any partition */
  PetscReal       *interiorWork2;       /* A duplicate of interiorWork[] */
  int              interiorWorkLen;     /* The size of interiorWork[] */

  /* Parallel variables */
  int              globalRank;          /* The rank of B, or total number of range space vectors */
  int              localRank;           /* The rank of B_p */
  int              numNullCols;         /* The total number of columns corresponding to zero singular values */
  int             *firstNullCol;        /* The first column with a zero singular value in each domain */
  int              numLocInterfaceRows; /* The number of rows with nonzeros in another domain */
  int              numInterfaceRows;    /* The global number of interface rows */
  int             *interfaceRows;       /* The row numbers of interface rows in the original matrix */
  int             *firstInterfaceRow;   /* The first interface row in each domain */
  int              numInteriorRows;     /* The number of rows which are not interface rows in this domain */
  int             *interiorRows;        /* The row numbers of interior rows in the original matrix */
#ifdef PETSC_HAVE_PLAPACK
  PLA_Template     PLAlayout;           /* Layout of PLA structures across processors */
  PLA_Obj          interfaceQR;         /* The QR factorization of the interface matrix B_I */
  PLA_Obj          interfaceTAU;        /* The elementary reflection multipliers for interfaceQR */
  PLA_Obj          PLArhsD;             /* The work vector for applying D size matrices */
  PLA_Obj          PLArhsP;             /* The work vector for applying P size matrices */
#else
  PetscReal       *interfaceQR;         /* The QR factorization of the interface matrix B_I */
  PetscReal       *interfaceTAU;        /* The elementary reflection multipliers for interfaceQR */
#ifdef CHOLESKY_SCHEME
  PetscReal       *interfaceGrad;       /* The interface matrix B_I */
#endif
#endif
  GVec            *interfaceV;          /* The column of V^i corresponding to the zero singular value in that domain */
  Mat              interfaceB;          /* The interface rows of B for each domain */
  Vec              interiorRhs;         /* The rhs for solves in the interior */
  VecScatter       interiorScatter;     /* The map from a global solution vector to an interior vector */
#ifndef PETSC_HAVE_PLAPACK
  Vec              locInterfaceRhs;     /* The rhs for reduction of the local interface */
  VecScatter       locInterfaceScatter; /* The map from a global solution vector to a local interface vector */
#endif
  Vec              interfaceRhs;        /* The rhs for solves on the interface: only on proc 1 */
  VecScatter       interfaceScatter;    /* The map from a global solution vector to an interface vector: only on proc 1 */
  Vec              interfaceColRhs;     /* The rhs for solves on the column interface: only on proc 1 */
  VecScatter       interfaceColScatter; /* The map from a global column vector to a column interface vector: only on proc 1 */
  VecScatter       reduceScatter;       /* The map from a parallel column vector to a serial column vector on each proc */
  GVec             colWorkVec;          /* A Work vector in the column space of B */
  Vec              globalColWorkVec;    /* A work vector the size of colWorkVec, but duplicated on each processor */

  /* Work space for LAPACK routines */
  int          workLen;
  PetscScalar *work;
  int          svdWorkLen;
  PetscScalar *svdWork;
} PC_Multilevel;

/* mlSerial prototypes */
extern PetscTruth PCMultiLevelDoQR_Private(PC, int, int);
extern int PCMultiLevelCalcVertexCoords_Private(PC, int, int, int *);
extern int PCMultiLevelFilterVertexCoords_Tutte(PC, int, double *, int **);
extern int PCMultiLevelReduceView(PC, int, int, int *, int, int *);
extern int PCMultiLevelReduce(PC);
extern int PCMultiLevelView_Private(PC, PetscViewer, int, int, int, int **);

/* mlApply prototypes */
extern int PCMultiLevelApplyDInv(PC, GVec, GVec);
extern int PCMultiLevelApplyDInvTrans(PC, GVec, GVec);
extern int PCMultiLevelApplyV(PC, GVec, GVec);
extern int PCMultiLevelApplyVTrans(PC, GVec, GVec);
extern int PCMultiLevelApplyP(PC, GVec, GVec);
extern int PCMultiLevelApplyPTrans(PC, GVec, GVec);
extern int PCMultiLevelApplyP1(PC, GVec, GVec);
extern int PCMultiLevelApplyP1Trans(PC, GVec, GVec);
extern int PCMultiLevelApplyP2(PC, GVec, GVec);
extern int PCMultiLevelApplyP2Trans(PC, GVec, GVec);

/* mlCheck prototypes */
extern int PCValidQ_Multilevel(PC);
extern int PCDebug_Multilevel(PC);
extern int PCDebugPrintMat_Multilevel(PC, const char []);

#endif
