#ifndef HPCG_NO_MPI
#include <mpi.h>
#endif
#include <fstream>
#include <iostream>
#include <cfloat>
using std::endl;
#include <vector>
#include <cmath>

#include "hpcg.hpp"

#include "ComputeSPMV.hpp"
#include "ComputeMG.hpp"
#include "ComputeDotProduct.hpp"
#include "ComputeResidual.hpp"
#include "Geometry.hpp"
#include "SparseMatrix.hpp"
#include "TestSymmetry.hpp"

#ifdef HPCG_WITH_MORPHEUS
#include "MorpheusUtils.hpp"
#endif  

int TestSymmetry(SparseMatrix& A, Vector& b, Vector& xexact,
                 TestSymmetryData& testsymmetry_data) {
  local_int_t nrow = A.localNumberOfRows;
  local_int_t ncol = A.localNumberOfColumns;

  Vector x_ncol, y_ncol, z_ncol;
  InitializeVector(x_ncol, ncol);
  InitializeVector(y_ncol, ncol);
  InitializeVector(z_ncol, ncol);

  double t4 = 0.0;  
  testsymmetry_data.count_fail = 0;

  FillRandomVector(x_ncol);
  FillRandomVector(y_ncol);

  double xNorm2, yNorm2;
  double ANorm = 2 * 26.0;

#ifdef HPCG_WITH_MORPHEUS
  MorpheusInitializeVector(x_ncol);
  MorpheusInitializeVector(y_ncol);
  MorpheusInitializeVector(z_ncol);

  MorpheusOptimizeVector(x_ncol);
  MorpheusOptimizeVector(y_ncol);
  MorpheusOptimizeVector(z_ncol);

  HPCG_Morpheus_Vec* xncolopt  = (HPCG_Morpheus_Vec*)x_ncol.optimizationData;
  HPCG_Morpheus_Vec* yncolopt  = (HPCG_Morpheus_Vec*)y_ncol.optimizationData;
  HPCG_Morpheus_Vec* zncolopt  = (HPCG_Morpheus_Vec*)z_ncol.optimizationData;
  HPCG_Morpheus_Vec* xexactopt = (HPCG_Morpheus_Vec*)xexact.optimizationData;
#endif

  ComputeDotProduct(nrow, y_ncol, y_ncol, yNorm2, t4, A.isDotProductOptimized);
  int ierr = ComputeSPMV(A, y_ncol, z_ncol);  
  if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
  double xtAy = 0.0;
  ierr        = ComputeDotProduct(nrow, x_ncol, z_ncol, xtAy, t4,
                                  A.isDotProductOptimized);  
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  ComputeDotProduct(nrow, x_ncol, x_ncol, xNorm2, t4, A.isDotProductOptimized);
  ierr = ComputeSPMV(A, x_ncol, z_ncol);  
  if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
  double ytAx = 0.0;
  ierr        = ComputeDotProduct(nrow, y_ncol, z_ncol, ytAx, t4,
                                  A.isDotProductOptimized);  
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  testsymmetry_data.depsym_spmv =
      std::fabs((long double)(xtAy - ytAx)) /
      ((xNorm2 * ANorm * yNorm2 + yNorm2 * ANorm * xNorm2) * (DBL_EPSILON));
  if (testsymmetry_data.depsym_spmv > 1.0)
    ++testsymmetry_data.count_fail;  
  if (A.geom->rank == 0)
    HPCG_fout
        << "Departure from symmetry (scaled) for SpMV abs(x'*A*y - y'*A*x) = "
        << testsymmetry_data.depsym_spmv << endl;

#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(yncolopt->dev, yncolopt->host);
#endif

  ierr = ComputeMG(A, y_ncol, z_ncol);  
  if (ierr) HPCG_fout << "Error in call to MG: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(zncolopt->host, zncolopt->dev);
#endif
  double xtMinvy = 0.0;
  ierr           = ComputeDotProduct(nrow, x_ncol, z_ncol, xtMinvy, t4,
                                     A.isDotProductOptimized);  
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(xncolopt->dev, xncolopt->host);
#endif

  ierr = ComputeMG(A, x_ncol, z_ncol);  
  if (ierr) HPCG_fout << "Error in call to MG: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(zncolopt->host, zncolopt->dev);
#endif
  double ytMinvx = 0.0;
  ierr           = ComputeDotProduct(nrow, y_ncol, z_ncol, ytMinvx, t4,
                                     A.isDotProductOptimized);  
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  testsymmetry_data.depsym_mg =
      std::fabs((long double)(xtMinvy - ytMinvx)) /
      ((xNorm2 * ANorm * yNorm2 + yNorm2 * ANorm * xNorm2) * (DBL_EPSILON));
  if (testsymmetry_data.depsym_mg > 1.0)
    ++testsymmetry_data.count_fail;  
  if (A.geom->rank == 0)
    HPCG_fout << "Departure from symmetry (scaled) for MG abs(x'*Minv*y - "
                 "y'*Minv*x) = "
              << testsymmetry_data.depsym_mg << endl;

#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(xexactopt->host, xncolopt->host, 0, xexactopt->host.size());
#else
  CopyVector(xexact, x_ncol);  
#endif

  int numberOfCalls = 2;
  double residual   = 0.0;
  for (int i = 0; i < numberOfCalls; ++i) {
    ierr = ComputeSPMV(A, x_ncol, z_ncol);  
    if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS

    Morpheus::copy(zncolopt->dev, zncolopt->host);
#endif
    if ((ierr = ComputeResidual(A.localNumberOfRows, b, z_ncol, residual)))
      HPCG_fout << "Error in call to compute_residual: " << ierr << ".\n"
                << endl;
    if (A.geom->rank == 0)
      HPCG_fout << "SpMV call [" << i << "] Residual [" << residual << "]"
                << endl;
  }
  DeleteVector(x_ncol);
  DeleteVector(y_ncol);
  DeleteVector(z_ncol);

  return 0;
}
//@HEADER
// ***************************************************
//
// HPCG: High Performance Conjugate Gradient Benchmark
//
// Contact:
// Michael A. Heroux ( maherou@sandia.gov)
// Jack Dongarra     (dongarra@eecs.utk.edu)
// Piotr Luszczek    (luszczek@eecs.utk.edu)
//
// ***************************************************
//@HEADER

/*!
 @file TestSymmetry.cpp

 HPCG routine
 */

// The MPI include must be first for Windows platforms
#ifndef HPCG_NO_MPI
#include <mpi.h>
#endif
#include <fstream>
#include <iostream>
#include <cfloat>
using std::endl;
#include <vector>
#include <cmath>

#include "hpcg.hpp"

#include "ComputeSPMV.hpp"
#include "ComputeMG.hpp"
#include "ComputeDotProduct.hpp"
#include "ComputeResidual.hpp"
#include "Geometry.hpp"
#include "SparseMatrix.hpp"
#include "TestSymmetry.hpp"

#ifdef HPCG_WITH_MORPHEUS
#include "MorpheusUtils.hpp"
#endif  // HPCG_WITH_MORPHEUS

/*!
  Tests symmetry-preserving properties of the sparse matrix vector multiply and
  multi-grid routines.

  @param[in]    geom   The description of the problem's geometry.
  @param[in]    A      The known system matrix
  @param[in]    b      The known right hand side vector
  @param[in]    xexact The exact solution vector
  @param[inout] testsymmetry_data The data structure with the results of the CG
  symmetry test including pass/fail information

  @return returns 0 upon success and non-zero otherwise

  @see ComputeDotProduct
  @see ComputeDotProduct_ref
  @see ComputeSPMV
  @see ComputeSPMV_ref
  @see ComputeMG
  @see ComputeMG_ref
*/
int TestSymmetry(SparseMatrix& A, Vector& b, Vector& xexact,
                 TestSymmetryData& testsymmetry_data) {
  local_int_t nrow = A.localNumberOfRows;
  local_int_t ncol = A.localNumberOfColumns;

  Vector x_ncol, y_ncol, z_ncol;
  InitializeVector(x_ncol, ncol);
  InitializeVector(y_ncol, ncol);
  InitializeVector(z_ncol, ncol);

  double t4 = 0.0;  // Needed for dot-product call, otherwise unused
  testsymmetry_data.count_fail = 0;

  // Test symmetry of matrix

  // First load vectors with random values
  FillRandomVector(x_ncol);
  FillRandomVector(y_ncol);

  double xNorm2, yNorm2;
  double ANorm = 2 * 26.0;

#ifdef HPCG_WITH_MORPHEUS
  MorpheusInitializeVector(x_ncol);
  MorpheusInitializeVector(y_ncol);
  MorpheusInitializeVector(z_ncol);

  MorpheusOptimizeVector(x_ncol);
  MorpheusOptimizeVector(y_ncol);
  MorpheusOptimizeVector(z_ncol);

  HPCG_Morpheus_Vec* xncolopt  = (HPCG_Morpheus_Vec*)x_ncol.optimizationData;
  HPCG_Morpheus_Vec* yncolopt  = (HPCG_Morpheus_Vec*)y_ncol.optimizationData;
  HPCG_Morpheus_Vec* zncolopt  = (HPCG_Morpheus_Vec*)z_ncol.optimizationData;
  HPCG_Morpheus_Vec* xexactopt = (HPCG_Morpheus_Vec*)xexact.optimizationData;
#endif

  // Next, compute x'*A*y
  ComputeDotProduct(nrow, y_ncol, y_ncol, yNorm2, t4, A.isDotProductOptimized);
  int ierr = ComputeSPMV(A, y_ncol, z_ncol);  // z_nrow = A*y_overlap
  if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
  double xtAy = 0.0;
  ierr        = ComputeDotProduct(nrow, x_ncol, z_ncol, xtAy, t4,
                                  A.isDotProductOptimized);  // x'*A*y
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  // Next, compute y'*A*x
  ComputeDotProduct(nrow, x_ncol, x_ncol, xNorm2, t4, A.isDotProductOptimized);
  ierr = ComputeSPMV(A, x_ncol, z_ncol);  // b_computed = A*x_overlap
  if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
  double ytAx = 0.0;
  ierr        = ComputeDotProduct(nrow, y_ncol, z_ncol, ytAx, t4,
                                  A.isDotProductOptimized);  // y'*A*x
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  testsymmetry_data.depsym_spmv =
      std::fabs((long double)(xtAy - ytAx)) /
      ((xNorm2 * ANorm * yNorm2 + yNorm2 * ANorm * xNorm2) * (DBL_EPSILON));
  if (testsymmetry_data.depsym_spmv > 1.0)
    ++testsymmetry_data.count_fail;  // If the difference is > 1, count it wrong
  if (A.geom->rank == 0)
    HPCG_fout
        << "Departure from symmetry (scaled) for SpMV abs(x'*A*y - y'*A*x) = "
        << testsymmetry_data.depsym_spmv << endl;

    // Test symmetry of multi-grid
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(yncolopt->dev, yncolopt->host);
#endif
  // Compute x'*Minv*y
  ierr = ComputeMG(A, y_ncol, z_ncol);  // z_ncol = Minv*y_ncol
  if (ierr) HPCG_fout << "Error in call to MG: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(zncolopt->host, zncolopt->dev);
#endif
  double xtMinvy = 0.0;
  ierr           = ComputeDotProduct(nrow, x_ncol, z_ncol, xtMinvy, t4,
                                     A.isDotProductOptimized);  // x'*Minv*y
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(xncolopt->dev, xncolopt->host);
#endif
  // Next, compute z'*Minv*x
  ierr = ComputeMG(A, x_ncol, z_ncol);  // z_ncol = Minv*x_ncol
  if (ierr) HPCG_fout << "Error in call to MG: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(zncolopt->host, zncolopt->dev);
#endif
  double ytMinvx = 0.0;
  ierr           = ComputeDotProduct(nrow, y_ncol, z_ncol, ytMinvx, t4,
                                     A.isDotProductOptimized);  // y'*Minv*x
  if (ierr) HPCG_fout << "Error in call to dot: " << ierr << ".\n" << endl;

  testsymmetry_data.depsym_mg =
      std::fabs((long double)(xtMinvy - ytMinvx)) /
      ((xNorm2 * ANorm * yNorm2 + yNorm2 * ANorm * xNorm2) * (DBL_EPSILON));
  if (testsymmetry_data.depsym_mg > 1.0)
    ++testsymmetry_data.count_fail;  // If the difference is > 1, count it wrong
  if (A.geom->rank == 0)
    HPCG_fout << "Departure from symmetry (scaled) for MG abs(x'*Minv*y - "
                 "y'*Minv*x) = "
              << testsymmetry_data.depsym_mg << endl;

#ifdef HPCG_WITH_MORPHEUS
  Morpheus::copy(xexactopt->host, xncolopt->host, 0, xexactopt->host.size());
#else
  CopyVector(xexact, x_ncol);  // Copy exact answer into overlap vector
#endif

  int numberOfCalls = 2;
  double residual   = 0.0;
  for (int i = 0; i < numberOfCalls; ++i) {
    ierr = ComputeSPMV(A, x_ncol, z_ncol);  // b_computed = A*x_overlap
    if (ierr) HPCG_fout << "Error in call to SpMV: " << ierr << ".\n" << endl;
#ifdef HPCG_WITH_MORPHEUS
    // Needed for computing residual on host
    Morpheus::copy(zncolopt->dev, zncolopt->host);
#endif
    if ((ierr = ComputeResidual(A.localNumberOfRows, b, z_ncol, residual)))
      HPCG_fout << "Error in call to compute_residual: " << ierr << ".\n"
                << endl;
    if (A.geom->rank == 0)
      HPCG_fout << "SpMV call [" << i << "] Residual [" << residual << "]"
                << endl;
  }
  DeleteVector(x_ncol);
  DeleteVector(y_ncol);
  DeleteVector(z_ncol);

  return 0;
}
