/***************************************************************************
 * $Id: schrodinger1.cpp,v 1.7 2004/02/21 09:51:20 troyer Exp $
 *
 * An example of the Lanczos method with Schrodinger wave equation and 
 * using stencil.
 *
 * Copyright (C) 2001-2003 by Prakash Dayal <prakash@comp-phys.org>
 *                            Matthias Troyer <troyer@comp-phys.org>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 **************************************************************************/

#include <ietl/lanczos.h>
#include <ietl/vectorspace.h>
#include <ietl/interface/blitz.h>
#include <boost/random/lagged_fibonacci.hpp>
#include <blitz/range.h>
#include <cmath>
#include <limits>
#include <iterator>
#include <iostream>

double m;
double grid_dl=0.5;

BZ_DECLARE_STENCIL2(schrodinger_stencil,y,x)
  y = -blitz::Laplacian2D(x)/(2*m)/grid_dl/grid_dl;
BZ_END_STENCIL

namespace ietl {
template< class S, class T, int N>
inline void mult (const S& s, const blitz::Array<T,N>& x, blitz::Array<T,N>& y) {
  y=0.;
  applyStencil(s,y,x);
} 

template<class VectorType>
class vectorspace2d {
public:
  typedef VectorType vector_type; 
  typedef typename  VectorType::T_numtype scalar_type;
  typedef int size_type;  
  vectorspace2d(size_type n, size_type m):n_(n), m_(m){}
  inline size_type vec_dimension() const {
    return n_ * m_;
  }
  void project(vector_type& vector) const {
    vector(blitz::Range(0,n_-1),0) = 0.;
    vector(blitz::Range(0,n_-1),m_-1) = 0.;
    vector(0,blitz::Range(0,m_-1)) = 0.;
    vector(n_-1,blitz::Range(0,m_-1)) = 0.;
  }
  vector_type new_vector()const {
    return VectorType(n_,m_);
  }  
private:
  size_type n_; 
  size_type m_;
};

}

int main() {
  const int N=4;
  m=1.;
  typedef blitz::Array<double,2> Vector;
  typedef ietl::vectorspace2d<Vector> Vecspace;
  typedef boost::lagged_fibonacci607 Gen;
  
  Vecspace vec(N,N);
  Gen mygen;
  schrodinger_stencil s;
  ietl::lanczos<schrodinger_stencil,Vecspace> mat(s,vec);

  // Creation of an iteration object:    
  int max_iter = 4*N*N;  
  double rel_tol = 500*std::numeric_limits<double>::epsilon();
  double abs_tol = std::pow(std::numeric_limits<double>::epsilon(),2./3);
  int n_lowest_eigenval = 3;
  ietl::lanczos_iteration_nlowest<double> 
        iter(max_iter, n_lowest_eigenval, rel_tol, abs_tol);
  std::vector<double> eigen;
  std::vector<double> err;
  std::vector<int> multiplicity;
  std::cout << "Computation of 3 lowest converged eigenvalues\n\n";
  std::cout << "-----------------------------------\n\n";  
  try {
    mat.calculate_eigenvalues(iter,mygen);
    eigen = mat.eigenvalues();
    err = mat.errors();
    multiplicity = mat.multiplicities();
  }
  catch (std::runtime_error& e) {
    std::cout << e.what() << "\n";
  } 
  std::cout << "#        eigenvalue            error         multiplicity\n";
  std::cout.precision(10);
  for (int i=0;i<eigen.size();++i) 
    std::cout << i << "\t" << eigen[i] << "\t" << err[i] << "\t" 
	      << multiplicity[i] << "\n";
  
  // call of eigenvectors function follows: 
  std::cout << "\nEigen vectors computations for 3 lowest eigenvalues:\n\n";  
  std::vector<double>::iterator start = eigen.begin();
  std::vector<double>::iterator end = eigen.begin()+3;
  std::vector<Vector> eigenvectors; // for storing the eigen vectors. 
  ietl::Info<double> info; // (m1, m2, ma, eigenvalue, residual, status).
  
  try {
    mat.eigenvectors(start, end, std::back_inserter(eigenvectors),info, mygen); 
  }
  catch (std::runtime_error& e) {
    std::cout << e.what() << std::endl;
  }  
  
  std::cout << "Printing eigen Vectors:\n\n"; 
  for(std::vector<Vector>::iterator it = eigenvectors.begin(); it != eigenvectors.end(); it++){
    std::cout << *it << "\n\n";
    //std::copy((it)->begin(),(it)->end(),std::ostream_iterator<double>(std::cout,"\n"));
    //std::cout << "\n\n";
  }
  std::cout << " Information about the eigenvector computations:\n\n";
  for(int i = 0; i < info.size(); i++)
    std::cout << " m1(" << i+1 << "): " << info.m1(i) << ", m2(" << i+1 << "): "
	      << info.m2(i) << ", ma(" << i+1 << "): " << info.ma(i) << " eigenvalue("
	      << i+1 << "): " << info.eigenvalue(i) << " residual(" << i+1 << "): "
	      << info.residual(i) << " error_info(" << i+1 << "): "
	      << info.error_info(i) <<"\n\n";
  return 0;
}