#ifdef PETSC_RCS_HEADER
static char vcid[] = "$Id: mlInterface.c,v 1.1 1999/06/09 20:24:12 knepley Exp $";
#endif
/*
   Defines the interface for the  multilevel preconditioner due to
   Vivek Sarin for AIJ matrices.
*/
#include "src/sles/pc/pcimpl.h"    /*I "pc.h" I*/
#include "petscsnes.h"
#include "ml.h"

#undef  __FUNCT__
#define __FUNCT__ "PCMultiLevelSetFields"
/*@C PCMultiLevelSetFields
	This function sets the shape and test function fields for the
  gradient operator. It must be called prior to PCSetUp().

  Collective on PC

  Input Parameters:
+ pc     - The preconditioner context
. sField - The shape function field
- tField - The test function field

  Level: intermediate

.keywords multilevel, fields
.seealso PCMultiLevelSetGradientOperator, PCMultiLevelSetNonlinearIterate
@*/
int PCMultiLevelSetFields(PC pc, int sField, int tField)
{
  PC_Multilevel *ml;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  ml = (PC_Multilevel *) pc->data;
  ml->sField = sField;
  ml->tField = tField;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCMultiLevelSetNonlinearIterate"
/*@C PCMultiLevelSetNonlinearIterate
	This function sets the solution vector used in a Newton iteration. If
  the inner solve is preconditioned with ML, the current iterate is
  necessary to perform the reduction of the system.

  Collective on PC

  Input Parameters:
+ pc - The preconditioner context
- u  - The current iterate in a Newton solve

  Level: intermediate

.keywords multilevel, nonlinear
.seealso PCMultiLevelSetFields, PCMultiLevelSetGradientOperator
@*/
int PCMultiLevelSetNonlinearIterate(PC pc, GVec u)
{
  PC_Multilevel *ml;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  PetscValidHeaderSpecific(u,  VEC_COOKIE);
  ml = (PC_Multilevel *) pc->data;
  ml->nonlinearIterate = u;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCMultiLevelSetGradientOperator"
/*@C PCMultiLevelSetGradientOperator
	This function sets the operator to use as the gradient and its transpose
  in the multilevel scheme. It must be called prior to PCSetUp().

  Collective on PC

  Input Parameters:
+ pc    - The preconditioner context
. grad  - The gradient operator
. div   - The divergence operator
- alpha - The scalar multiplier for the gradient operator

  Level: intermediate

.keywords multilevel, operator
.seealso PCMultiLevelSetFields, PCMultiLevelSetNonlinearIterate
@*/
int PCMultiLevelSetGradientOperator(PC pc, int Grad, int Div, PetscScalar alpha)
{
  PC_Multilevel *ml;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  ml = (PC_Multilevel *) pc->data;
  ml->gradOp    = Grad;
  ml->divOp     = Div;
  ml->gradAlpha = alpha;
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCMultiLevelBuildSolution"
/*@C PCMultiLevelBuildSolution
	This function recovers the solution of the global problem
  from the solution of the reduced system.

  Collective on GVec

  Input Parameters:
+ pc - The preconditioner context
- u  - The solution of the reduced system

  Output Parameter:
. u  - The global solution

  Level: beginner

.keywords multilevel, solution
.seealso PCMultiLevelApplyGradientGetMultiplier
@*/
int PCMultiLevelBuildSolution(PC pc, GVec u)
{
  PC_Multilevel *ml;
  Grid           grid;
  PetscScalar    alpha;
  int            ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  PetscValidHeaderSpecific(u,  VEC_COOKIE);
  if (pc->setupcalled < 2) {
    ierr = PCSetUp(pc);                                                                                  CHKERRQ(ierr);
  }

  ml = (PC_Multilevel *) pc->data;
  if (ml->reduceSol != PETSC_NULL) {
    ierr  = GVecGetGrid(u, &grid);                                                                       CHKERRQ(ierr);
    ierr  = GridGetBCMultiplier(grid, &alpha);                                                           CHKERRQ(ierr);
    alpha = -alpha;
    ierr  = VecAXPY(&alpha, ml->reduceSol, u);                                                           CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

#undef  __FUNCT__
#define __FUNCT__ "PCMultiLevelGetMultiplier"
/*@C PCMultiLevelGetMultiplier
	This function recovers the multiplier of the global problem
  from the solution of the reduced system.

  Collective on GVec

  Input Parameters:
+ pc - The preconditioner context
- u  - The solution of the reduced system

  Output Parameter:
. p  - The global multiplier

  Level: beginner

.keywords multilevel, solution, multiplier
.seealso PCMultiLevelBuildSolution
@*/
int PCMultiLevelGetMultiplier(PC pc, GVec u, GVec p)
{
  GVec        tempVec;
  PetscScalar minusOne = -1.0;
  PetscScalar zero     =  0.0;
  int         ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_COOKIE);
  PetscValidHeaderSpecific(u,  VEC_COOKIE);
  PetscValidHeaderSpecific(p,  VEC_COOKIE);
  if ((pc->vec == PETSC_NULL) || (pc->mat == PETSC_NULL)) {
    ierr = VecSet(&zero, p);                                                                             CHKERRQ(ierr);
    PetscFunctionReturn(0);
  }
  if (pc->setupcalled < 2) {
    ierr = PCSetUp(pc);                                                                                  CHKERRQ(ierr);
  }

  /* Multiply P_2 (P^{-1} u)_2 by A */
  ierr = GVecDuplicate(u, &tempVec);                                                                     CHKERRQ(ierr);
  ierr = MatMult(pc->mat, u, tempVec);                                                                   CHKERRQ(ierr);

  /* Let A P_2 (P^{-1} u)_2 = f - A P_1 D^{-1} Z^T g - A P_2 (P^{-1} u)_2 --
       This is really unstaisfactory, but how do I get the Rhs after solution of
       the system. Perhaps I should have a flag for storing the multiplier.
  */
  ierr = VecAYPX(&minusOne, pc->vec, tempVec);                                                           CHKERRQ(ierr);

  /* Apply P^T_1 */
  ierr = PCMultiLevelApplyP1Trans(pc, tempVec, p);                                                       CHKERRQ(ierr);

  /* Apply Z D^{-1} */
  ierr = PCMultiLevelApplyDInv(pc, p, p);                                                                CHKERRQ(ierr);
  ierr = PCMultiLevelApplyV(pc, p, p);                                                                   CHKERRQ(ierr);

  /* Cleanup */
  ierr = VecDestroy(tempVec);                                                                            CHKERRQ(ierr);
  PetscFunctionReturn(0);
}