//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: my_interface.h,v $  $Revision: 1.8 $  $State: Exp $ 
//     $Author: llee $  $Date: 2001/10/23 16:15:01 $ 
//     $Locker:  $ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//              
// This is to provide an example interface to where matrix is implemented 
// as a vector of vector (column-wise) and vector is an stl vector.
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
// Copyright 1997, University of Notre Dame.
// Authors: Andrew Lumsdaine, Lie-Quan Lee
//
// This file is part of the Iterative Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software;  see the
// file LICENSE.  If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN  46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//---------------------------------------------------------------------------
//                                                                           
// REVISION HISTORY:                                                         
//                                                                           
// $Log: my_interface.h,v $
// Revision 1.8  2001/10/23 16:15:01  llee
// *** empty log message ***
//
// Revision 1.7  2001/10/18 16:57:11  llee
// *** empty log message ***
//
// Revision 1.6  2000/07/27 15:54:43  llee1
// *** empty log message ***
//
// Revision 1.5  2000/07/26 21:34:12  llee1
// *** empty log message ***
//
// Revision 1.4  2000/07/19 19:58:00  llee1
// *** empty log message ***
//
// Revision 1.3  2000/07/18 18:00:32  llee1
// *** empty log message ***
//
// Revision 1.2  2000/07/18 14:30:44  llee1
// *** empty log message ***
//
// Revision 1.1  2000/07/17 15:44:06  llee1
// *** empty log message ***
//
//
//===========================================================================
#ifndef MY_INTERFACE_H
#define MY_INTERFACE_H



#include <vector>
#include <complex>
#include <algorithm>
#include "scaled1D.h"
#include <itl/itl_tags.h>
#include <itl/number_traits.h>

namespace std {

  inline float conj(float a) {
    return a; 
  }

  inline double conj(double a) {
    return a; 
  }

  inline int conj(int a) {
    return a; 
  }

}


namespace itl {
  
  //: The vector type used inside of the ITL routines for work space
  template <class Vec>
  struct itl_traits {
    typedef non_referencing_object_tag vector_category;
    typedef typename Vec::value_type value_type;
    typedef typename Vec::size_type size_type;
  };


  template<class T>
  inline typename std::vector<std::vector<T> >::size_type 
  nrows(const std::vector<std::vector<T> >& A)
  {
    return A.size();
  }

  template<class T>
  inline typename std::vector<std::vector<T> >::size_type 
  ncols(const std::vector<std::vector<T> >& A)
  {
    return A[0].size();
  }

  template<class Matrix>
  inline typename Matrix::size_type
  nrows(const Matrix& A) {
    return A.nrows();
  }

  template<class Matrix>
  inline typename Matrix::size_type
  ncols(const Matrix& A) {
    return A.ncols();
  }


  template <class Vec>
  inline typename number_traits< typename Vec::value_type >::magnitude_type
  two_norm(const Vec& v) {
    typename number_traits< typename Vec::value_type >::magnitude_type d=0;
    for ( typename Vec::const_iterator i = v.begin(); i != v.end(); ++i )
      d += *i * *i;
    return std::sqrt(d);
  }
 
  template <class VecA, class VecB>
  inline void copy(const VecA& a, VecB& b) {
    std::copy(a.begin(), a.end(), b.begin());
  }
  
  template <class Matrix, class VecX, class VecY>
  inline void mult(const Matrix& A, const VecX& x, VecY& y) {
    typedef typename Matrix::size_type size_type;
    size_type s = nrows(A);
    for (size_type i=0; i<s; i++) {
      typename VecX::value_type tmp = 0;
      for ( size_type j=0; j<ncols(A); j++)
	tmp += A[j][i] * x[j];
      y[i] = tmp;
    }

  }

  template <class Matrix, class VecX, class VecY>
  inline void trans_mult(const Matrix& A, const VecX& x, VecY& y) {
    typedef typename Matrix::size_type size_type;
    for ( size_type i=0; i<ncols(A); i++) {
      typename VecX::value_type tmp = 0;
      for (size_type j=0; j<nrows(A); j++)
	tmp += A[i][j] * x[j];
      y[i] = tmp;
    }

  }

  template <class Matrix, class VecX, class VecY, class VecR>
  inline void mult(const Matrix& A, const VecX& x, const VecY& y, VecR& r) {
    typedef typename Matrix::size_type size_type;
    size_type s = nrows(A);
    for (size_type i=0; i<s; i++) {
      typename VecX::value_type tmp = 0;
      for (size_type j=0; j<ncols(A); j++)
	tmp += A[j][i] * x[j];
      r[i] = tmp + y[i];
    }

  }
  
  template <class VecA, class VecB>
  inline typename VecA::value_type dot(const VecA& a, const VecB& b) {
    typename VecA::value_type tmp = 0;

    typename VecA::const_iterator a_i = a.begin();
    typename VecB::const_iterator b_i = b.begin();
    while ( a_i != a.end() ) {
      tmp += *a_i * *b_i;
      ++a_i;
      ++b_i;
    }
    return tmp;
  }

  template <class VecA, class VecB>
  inline typename VecA::value_type dot_conj(const VecA& a, const VecB& b) {
    typename VecA::value_type tmp = 0;

    typename VecA::const_iterator a_i = a.begin();
    typename VecB::const_iterator b_i = b.begin();
    while ( a_i != a.end() ) {
      tmp += *a_i * std::conj(*b_i);
      ++a_i;
      ++b_i;
    }
    return tmp;
  }

  template <class VecX, class VecY, class VecZ>
  inline void add(const VecX& x, const VecY& y, VecZ& z) {
    typename VecX::const_iterator x_i = x.begin();
    typename VecY::const_iterator y_i = y.begin();
    typename VecZ::iterator z_i = z.begin();

    while ( x_i != x.end() ) {
      *z_i = *x_i + *y_i;
      ++x_i; ++y_i; ++z_i;
    }
  }

  template <class VecX, class VecY>
  inline void add(const VecX& x, VecY& y) {
    add(y, x, y);
  }

  template <class VecX, class VecY, class VecZ, class VecR>
  inline void add(const VecX& x, const VecY& y, const VecZ& z, VecR& r ) {
    typename VecX::const_iterator x_i = x.begin();
    typename VecY::const_iterator y_i = y.begin();
    typename VecZ::const_iterator z_i = z.begin();
    typename VecR::iterator r_i = r.begin();
    
    while ( x_i != x.end() ) {
      *r_i = *x_i + *y_i + *z_i;
      ++r_i; ++x_i; ++y_i; ++z_i;
    }
  }

  template <class Vec, class T>
  inline void scale(Vec& v, T t) {
    typename Vec::iterator v_i = v.begin();
    while ( v_i != v.end() ) {
      *v_i *= t;
      ++v_i;
    }
  }

  template <class Vec>
  inline typename Vec::size_type 
  size(const Vec& v) {
    return v.size();
  }

  template <class Vec, class Size>
  inline void
  resize(Vec& v, const Size& s) {
    return v.resize(s);
  }

  template <class T>
  class internal_matrix {
  public:
    typedef typename std::vector<std::vector<T> >::size_type size_type;
    typedef std::vector<T> OneD;
    internal_matrix(size_type m, size_type n)
      : A(std::vector<std::vector<T> >(n, std::vector<T>(m, 0))) {}
    
    std::vector<T>& operator[](size_type i) {
      return A[i];
    }
    const std::vector<T>& operator[](size_type i) const {
      return A[i];
    }

    T& operator()(size_type i, size_type j) {
      return A[j][i];
    }

    T operator()(size_type i, size_type j) const {
      return A[j][i];
    }
    
    size_type nrows() const { return A[0].size(); }
    size_type ncols() const { return A.size(); }

    struct submatrix {
      typedef typename std::vector<std::vector<T> >::size_type size_type;
      submatrix( std::vector<std::vector<T> >& _A, 
		 size_type _nrows, size_type _ncols) 
	: A(_A), nrows_(_nrows), ncols_(_ncols) {}
      std::vector<std::vector<T> >& A;
      size_type nrows() const { return nrows_;; }
      size_type ncols() const { return ncols_; }

      std::vector<T>& operator[](size_type i) {
	return A[i];
      }

      const std::vector<T>& operator[](size_type i) const {
	return A[i];
      }

      T& operator()(size_type i, size_type j) {
	return A[j][i];
      }
      T operator()(size_type i, size_type j) const {
	return A[j][i];
      }
  
      size_type nrows_;
      size_type ncols_;
    };

    submatrix 
    sub_matrix(size_type, size_type en, size_type, size_type em) {
      return submatrix(A, en, em);
    }

    submatrix 
    sub_matrix(size_type, size_type en, size_type, size_type em) const {
      return submatrix(const_cast<std::vector<std::vector<T> >&>(A), en, em);
    }

  private:
    std::vector<std::vector<T> > A;
  };

  template <class T>
  struct internal_matrix_traits {
    typedef internal_matrix<T> Matrix;
  };

  template <class Hessenberg, class Vec>
  void upper_tri_solve(const Hessenberg& hh, Vec& rs, int i) {
    i--;
    typedef typename Vec::value_type T;
    rs[i] /= hh(i, i);
    for (int ii = 1; ii <= i; ii++) {
      int k  = i - ii;
      int k1 = k + 1;
      T t  = rs[k];
      for (int j = k1; j <= i; j++) 
	t -= hh(k, j) * rs[j];
      rs[k] = t / hh(k, k);
    }
  }

}

#endif