#define _USE_MATH_DEFINES
static char help[] = "Solves a linear system in parallel with KSP.\n\n";
/*
   Description: Solves a real linear system in parallel with KSP.

   The model problem:
      Solve Helmholtz equation on the unit square: (0,1) x (0,1)
          -delta phi - k2*phi = f,
           where delta = Laplace operator
      Dirichlet b.c.'s on all sides
      Use the 2-D, five-point finite difference stencil.
*/
#include <iostream>
#include <cstdio>
#include <ctime>


#include <petscksp.h>
// #include <iostream>
#include <math.h>

int main(int argc,char **args)
{
  // Timer definition
  std::clock_t start;
  double duration;
  start = std::clock();


  Vec            phi_num,b,phi, phi_c, phi_numc;    /* approx solution, RHS, exact solution */
  Mat            A;                /* finite difference matrix */
  KSP            ksp;              /* linear solver */
  PetscInt       dim,i,j,Ii,J,Istart,Iend,n = 200,its,size;
  PetscErrorCode ierr;
  PetscScalar    v,*xa,*array,norm,dx2,k2 = 1.0;
  PetscBool      flg = PETSC_FALSE;

  ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
  ierr = PetscOptionsGetReal(NULL,NULL,"-k2",&k2,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
  dim  = n*n;

  ierr = PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  /*
     Create parallel matrix, specifying only its global dimensions.
     When using MatCreate(), the matrix format can be specified at
     runtime. Also, the parallel partitioning of the matrix is
     determined by PETSc at runtime.
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);

  /*
     Currently, all PETSc parallel matrix formats are partitioned by
     contiguous chunks of rows across the processors.  Determine which
     rows of the matrix are locally owned.
  */
  ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);

  /*
     Set matrix elements in parallel.
      - Each processor needs to insert only elements that it owns
        locally (but any non-local elements will be sent to the
        appropriate processor during matrix assembly).
      - Always specify global rows and columns of matrix entries.
  */

  dx2 = 1.0/((n-1)*(n-1));
  for (Ii=Istart; Ii<Iend; Ii++) {
    v = -1.0; i = Ii/n; j = Ii - i*n;
    if (i>0) {
      J = Ii-n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
    }
    if (i<n-1) {
      J = Ii+n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
    }
    if (j>0) {
      J = Ii-1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
    }
    if (j<n-1) {
      J = Ii+1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
    }
    v    = 4.0 - k2*dx2;
    ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
  }

  /*
     Assemble matrix, using the 2-step process:
       MatAssemblyBegin(), MatAssemblyEnd()
     Computations can be done while messages are in transition
     by placing code between these two statements.
  */
  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /*
     Create parallel vectors.
      - When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
      we specify only the vector's global
        dimension; the parallel partitioning is determined at runtime.
      - Note: We form 1 vector from scratch and then duplicate as needed.
  */

  ierr = VecCreate(PETSC_COMM_WORLD,&phi);CHKERRQ(ierr);
  ierr = VecSetSizes(phi,PETSC_DECIDE,dim);CHKERRQ(ierr);
  ierr = VecSetFromOptions(phi);CHKERRQ(ierr);
  ierr = VecDuplicate(phi,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&phi_num);CHKERRQ(ierr);

  /*
     Set exact solution; then compute right-hand-side vector.
  */

  for (int col = 0; col < n; ++col) {
    for (int row = 0; row < n; ++row) {
      PetscInt indx = row*n+col;
      PetscScalar value = sin(M_PI/2.0*col/(n-1.0))*cos(M_PI/2.0*row/(n-1.0))+5.0;
      // PetscScalar value = col;
      ierr = VecSetValues(phi,1,&indx,&value,INSERT_VALUES);CHKERRQ(ierr);

      // set phi_num bcs
      if ((row == 0)||(row == n-1)||(col == 0)||(col == n-1)) {
        ierr = VecSetValues(phi_num,1,&indx,&value,INSERT_VALUES);CHKERRQ(ierr);
      }

    }
  }

  ierr = VecGetSize(phi,&size);CHKERRQ(ierr);
  ierr = VecGetArray(phi,&array);CHKERRQ(ierr);

  ierr = MatMult(A,phi,b);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                Create the linear solver and set various options
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  ierr = KSPSetOptionsPrefix(ksp,"ksp_atol 1e-12");
  ierr = KSPSetOptionsPrefix(ksp,"ksp_rtol 1e-9");

  /*
    Set runtime options, e.g.,
        -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
  */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = KSPSolve(ksp,b,phi_num);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Check solution and clean up
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  // std::cout << "phi mat";
  // VecView(phi_num,PETSC_VIEWER_STDOUT_WORLD);

  /*
      Print the first 3 entries of phi_num; this demonstrates extraction of the
      real and imaginary components of the complex vector, phi_num.
  */
  flg  = PETSC_TRUE;
  ierr = PetscOptionsGetBool(NULL,NULL,"-print_x3",&flg,NULL);CHKERRQ(ierr);
  if (flg) {
    ierr = VecGetArray(phi_num,&xa);CHKERRQ(ierr);
    //output ierr = PetscPrintf(PETSC_COMM_WORLD,"The first three entries of phi_num are:\n");CHKERRQ(ierr);
    for (i=0; i<10; i++) {
      //output ierr = PetscPrintf(PETSC_COMM_WORLD,"phi_num[%D] = %g \n",i,(double)PetscRealPart(xa[i]));CHKERRQ(ierr);
    }
    ierr = VecRestoreArray(phi_num,&xa);CHKERRQ(ierr);
  }

  /*
     Check the error
  */


  ierr = VecDuplicate(phi,&phi_c);CHKERRQ(ierr);
  ierr = VecDuplicate(phi_num,&phi_numc);CHKERRQ(ierr);
  ierr = VecCopy(phi,phi_c);CHKERRQ(ierr);
  ierr = VecCopy(phi_num,phi_numc);CHKERRQ(ierr);

  // // Prev
  // ierr = VecAXPY(phi_c,none,phi_numc);CHKERRQ(ierr); // phi_c = phi_c-phi_num
  // ierr = VecAbs(phi_c);CHKERRQ(ierr);               // phi_c = |phi_c|
  // ierr = VecPointwiseDivide(phi_c,phi_c,phi);CHKERRQ(ierr); //phi_c = |phi_c-phi_num|./phi
  // ierr = VecSum(phi_c,&norm);

  // New
  ierr = VecAXPY(phi_numc,-1.0,phi_c);CHKERRQ(ierr);
  ierr = VecNorm(phi_numc,NORM_2,&norm);CHKERRQ(ierr);

  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
  //output ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm/(n*n),its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"%g\n",(double)norm/(n*n),its);CHKERRQ(ierr);

  /*
     Free work space
  */
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = VecDestroy(&phi);CHKERRQ(ierr); ierr = VecDestroy(&phi_num);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);   ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = PetscFinalize();


// Print time usage
  duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
  //output std::cout<<"Elapsed: "<< duration <<"s\n";
  return ierr;
}