/*
  In this header file, I have defined a simplified function call to
  the ARPACK solver routine for a complex, asymmetric eigenvalue
  problem Av = lv.  The looping procedure and final extraction of
  eigenvalues and vectors is handled automatically.  Most of the
  parameters to the FORTRAN functions are hidden from the user, since
  most of them are determined from user input anyway.
  
  The remaining parameters to the function calls are as follows:
  
    dsaupd(int n, int nev, complex<double> *Evals)
    dsaupd(int n, int nev, complex<double> *Evals, double **Evecs)

    n: the order of the square matrix A
    nev: the number of eigenvalues to be found, starting at the
         bottom.  Note that the highest eigenvalues, or some
	 other choices, can be found.  For now, the choice of
	 the lowest nev eigenvalues is hard-coded.
    Evals: a one-dimensional array of length nev to hold the
           eigenvalues.
    Evecs: a two-dimensional array of size nev by n to hold the
           eigenvectors.  If this argument is not provided, the
	   eigenvectors are not calculated.  Note that the
	   vectors are stored as columns of Evecs, so that the
	   elements of vector i are the values Evecs[j][i].

  The function is overdefined, so that you can call it in either
  fashion and the appropriate version will be run.

  To use these function calls, there must be a function
  defined in the calling program of the form

    av(int n, complex<double> *in, complex<double> *out)

  where the function av finds out = A.in, and n is the order of the
  matrix A.  This function must be defined before the statement that
  includes this header file, as it needs to know about this function.
  It is used in the looping procedure.

  Note that 0 < nev < n-1.

  Scot Shaw
  30 August 1999 */

#include <math.h>
#include <complex.h>

extern "C" void znaupd_(int *ido, char *bmat, int *n, char *which,
			int *nev, double *tol, complex<double> *resid,
			int *ncv, complex<double> *v, int *ldv,
			int *iparam, int *ipntr, complex<double> *workd,
			complex<double> *workl, int *lworkl,
			double *rwork, int *info);

extern "C" void zneupd_(int *rvec, char *All, int *select,
			complex<double> *d, complex<double> *v, int *ldv,
			double *sigma, complex<double> *workev, char *bmat,
			int *n, char *which, int *nev, double *tol,
			complex<double> *resid, int *ncv,
			complex<double> *v, int *ldv, int *iparam,
			int *ipntr, complex<double> *workd,
			complex<double> *workl, int *lworkl,
			double *rwork, int *ierr);

void znaupd(int n, int nev, complex<double> *Evals)
{
  int ido = 0; /* Initialization of the reverse communication
		  parameter. */

  char bmat[2] = "I"; /* Specifies that the right hand side matrix
			 should be the identity matrix; this makes
			 the problem a standard eigenvalue problem.
			 Setting bmat = "G" would have us solve the
			 problem Av = lBv (this would involve using
			 some other programs from BLAS, however). */

  char which[3] = "SM"; /* Ask for the nev eigenvalues of smallest
			   magnitude.  The possible options are
			   LM: largest magnitude
			   SM: smallest magnitude
			   LA: largest real component
			   SA: smallest real compoent
			   LI: largest imaginary component
			   SI: smallest imaginary component */

  double tol = 0.0; /* Sets the tolerance; tol<=0 specifies 
		       machine precision */

  complex<double> *resid;
  resid = new complex<double>[n];

  int ncv = 4*nev; /* The largest number of basis vectors that will
		      be used in the Implicitly Restarted Arnoldi
		      Process.  Work per major iteration is
		      proportional to N*NCV*NCV. */
  if (ncv>n) ncv = n;

  complex<double> *v;
  int ldv = n;
  v = new complex<double>[ldv*ncv];

  int *iparam;
  iparam = new int[11]; /* An array used to pass information to the routines
			   about their functional modes. */
  iparam[0] = 1;   // Specifies the shift strategy (1->exact)
  iparam[2] = 3*n; // Maximum number of iterations
  iparam[6] = 1;   /* Sets the mode of dsaupd.
		      1 is exact shifting,
		      2 is user-supplied shifts,
		      3 is shift-invert mode,
		      4 is buckling mode,
		      5 is Cayley mode. */

  int *ipntr;
  ipntr = new int[14]; /* Indicates the locations in the work array workd
			  where the input and output vectors in the
			  callback routine are located. */

  complex<double> *workd;
  workd = new complex<double>[3*n];

  int lworkl = 3*ncv*ncv + 5*ncv; /* Length of the workl array */
  complex<double> *workl;
  workl = new complex<double>[lworkl];

  double *rwork;
  rwork = new double[ncv];

  int info = 0; /* Passes convergence information out of the iteration
		   routine. */

  int rvec = 0; /* Specifies that eigenvectors should not be calculated */

  int *select;
  select = new int[ncv];
  complex<double> *d;
  d = new complex<double>[ncv]; /* This vector will return the
				     eigenvalues from the second routine,
				     dseupd. */
  double sigma;
  
  complex<double> *workev;
  workev = new complex<double>[3*ncv]; // I don't know what this is used for

  int ierr;

  /* Here we enter the main loop where the calculations are
     performed.  The communication parameter ido tells us when
     the desired tolerance is reached, and at that point we exit
     and extract the solutions. */

  do {
    znaupd_(&ido, bmat, &n, which, &nev, &tol, resid, 
	    &ncv, v, &ldv, iparam, ipntr, workd, workl,
	    &lworkl, rwork, &info);
    
    if ((ido==1)||(ido==-1)) av(n, workd+ipntr[0]-1, workd+ipntr[1]-1);
  } while ((ido==1)||(ido==-1));

  /* From those results, the eigenvalues and vectors are
     extracted. */

  if (info<0) {
         cout << "Error with znaupd, info = " << info << "\n";
         cout << "Check documentation in dsaupd\n\n";
  } else {
    zneupd_(&rvec, "All", select, d, v, &ldv, &sigma, workev,
	    bmat, &n, which, &nev, &tol, resid, &ncv, v, &ldv,
	    iparam, ipntr, workd, workl, &lworkl, rwork, &ierr);

    if (ierr!=0) {
      cout << "Error with zneupd, info = " << ierr << "\n";
      cout << "Check the documentation of zneupd.\n\n";
    } else if (info==1) {
      cout << "Maximum number of iterations reached.\n\n";
    } else if (info==3) {
      cout << "No shifts could be applied during implicit\n";
      cout << "Arnoldi update, try increasing NCV.\n\n";
    }

    /* Before exiting, we copy the solution information over to
       the arrays of the calling program, then clean up the
       memory used by this routine.  For some reason, when I
       don't find the eigenvectors I need to reverse the order of
       the values. */

    int i;
    for (i=0; i<nev; i++) Evals[i] = d[nev-1-i];

    // Sort the energies by real part

    complex<double> temp;
    for (int i=0; i<nev; i++) {
      for (int j=i; j<nev; j++) {
	if (Evals[j].real() > Evals[i].real()) {
	  temp = Evals[j];
	  Evals[j] = Evals[i];
	  Evals[i] = temp;
	}
      }
    }
    
    delete resid;
    delete v;
    delete iparam;
    delete ipntr;
    delete workd;
    delete workl;
    delete select;
    delete d;
  }
}


void znaupd(int n, int nev, complex<double> *Evals, complex<double> **Evecs)
{
  int ido = 0;
  char bmat[2] = "I";
  char which[3] = "SM";
  double tol = 0.0;
  complex<double> *resid;
  resid = new complex<double>[n];
  int ncv = 4*nev;
  if (ncv>n) ncv = n;
  complex<double> *v;
  int ldv = n;
  v = new complex<double>[ldv*ncv];
  int *iparam;
  iparam = new int[11];
  iparam[0] = 1;
  iparam[2] = 3*n;
  iparam[6] = 1;
  int *ipntr;
  ipntr = new int[14];
  complex<double> *workd;
  workd = new complex<double>[3*n];
  int lworkl = 3*ncv*ncv + 5*ncv;
  complex<double> *workl;
  workl = new complex<double>[lworkl];
  double *rwork;
  rwork = new double[ncv];
  int info = 0;
  int rvec = 1;  // Changed from above
  int *select;
  select = new int[ncv];
  complex<double> *d;
  d = new complex<double>[ncv];
  double sigma;
  complex<double> *workev;
  workev = new complex<double>[3*ncv];
  int ierr;

  do {
    znaupd_(&ido, bmat, &n, which, &nev, &tol, resid, 
	    &ncv, v, &ldv, iparam, ipntr, workd, workl,
	    &lworkl, rwork, &info);
    
    if ((ido==1)||(ido==-1)) av(n, workd+ipntr[0]-1, workd+ipntr[1]-1);
  } while ((ido==1)||(ido==-1));

  if (info<0) {
         cout << "Error with znaupd, info = " << info << "\n";
         cout << "Check documentation in dsaupd\n\n";
  } else {
    zneupd_(&rvec, "All", select, d, v, &ldv, &sigma, workev,
	    bmat, &n, which, &nev, &tol, resid, &ncv, v, &ldv,
	    iparam, ipntr, workd, workl, &lworkl, rwork, &ierr);

    if (ierr!=0) {
      cout << "Error with zneupd, info = " << ierr << "\n";
      cout << "Check the documentation of zneupd.\n\n";
    } else if (info==1) {
      cout << "Maximum number of iterations reached.\n\n";
    } else if (info==3) {
      cout << "No shifts could be applied during implicit\n";
      cout << "Arnoldi update, try increasing NCV.\n\n";
    }
    
    int i, j;
    for (i=0; i<nev; i++) Evals[i] = d[i];
    for (i=0; i<nev; i++) for (j=0; j<n; j++) Evecs[j][i] = v[i*n+j];

    complex<double> temp;
    for (int i=0; i<nev; i++) {
      for (int j=i; j<nev; j++) {
	if (Evals[j].real() > Evals[i].real()) {
	  temp = Evals[j];
	  Evals[j] = Evals[i];
	  Evals[i] = temp;
	  for (int k=0; k<n; k++) {
	    temp = Evecs[k][i];
	    Evecs[k][i] = Evecs[k][j];
	    Evecs[k][j] = temp;
	  }
	}
      }
    }
    
    delete resid;
    delete v;
    delete iparam;
    delete ipntr;
    delete workd;
    delete workl;
    delete select;
    delete d;
  }
}