#ifdef LAPACK
#include "../include/solve.hpp"

namespace ngsolve
{
  using namespace ngsolve;


  /* ***************************** Numproc EVP-Lapack *************************** */

  ///
  class NumProcEVPLapack : public NumProc
  {
  protected:
    BilinearForm * bfa;
    BilinearForm * bfm;

    GridFunction * gfu;
    int num;

    double prec, shift;
    bool print;
  public:
    NumProcEVPLapack (PDE & apde, const Flags & flags);
    virtual ~NumProcEVPLapack();

    static NumProc * Create (PDE & pde, const Flags & flags)
    {
      return new NumProcEVPLapack (pde, flags);
    }

    virtual void Do();

    virtual string GetClassName () const
    {
      return " Eigenvalue Problem (Lapack)";
    }

    virtual void PrintReport (ostream & ost)
    {
      ost << GetClassName() << endl
	  << "Bilinear-form A = " << bfa->GetName() << endl
	  << "Bilinear-form M = " << bfm->GetName() << endl
	  << "Gridfunction  = " << gfu->GetName() << endl;
    }
  };


  NumProcEVPLapack :: NumProcEVPLapack (PDE & apde, const Flags & flags)
    : NumProc (apde)
  {
    bfa = pde.GetBilinearForm (flags.GetStringFlag ("bilinearforma", ""));
    bfm = pde.GetBilinearForm (flags.GetStringFlag ("bilinearformm", ""));
    gfu = pde.GetGridFunction (flags.GetStringFlag ("gridfunction", ""));
    num = int(flags.GetNumFlag ("num", 500));
    shift = flags.GetNumFlag ("shift",-1.); 
  }

  NumProcEVPLapack :: ~NumProcEVPLapack()
  {
    ;
  }


  void NumProcEVPLapack :: Do()
  {
    cout << "solve evp with Lapack" << endl;

    int dim = bfa->GetFESpace().GetDimension();
    int size = bfa->GetMatrix().Height();
    bool iscomplex = bfa->GetFESpace().IsComplex();

    
    if (!iscomplex)
      {
	int n = dim*size;
	VVector<> hv(n), hva(n), hvm(n);
	Matrix<> mata(n), matm(n);
	
	for (int i = 0; i < n; i++)
	  {
	    hv = 0.0;
	    hv(i) = 1.0;
	    hva = bfa->GetMatrix() * hv;
	    hvm = bfm->GetMatrix() * hv;

	    for (int j = 0; j < n; j++)
	      {
		mata(j,i) = hva(j);
		matm(j,i) = hvm(j);
	      }
	  }

	Vector<> lami(n);
	Matrix<> evecs(n);

	LapackGHEPEPairs (n, &mata(0,0), &matm(0,0), &lami(0));
	cout.precision(12);
	cout << "lami = " << lami << endl;

	// LaEigNSSolve (n, &mata(0,0), &matm(0,0), &lami(0), 1, &evecs(0,0), 0, 'N');

	cout << "eigensystem done" << endl;
	;
      }
    else
      {
	int n = dim*size;
	VVector<Complex> hv(n), hva(n), hvm(n), hwa(n), hwm(n), w(n);
	/*
	  Matrix<Complex> mata(n), matm(n);

	  for (int i = 0; i < n; i++)
	  {
	  hv = 0.0;
	  hv(i) = 1.0;
	  hva = bfa->GetMatrix() * hv;
	  hvm = bfm->GetMatrix() * hv;

	  for (int j = 0; j < n; j++)
	  {
	  mata(j,i) = hva(j);
	  matm(j,i) = hvm(j);
	  }
	  }

	  cout << "matrix copied" << endl;
	*/

	// Arnoldi:
	int m = min2 (num, n);

	Matrix<Complex> matV (m,n), matVt(n,m);
	Matrix<Complex> matH(m), matHt(m), matI(m);
	
	const BaseSparseMatrix&  mata1 =
	  dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix());
	
	BaseSparseMatrix&  mata =
	  dynamic_cast<BaseSparseMatrix&> (*mata1.CreateMatrix());
	
	const BaseSparseMatrix& matm = 
	  dynamic_cast<const BaseSparseMatrix&> (bfm->GetMatrix()); 
	
	mata.AsVector() = mata1.AsVector() - shift*matm.AsVector();  
	BaseMatrix * inva = mata.InverseMatrix();	    
	
	// BaseMatrix * inva = dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix()).InverseMatrix();

	hv.SetRandom();
	matV = Complex(0.0);
	matH = Complex(0.0);
	matI = Complex(0.0);
	for (int i = 0; i < m; i++)
	  matI(i,i) = 1.0;

	Complex len = sqrt (S_InnerProduct<Complex> (hv, hv));
	hv /= len;


	for (int i = 0; i < m; i++)
	  {
	    cout << "i = " << i << endl;
	    for (int j = 0; j < n; j++)
	      matV(i,j) = hv(j);
	    
	    hva = bfm->GetMatrix() * hv;
	    hvm = *inva * hva;

	    for (int j = 0; j <= i; j++)
	      {
		Complex sum = 0.0;
		for (int k = 0; k < n; k++)
		  sum += hvm(k) * matV(j,k);
		matH(j,i) = sum;
		for (int k = 0; k < n; k++)
		  hvm(k) -= sum * matV(j,k);
	      }
	    
	    hv = hvm;
	    Complex len = sqrt (S_InnerProduct<Complex> (hv, hv));
	    if (i<m-1)
	      matH(i+1,i) = len; 

	    //  cout << ", len = " << len << endl;
	    hv /= len;
	  }

	
	//	(*testout) << "matH = " << endl << matH << endl;
	//	(*testout) << "V^T V = " << endl << matV * Trans (matV) << endl;
	cout << "has Hessenberg" << endl;


	Vector<Complex> lami(m);
	Matrix<Complex> evecs(m);
	Matrix<Complex> levecs(n,m);
	matHt = Trans (matH);

	evecs = Complex (0.0);
	lami = Complex (0.0);

	// LaEigNSSolve (mata.Height(), &mata(0,0), &matm(0,0), &lami(0), 1, &evecs(0,0), 0, 'N');
	// LaEigNSSolve (matH.Height(), &matHt(0,0), &matI(0,0), &lami(0), 1, &evecs(0,0), 0, 'N');

	LapackHessenbergEP (matH.Height(), &matHt(0,0), &lami(0), &evecs(0,0));

	cout << "hessenberg is back" << endl;

	// levecs = Trans(matV) * Trans(evecs);
	matVt = Trans (matV);
	levecs = Trans(evecs * Trans (matVt));
	
	for (int i = 0; i < m; i++)
	  lami(i) =  1.0 / lami(i) + shift;

	cout << "eigensystem done" << endl;


	// sort by angle
	BitArray used(m);
	ARRAY<int> reorder;
	used.Clear();
	for (int i = 0; i < m; i++)
	  {
	    int mini = -1;
	    double val = 0;
	    for (int j = 0; j < m; j++)
	      if (!used[j] && lami(j).real() < 1000 && lami(j) != Complex(100,100))
		{
		  if (mini < 0 ||
		      lami(j).imag() / lami(j).real() > val)
		    {
		      mini = j;
		      val = lami(j).imag() / lami(j).real();
		    }
		}
	    if (mini != -1)
	      {
		reorder.Append (mini);
		used.Set(mini);
	      }
	  }
	for (int i = 0; i < m; i++)
	  {
	    if (!used[i] && lami(i) != Complex(100,100))
	      reorder.Append(i);
	  }

	for (int i = 0; i < reorder.Size(); i++)
	  cout << "lam(" << i << ") = " 
	       << lami(reorder[i]).real()/M_PI/M_PI << "+i" << lami(reorder[i]).imag()/M_PI/M_PI << endl;
	//	cout << "lami = " << endl << lami << endl;
	//	cout << "reorder = " << reorder << endl;
	ofstream out ("eigenvalues.dat");
	//	for (int i = 0; i < reorder.Size(); i++)
	//	  out << i << " " << lami(reorder[i]).real() << " " << lami(reorder[i]).imag() << endl;
	
	/*
	  (*testout) << "mata = " << bfa->GetMatrix() << endl
	  << " = " << endl << mata << endl;
	*/

	for (int i = 0; i < min2 (gfu->GetMultiDim(), reorder.Size()); i++)
	  {
	    FlatVector<Complex> vu = gfu->GetVector(i).FVComplex();
	    for (int j = 0; j < n; j++)
	      vu(j) = levecs(j, reorder[i]);
	  }
	
	int postit = 1; 
	// postiterate:
	Mat<2,2,Complex> sma, smm;
	Vec<2,Complex> smlam;
	for (int i = 0; i < min2 (gfu->GetMultiDim(), reorder.Size()); i++)
	  {
	    for (int k = 0; k < postit; k++)
	      {
		//	cout << "lam = " << lami(reorder[i]) << endl;
	       	BaseVector & vu = gfu->GetVector(i);
		//SZ hva = bfa->GetMatrix() * vu;
		hva = mata * vu;
	       	hvm = bfm->GetMatrix() * vu;
		hv = hva - ( - shift +lami(reorder[i])) * hvm;
		// cout << "def = " << hv.L2Norm() << endl;
		w = (*inva) * hv;
		
		//SZ	hwa = bfa->GetMatrix() * w;
		hwa = mata * w;
		hwm = bfm->GetMatrix() * w;

		sma(0,0) = S_InnerProduct<Complex> (hva, vu);
		sma(0,1) = sma(1,0) = S_InnerProduct<Complex> (hva, w);
		sma(0,0) = S_InnerProduct<Complex> (hwa, w);

		smm(0,0) = S_InnerProduct<Complex> (hvm, vu);
		smm(0,1) = sma(1,0) = S_InnerProduct<Complex> (hvm, w);
		smm(0,0) = S_InnerProduct<Complex> (hwm, w);

		LaEigNSSolve(2, &sma(0,0), &smm(0,0), &smlam(0), 0, 0, 0, 'N');
		//	cout << "smlami = " << smlam << endl;
		
		out << i << " " << lami(reorder[i]).real() << " " << lami(reorder[i]).imag() 
		<< " " << smlam(0).real()-1 << " " << smlam(0).imag() << endl;
	      }
	  }
	
	
      }
    
  }




  


  namespace
  {
    class Init
    { 
    public: 
      Init ();
    };
    
    Init::Init()
    {
      GetNumProcs().AddNumProc ("evplapack", NumProcEVPLapack::Create, NumProcEVPLapack::PrintDoc);
    }
    
    
    Init init;
    
  }
  


}
#endif