#include <stdio.h>
#include <stdlib.h>
#include <slepceps.h>

// Minimal example for EPS using ScaLAPACK
int main(int argc, char** argv){
  PetscInt N = 16000, nconv, i;
  Mat mat;
  EPS eps;
  PetscScalar eval_r;

  PetscCall(SlepcInitialize(&argc, &argv, NULL, NULL));


  // Generate random matrix
  PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, NULL, &mat));
  PetscCall(MatSetRandom(mat, NULL));


  // Set up and run EPS
  PetscCall(EPSCreate(PETSC_COMM_WORLD, &eps));
  PetscCall(EPSSetOperators(eps, mat, NULL));
  PetscCall(EPSSetType(eps, EPSSCALAPACK));
  PetscCall(EPSSetProblemType(eps, EPS_HEP));
  PetscCall(EPSSetWhichEigenpairs(eps, EPS_ALL));
  PetscCall(EPSSolve(eps));
  PetscCall(EPSGetConverged(eps, &nconv));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%d eigenpairs converged\n", nconv));
  PetscCall(MatDestroy(&mat));


  // Report a few eigenvalues
  for (i = 0; i < 5; i++){
    PetscCall(EPSGetEigenvalue(eps, i, &eval_r, NULL));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "lamda_%d = %e + %ei\n", i, eval_r));
  }


  PetscCall(EPSDestroy(&eps));
  PetscCall(SlepcFinalize());


  return(0);
  }