/* 
 * This example code is modified from ./didasko/examples/aztecoo/ex2.cpp in 
 * Trilinos-8.0.8 for testing purpose only.
 *
 * The original header is attached below, in case of any copyright issue.
 */

// @HEADER
// // ***********************************************************************
// // 
// //                      Didasko Tutorial Package
// //                 Copyright (2005) Sandia Corporation
// // 
// // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// // license for use of this work by or on behalf of the U.S. Government.
// // 
// // This library is free software; you can redistribute it and/or modify
// // it under the terms of the GNU Lesser General Public License as
// // published by the Free Software Foundation; either version 2.1 of the
// // License, or (at your option) any later version.
// //  
// // This library is distributed in the hope that it will be useful, but
// // WITHOUT ANY WARRANTY; without even the implied warranty of
// // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// // Lesser General Public License for more details.
// //  
// // You should have received a copy of the GNU Lesser General Public
// // License along with this library; if not, write to the Free Software
// // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// // USA
// //
// // Questions about Didasko? Contact Marzio Sala (marzio.sala _AT_ gmail.com)
// // 
// // ***********************************************************************
// // @HEADER
//
// //  Solve a linear system with Aztec, using Aztec as a preconditioner
// // (recursive way)
// //
// // NOTE: This example first builds the matrix, then solves it with AztecOO
// //
// // NOTE2: this example implemenets minor modifications to one of the
// // examples included in the AztecOO package. Please give a look  to file
// // ${TRILINOS_HOME}/packages/aztecoo/examples/AztecOO_RecursiveCall/cxx_main.cpp
// // for more details.
//
//

#define HAVE_MPI

#include "Epetra_MpiComm.h"
#include "mpi.h"

#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "AztecOO_Operator.h"

/* this example creates a tridiagonal matrix of type
 */

int main(int argc, char *argv[])
{

  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);

  int NumGlobalElements = 50;
  Epetra_Map Map(NumGlobalElements,0,Comm);
  int NumMyElements = Map.NumMyElements();
  int * MyGlobalElements = new int [NumMyElements];
  Map.MyGlobalElements( MyGlobalElements );

  Epetra_CrsMatrix A(Copy,Map,3);

  double *Values = new double[2];
  Values[0] = -1.0; Values[1] = -1.0;
  int *Indices = new int[2];
  double two = 2.0;
  int NumEntries;

  for( int i=0 ; i<NumMyElements; ++i ) {
    if (MyGlobalElements[i]==0) {                  
	Indices[0] = 1;
	NumEntries = 1;
    } else if (MyGlobalElements[i] == NumGlobalElements-1) { 
      Indices[0] = NumGlobalElements-2;
      NumEntries = 1;
    } else {                                               
      Indices[0] = MyGlobalElements[i]-1;
      Indices[1] = MyGlobalElements[i]+1;
      NumEntries = 2;
    }
    A.InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
    A.InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i);
  }
  
  // Finish up
  A.FillComplete();

  Epetra_Vector x(Map);
  Epetra_Vector b(Map);

  x.Random();
  Epetra_LinearProblem A_Problem(&A, &x, &b);
  AztecOO A_Solver(A_Problem);
  Epetra_CrsMatrix P(A);

  Epetra_LinearProblem P_Problem;
  P_Problem.SetOperator(&P);
  AztecOO P_Solver(P_Problem);
  P_Solver.SetAztecOption(AZ_precond, AZ_Jacobi);
  P_Solver.SetAztecOption(AZ_output, AZ_none);
  P_Solver.SetAztecOption(AZ_solver, AZ_cg);
  AztecOO_Operator P_Operator(&P_Solver, 10);
  A_Solver.SetPrecOperator(&P_Operator);
  
  // solve
  
  int Niters=1000;
  A_Solver.SetAztecOption(AZ_kspace, Niters);
  A_Solver.SetAztecOption(AZ_solver, AZ_GMRESR);

  A_Solver.Iterate(Niters, 1.0E-12);

  MPI_Finalize() ;

  return 0;
}