# include <cstdlib>
# include <iostream>

using namespace std;
//
//  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 );

//****************************************************************************80

int main ( int argc, char **argv )

//****************************************************************************80
//
//  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 )
  {
    cout << "\n";
    cout << "Divergence because of indefinite preconditioner;\n";
    cout << "Run the executable again but with -pc_ilu_shift option.\n";
  }
  else if ( reason < 0 )
  {
    cout << "\n";
    cout << "Other kind of divergence: this should not happen.\n";
  }
  else
  {
    KSPGetIterationNumber ( solver, &its );
    cout << "\n";
    cout << "Convergence in " << ( int ) its << " iterations.\n";
  }
  cout << "\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 );
}