/*
static char help[] = "Solves a linear system in parallel with KSP.\n\
Input parameters include:\n\
  -random_exact_sol : use a random exact solution vector\n\
  -view_exact_sol   : write exact solution vector to stdout\n\
  -m <mesh_x>       : number of mesh points in x-direction\n\
  -n <mesh_n>       : number of mesh points in y-direction\n\n";

*/
/*
Input parameters include:\n\
  -random_exact_sol : use a random exact solution vector\n\
  -view_exact_sol   : write exact solution vector to stdout\n\
  -m <mesh_x>       : number of mesh points in x-direction\n\
  -n <mesh_n>       : number of mesh points in y-direction\n\n";
*/

/*T
   Concepts: KSP^basic parallel example;
   Concepts: KSP^Laplacian, 2d
   Concepts: Laplacian, 2d
   Processors: n
T*/

/*
  Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  automatically includes:
     petscsys.h       - base PETSc routines   petscvec.h - vectors
     petscmat.h - matrices
     petscis.h     - index sets            petscksp.h - Krylov subspace methods
     petscviewer.h - viewers               petscpc.h  - preconditioners
*/
#include <petscksp.h>
#include <iostream>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
  Vec            x,b;  /* approx solution, RHS, exact solution */
  Mat            A;        /* linear system matrix */
  KSP            ksp;     /* linear solver context */
//  PetscRandom    rctx;     /* random number generator context */
//  PetscReal      norm;     /* norm of solution error */
  PetscInt       i,j,m = 3,n = 3,its;
  PetscErrorCode ierr;
//  PetscBool      flg = PETSC_FALSE;
  PetscScalar    v;
//#if defined(PETSC_USE_LOG)
//  PetscLogStage stage;
//#endif


  PetscInitialize(&argc,&args,(char*)0,NULL);

  // print value of an option 
  const char* pre = NULL;
  const char name[] = "-ksp_pc_side";
  const size_t len = 256;
  char str[len];
  PetscBool isset;

  ierr = PetscOptionsGetString(pre, name, str, len, &isset);
  CHKERRQ(ierr);

  if (isset)
  {
    std::cout << "option " << name << " was found" << std::endl;
    std::cout << "returned string is: " << str << std::endl;
  } else
  {
    std::cout << "option " << name << " not found" << std::endl;
  }

    

  
   
 
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr);

//  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);

//  ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr);
//  ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
  for (i=0; i<m; i++) {
    for (j=0; j<n; j++) {
      v = i*n + j + 1;
      std::cout << "inserting value " << v << " into position " << i << ", " << j << std::endl;
      MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES); CHKERRQ(ierr);
    }
  }

  // change one entry so it is non singular
  i = 2;
  j = 0;
  v = 8;

  std::cout << "inserting value " << v << " into position " << i << ", " << j << std::endl;
  MatSetValues(A, 1, &i, 1, &j, &v, INSERT_VALUES); CHKERRQ(ierr);

  // assemble matrix
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
//  ierr = PetscLogStagePop();CHKERRQ(ierr);


  MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

  /*
     Create parallel vectors.
      - We form 1 vector from scratch and then duplicate as needed.
      - When using VecCreate(), VecSetSizes and VecSetFromOptions()
        in this example, we specify only the
        vector's global dimension; the parallel partitioning is determined
        at runtime.
      - When solving a linear system, the vectors and matrices MUST
        be partitioned accordingly.  PETSc automatically generates
        appropriately partitioned matrices and vectors when MatCreate()
        and VecCreate() are used with the same communicator.
      - The user can alternatively specify the local vector and matrix
        dimensions when more sophisticated partitioning is needed
        (replacing the PETSC_DECIDE argument in the VecSetSizes() statement
        below).
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecSetSizes(b,PETSC_DECIDE,m);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
//  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

  // set rhs
  for (i=0; i < m; ++i)
  {
    v = i;  //convert to scalar
    ierr = VecSetValues(b, 1, &i, &v, INSERT_VALUES);CHKERRQ(ierr);
  }

  ierr = VecAssemblyBegin(b); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b); CHKERRQ(ierr);

  ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  ierr = KSPSetTolerances(ksp,1.e-12,1.e-50,PETSC_DEFAULT,
                          PETSC_DEFAULT);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  // solve
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);


  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); 
  ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %d\n", its);CHKERRQ(ierr);
  

  Mat B, C;
  ierr = MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &B); CHKERRQ(ierr);

  std::cout << "pointer C = " << C << std::endl;
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Mat B = \n"); CHKERRQ(ierr);
  ierr = MatView(B, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);


  // multiply matrices
  std::cout << "sizeof(MAT_INITIAL_MATRIX) = " << sizeof(MAT_INITIAL_MATRIX) << std::endl;

  ierr = MatMatMult(A, B, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C); CHKERRQ(ierr);

  // view result
  ierr = PetscPrintf(PETSC_COMM_WORLD, "Mat C = \n"); CHKERRQ(ierr);
  ierr = MatView(C, PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

  // cleanup
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = PetscFinalize();

  return 0;
}