#include <slepceps.h>
#include <chrono>

void
readMatA(Mat &A)
{
  char file[PETSC_MAX_PATH_LEN] = "";
  PetscBool flag;
  PetscOptionsGetString(NULL,NULL,"-f", file, sizeof(file), &flag);

  PetscPrintf(PETSC_COMM_WORLD, "\n Reading Matrix A from file: %s", file);

  PetscViewer viewer;
  PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);
  MatLoad(A,viewer);
  PetscViewerDestroy(&viewer);

  PetscInt size;
  MatGetSize(A, &size, NULL);
  PetscPrintf(PETSC_COMM_WORLD, "\n Read Matrix A of size = %d", size);

  PetscReal normA;
  MatNorm(A, NORM_INFINITY, &normA);
  PetscPrintf(PETSC_COMM_WORLD,"\n Matrix A norm = %e", normA);
}

int main(int argc, char **argv)
{
  SlepcInitialize(&argc,&argv,(char*)0,NULL);

  Mat A;
  EPS eps;
  ST st;
  KSP ksp;
  PC pc;
  PetscInt nev = 5;
  PetscOptionsGetInt(NULL,NULL,"-nev", &nev, NULL);

  MatCreate(PETSC_COMM_WORLD,&A);
  readMatA(A);

  std::chrono::high_resolution_clock::time_point ts, te;
  std::chrono::duration<double> time;
  ts = std::chrono::high_resolution_clock::now();

  PetscPrintf(PETSC_COMM_WORLD,"\n Computating Eigen Frequencies...\n");

  EPSCreate(PETSC_COMM_WORLD,&eps);
  EPSSetOperators(eps,A,NULL);
  EPSSetProblemType(eps,EPS_NHEP);

  EPSGetST(eps,&st);
  STSetType(st,STSINVERT);

  STGetKSP(st,&ksp);
  KSPSetType(ksp,KSPPREONLY);

  KSPGetPC(ksp,&pc);
  PCSetType(pc,PCLU);
  PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
  STGetOperator(st,NULL);
  PCFactorSetUpMatSolverType(pc);

  Mat F;
  PCFactorGetMatrix(pc,&F);
  MatMumpsSetCntl(F,3,1e-8);
  MatMumpsSetCntl(F,4,0);
  MatMumpsSetIcntl(F,13,1);

  EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE);
  EPSSetTarget(eps,0.0);

  EPSSetDimensions(eps, nev, PETSC_DEFAULT, PETSC_DEFAULT);
  EPSSetTolerances(eps, 1e-6, 100);

  EPSSetFromOptions(eps);
  EPSSolve(eps);

  PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
  EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
  PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

  EPSDestroy(&eps);
  MatDestroy(&A);

  te = std::chrono::high_resolution_clock::now();
  time = te - ts;

  PetscPrintf(PETSC_COMM_WORLD, "\n Computation finished in %f seconds\n", time.count());

  SlepcFinalize();
  return 0;
}