# include <stdlib.h>
/*
  When using PETSC, you usually only have to point to the highest level include
  file, and all the lower level ones will be included automatically.
*/
# include "petscksp.h"
/*
  I guess the following pair of lines are intended to tell PETSc 
  that the main program is called "main".
*/
# undef __FUNCT__
# define __FUNCT__ "main"

int main ( int argc, char **argv );

/*********************************************************************/

int main ( int argc, char **argv )

/*********************************************************************/
/*
  Purpose:

    MAIN is a demonstration of the use of PETSC on a linear system.

  Discussion:

    Test file for the PCILUSetShift routine or -pc_shift option.
    The test matrix is the example from Kershaw's paper 
    of a positive definite matrix for which ILU(0) will give a negative pivot.
    This means that the CG method will break down; the Manteuffel shift
    repairs this.

    Run the executable twice:

    1, with the call to PCFactorSetShiftPd turned off,
    the iterative method diverges because of an indefinite preconditioner;

    2, with the call to PCFactorSetShiftPd turned on,
    the method will now converge.

  Modified:

    01 November 2005

  Author:

    Victor Eijkhout

  Reference:

    David Kershaw,
    The Incomplete Cholesky-Conjugate Gradient Method for the Iterative
    Solution of Systems of Linear Equations,
    Journal of Computational Physics,
    Volume 26, pages 43-65, 1978.

    Tom Manteuffel,
    An Incomplete Factorization Technique for Positive Definite Linear Systems,
    Mathematics of Computation,
    Volume 34, pages 473-497, 1980.
 */
{
  Mat A;
  Vec B;
  MPI_Comm comm;
  Vec D;
  PetscInt i;
  PetscErrorCode ierr;
  PetscInt its;
  PetscInt j;
  Mat M;
  PC prec; 
  KSPConvergedReason reason;
  KSP solver;
  PetscScalar v; 
  Vec X;

  PetscFunctionBegin;
/*
  Initialize PETSc.  
  IF MPI has not been initialized already, this call will do that as well.
*/
  PetscInitialize ( &argc, &argv, 0, 0 );

  PetscOptionsSetValue ( "-options_left", PETSC_NULL );

  comm = MPI_COMM_SELF;
/*
  Construct the Kershaw matrix A
  and a suitable right hand side B and initial solution guess X.
*/
  MatCreateSeqAIJ ( comm, 4, 4, 4, 0, &A );

  VecCreateSeq ( comm, 4, &B );

  VecDuplicate ( B, &X );

  for ( i = 0; i < 4; i++ ) 
  {
    v = 3;
    MatSetValues ( A, 1, &i, 1, &i, &v, INSERT_VALUES );
    v = 1;
    VecSetValues ( B, 1, &i, &v, INSERT_VALUES );
    VecSetValues ( X, 1, &i, &v, INSERT_VALUES );
  }

  i = 0;
  v = 0;

  VecSetValues ( X, 1, &i, &v, INSERT_VALUES );

  for ( i = 0; i < 3; i++ ) 
  {
    v = -2;
    j = i + 1;

    MatSetValues ( A, 1, &i, 1, &j, &v, INSERT_VALUES );

    MatSetValues ( A, 1, &j, 1, &i, &v, INSERT_VALUES );
  }

  i = 0;
  j = 3;
  v = 2;

  MatSetValues ( A, 1, &i, 1, &j, &v, INSERT_VALUES );

  MatSetValues ( A, 1, &j, 1, &i, &v, INSERT_VALUES );

  MatAssemblyBegin ( A, MAT_FINAL_ASSEMBLY );
 
  MatAssemblyEnd ( A, MAT_FINAL_ASSEMBLY );

  VecAssemblyBegin ( B );

  VecAssemblyEnd ( B );

  printf ( "\n"); 
  printf ( "The Kershaw matrix A:\n"); 

  MatView ( A, 0 );
/*
  Set up a Conjugate Gradient method with ILU(0) preconditioning.
*/
  KSPCreate ( comm, &solver );

  KSPSetOperators ( solver, A, A, SAME_NONZERO_PATTERN );

  KSPSetType ( solver, KSPCG );

  KSPSetInitialGuessNonzero ( solver, PETSC_TRUE );
/*
  The ILU preconditioner will break down unless you invoke the PCFactorSetShiftPd
  routine or use the -pc_ilu_shift option.
*/
  KSPGetPC ( solver, &prec );

  PCSetType ( prec, PCILU );
/*
  Turning on this call will avoid the algorithmic breakdown
  caused by the negative pivot.
*/
  if ( 1 )
  {
    PCFactorSetShiftPd ( prec, PETSC_TRUE );
  }

  KSPSetFromOptions ( solver );

  KSPSetUp ( solver );
/*
  Now that the factorization is done, show the pivots;
  note that the last one is negative. This in itself is not an error,
  but it will make the iterative method diverge.
*/
  PCGetFactoredMatrix ( prec, &M );

  VecDuplicate ( B, &D );

  MatGetDiagonal ( M, D );

  printf ( "\n" );
  printf ( "Pivots:\n"); 

  VecView ( D, 0 );
/*
  Solve the system; without the shift this will diverge with
  an indefinite preconditioner
*/
  KSPSolve ( solver, B, X );
/*
  Retrieve the reason that the solver halted, and print it out.
*/
  KSPGetConvergedReason ( solver, &reason );

  if ( reason == KSP_DIVERGED_INDEFINITE_PC )
  {
    printf ( "\n" );
    printf ( "Divergence because of indefinite preconditioner;\n" );
    printf ( "Run the executable again but with -pc_ilu_shift option.\n" );
  }
  else if ( reason < 0 )
  {
    printf ( "\n" );
    printf ( "Other kind of divergence: this should not happen.\n" );
  }
  else
  {
    KSPGetIterationNumber ( solver, &its );
    printf ( "\n" );
    printf ( "Convergence in %d iterations.\n", ( int ) its );
  }
  printf ( "\n" );
/*
  Free up the memory.
*/
  VecDestroy ( X );
  VecDestroy ( B );
  VecDestroy ( D );
  MatDestroy ( A );
  KSPDestroy ( solver );
/*
  Inform PETSC that we are done with it.
  If MPI is not already shut down, this command will take care of that as well.
*/
  PetscFinalize ( );

  PetscFunctionReturn ( 0 );
}