// -*- c++ -*-
//
// $COPYRIGHT$
//

#ifndef ITL_MODIFIED_GRAM_SCHMIDT_ORTHOGONALIZATION_H
#define ITL_MODIFIED_GRAM_SCHMIDT_ORTHOGONALIZATION_H

#include <vector>

#include <itl/itl.h>
#include <itl/itl_tags.h>

namespace itl {

  template <class Vec>
  class modified_gram_schmidt {
  public:
    typedef typename itl_traits<Vec>::size_type size_type; 
    typedef typename itl_traits<Vec>::value_type value_type;

    template <class Size>
    modified_gram_schmidt(int restart, const Size& s) 
      : V(restart+1) {
      typedef typename itl_traits<Vec>::vector_category VecCat;
      do_initialize(VecCat(), s);
    }
    
    inline const Vec& operator[](size_type i) const {
      return V[i];
    }

    inline Vec& operator[](size_type i) {
      return V[i];
    }
  protected:
    //such as mtl::dense1D
    template <class Size>
    void do_initialize(referencing_object_tag, const Size& s) {
      for (size_type i=0; i<V.size(); ++i)
	V[i] = Vec(s);
    }
    
    //such as blitz array
    template <class Size>
    void do_initialize(non_referencing_object_tag, const Size& s) {
      for (size_type i=0; i<V.size(); ++i)
	itl::resize(V[i], s);
    }
    
    std::vector<Vec> V;
  };
  


  template <class Vec, class VecHi, class Size>
  void orthogonalize(modified_gram_schmidt<Vec>& V,  const VecHi& _Hi, Size i)
 {
    VecHi& Hi = const_cast<VecHi&>(_Hi);
    ++i;
    for (Size k = 0; k < i; k++) {
      Hi[k] = itl::dot_conj(V[i], V[k]);
      itl::add(itl::scaled(V[k], -Hi[k]), V[i]);
    }
  }
  
  template <class Vec, class VecS, class VecX, class Size>
  void combine(modified_gram_schmidt<Vec>& V, const VecS& s, VecX& x, Size i)
  {
    itl::copy(itl::scaled(V[0], s[0]), x);
    for (Size j = 1; j < i; ++j)
      itl::add(itl::scaled(V[j], s[j]), x);
  }
  
}

#endif