// -*- c++ -*-
//===========================================================================
//  CVS Information:                                                         
//                                                                           
//     $RCSfile: mtl_classical_gram_schmidt.h,v $  $Revision: 1.10 $  $State: Exp $ 
//     $Author: llee $  $Date: 2002/06/19 15:37:07 $ 
//     $Locker:  $ 
//---------------------------------------------------------------------------
//                                                                           
// DESCRIPTION                                                               
//
// This is a module for the last template in gmres function. It
// provides the basis stored in MTL dense column-wise matrix
// format. The class Gram-Schmidt orthogonalization with iterative
// refinement is performed.  Matrix-vector operations consist of the
// lowel level linear algebra operations.
//
//---------------------------------------------------------------------------
//                                                                           
// LICENSE AGREEMENT                                                         
//     $COPYRIGHT$
//---------------------------------------------------------------------------
//                                                                           
// REVISION HISTORY:                                                         
//                                                                           
// $Log: mtl_classical_gram_schmidt.h,v $
// Revision 1.10  2002/06/19 15:37:07  llee
// change the behavior of combine in OrthogonalBasis
//
// Revision 1.9  2002/06/15 15:49:47  llee
// *** empty log message ***
//
// Revision 1.8  2002/06/15 15:21:48  llee
// *** empty log message ***
//
// Revision 1.7  2002/06/02 04:33:36  llee
// *** empty log message ***
//
// Revision 1.2  2002/03/27 17:49:22  llee
// include mtl_amend.h
//
// Revision 1.1  2001/10/23 02:18:14  llee
// *** empty log message ***
//
// Revision 1.1  2001/10/18 14:08:32  llee
// re-organize the directory structures
//
// Revision 1.3  2001/07/05 22:28:58  llee1
// gcc 3.0 fix
//
// Revision 1.2  2000/07/25 22:57:54  llee1
// *** empty log message ***
//
//                                                                           
//===========================================================================


#ifndef ITL_CLASSICAL_GRAM_SCHMIDT_ORTHOGONALIZATION_H
#define ITL_CLASSICAL_GRAM_SCHMIDT_ORTHOGONALIZATION_H

#include <itl/itl.h>
#include <itl/itl_tags.h>
#include <mtl/matrix.h>
#include <itl/interface/detail/mtl_amend.h> //for linalg_traits of std::vector

#define SP_NAME(X) X##_
  void   SP_NAME(dgemv)(char* trans, int* m, int* n,  double*  alpha,
			double *da,  int* lda, double *dx, int* incx, 
			double* dbeta, double *dy, int* incy);


namespace itl {
  //in parallel envirnment, this is an all-reduce operation
  template<class VecX, class VecY>
  inline void all_reduce(const VecX& x, const VecY& yy, int k) {
    VecY& y = const_cast<VecY&>(yy);
    for (int i=0; i<k; i++)
      y[i] += x[i];
  }


  //This is to use mtl dense column-wise matrix for orthogonal vectors
  template <class Vec>
  class classical_gram_schmidt {
    typedef typename itl_traits<Vec>::value_type value_type;
    typedef typename mtl::matrix<value_type, mtl::rectangle<>,
				 mtl::dense<>,
				 mtl::column_major>::type ColumnMatrix;
  public:
    typedef size_t size_type;

    template <class Size>
    classical_gram_schmidt(int restart, const Size& s)
      : V(s, restart+1), ss(restart+1), w(s) {}
    
    
    typedef typename ColumnMatrix::OneD OneD;
    OneD operator[](size_type i) const {
      return V[i];
    }

    ColumnMatrix V;
    Vec ss;

    Vec w;
  };


    template <class Vec, class VecHi>
    void orthogonalize(classical_gram_schmidt<Vec>& KS,
		       const VecHi& _Hi, int i) {
      /*classical Gram-Schmidt with 1 iterative refinement.
      */
      VecHi& Hi = const_cast<VecHi&>(_Hi);
      char trans = 'T' ;
      char not_trans = 'N';
      int m = KS.V.nrows();
      double alpha = 1.0;
      double beta = 0.0;
      int one = 1;
      int i1 = i + 1; 
      

      for (int k=0; k<=i; k++) Hi[k] = 0.0;
      
      for (int ii=0; ii<2; ii++) {
#if 0 
	//V is a N x m+1  column-wise matrix
<<<<<<< mtl_classical_gram_schmidt.h
	itl::trans_mult(KS.V.sub_matrix(0, size(KS[i+1]), 0, i+1), 
			KS[i+1], KS.ss);
=======
	//this will only work for non complex element type
	//otherwise it should be M^H * V[i+1] instead of M^T * V[i+1]
	//(where M=V.V.sub_matrix(0, size(V[i+1])) 
	itl::trans_mult(V.V.sub_matrix(0, size(V[i+1]), 0, i+1), V[i+1], V.ss);
>>>>>>> 1.10
	
	//aggregate ss
	itl::all_reduce(KS.ss, Hi, i+1);
	
	itl::mult(KS.V.sub_matrix(0, size(KS[i+1]), 0, i+1), KS.ss, KS.w);
	itl::add(itl::scaled(KS.w, -1.0), KS[i+1]);
#else
	SP_NAME(dgemv)(&trans, &m, &i1, &alpha, KS.V.get_val(),
		       &m, KS[i+1].data(), &one, &beta, get_data(KS.ss), &one);
	
	//aggregate ss
	itl::all_reduce(KS.ss, Hi, i+1);
	
	//itl::mult(V.V.sub_matrix(0, size(V[i+1]), 0, i+1), V.ss, V.w);

	SP_NAME(dgemv)(&not_trans, &m, &i1, &alpha, KS.V.get_val(),
		       &m, get_data(KS.ss), &one, &beta, get_data(KS.w), &one);
#endif
      }

    }
      
    template <class Vec, class VecS, class VecX>
    void combine(classical_gram_schmidt<Vec>& KS,
		 const VecS& s, VecX& x, int i) {
      itl::mult(KS.V.sub_matrix(0, size(x), 0, i), s, x);
    }
  
}

#endif