added exercise 4

exercises/ex04_skeleton/build_alps.txt

+# INSTRUCTIONS
+# ============
+#
+# 0. Make sure that lapack and blas are available
+#    (e.g. apt-get install liblapack-dev libblas-dev)
+# 1. Download http://alps.comp-phys.org/static/software/releases/alps-2.2.b1-r7195-src-with-boost.tar.gz
+#    (cf. http://alps.comp-phys.org/mediawiki/index.php/Download_and_install_ALPS_2)
+# 2. Unzip it (you do NOT need to build or install anything!)
+# 3. Set the following environment variable to the path of the unzipped directory:
+
+export ALPSDIR=.../alps-2.2.b1-r7195-src-with-boost
+
+# Now you can compile the skeleton by using:
+
+c++ -I${ALPSDIR}/boost -I${ALPSDIR}/alps/src -o skeleton skeleton.cpp -llapack -lblas
+
+# It will fail linking, because you have not yet implemented the Hamiltonian :-)
+#
+# Happy coding!
+#
+# If you run into problems that you can't solve, please contact your teaching assistant!

exercises/ex04_skeleton/skeleton.cpp

+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <cassert>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/random.hpp>
+#include <boost/tuple/tuple.hpp>
+
+class Hamiltonian
+{
+public:
+  typedef double Scalar;
+  typedef boost::numeric::ublas::vector<Scalar> Vector;
+
+  // Dimension of Vector Space
+  unsigned dimension() const;
+
+  // Calculate y = H x
+  void multiply(const Vector& x, Vector& y) const;
+
+  // ...
+};
+
+namespace ietl
+{
+	/* here we include a "mult" function in the ietl name space specifically designed for use with our Hamiltonian class objects.  */
+  inline void mult( const Hamiltonian& h, const Hamiltonian::Vector& x, Hamiltonian::Vector& y )
+  {
+    h.multiply( x, y );
+	/*To minimize your memory requirements you could implement this method "on the fly", i.e. in a way that never constructs the hamitonian matrix explicitly. */
+  }
+}
+
+// ietl::mult() needs to be declared before including these
+#include <ietl/interface/ublas.h>
+#include <ietl/lanczos.h>
+
+std::pair< std::vector<double>, std::vector<int> >
+diagonalize( const Hamiltonian& h, unsigned nvals=1, unsigned maxiter=1000)
+{
+  typedef ietl::vectorspace<Hamiltonian::Vector> Vectorspace;
+  typedef ietl::lanczos<Hamiltonian,Vectorspace> Lanczos;
+  typedef ietl::lanczos_iteration_nlowest<double> Iteration;
+  typedef boost::mt19937 Generator; 
+
+  Vectorspace vspace(h.dimension());
+  Lanczos solver(h,vspace);
+  Iteration iter(maxiter,nvals); /* maxiter limits the number of Lanczos iterations and nvals tells the solver the number of eigenvalues to calculate.*/
+  solver.calculate_eigenvalues(iter,Generator()); /* Generator() creates a pseudo-random number generator (Mersenne-Twister) used for starting the Lanczos algorithm. */
+
+  if( iter.iterations() == maxiter )
+    std::cerr << "Lanczos did not converge in " << iter.iterations() << " iterations." << std::endl;
+  else
+    std::cout << "  Lanczos converged after " << iter.iterations() << " iterations." << std::endl;
+
+  return std::make_pair(solver.eigenvalues(),solver.multiplicities());
+}
+

exercises/ex04_skeleton/skeleton.py

+import math as m, numpy as np
+
+from scipy.sparse.linalg as la
+# Documentation at http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html
+
+d = # size of local basis
+L = # length of the chain
+
+
+def heisenberg(v):
+"""
+Given an input vector v, this function returns H*v.
+"""
+    pass
+
+def sparse_diag():
+	D = d**L
+	H = la.LinearOperator( (D,D), matvec=heisenberg, dtype='D')
+	(w,v) = la.eigsh(H, 10, which='SA')
+	return w
+