gram_schmidt.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */


#ifndef SLIP_GRAM_SCHMIDT_HPP
#define SLIP_GRAM_SCHMIDT_HPP
#include <cassert>
#include <cmath>
#include <numeric>
#include "DPoint2d.hpp"
#include "linear_algebra.hpp"
#include "linear_algebra_traits.hpp"

namespace slip
{

  template<class RandomAccessIterator, class InnerProduct>
  void gram_schmidt_normalization(RandomAccessIterator init_base_first,
                                                                  RandomAccessIterator init_base_end,
                                                                  RandomAccessIterator ortho_base_first,
                                                                  RandomAccessIterator ortho_base_end,
                                                                  InnerProduct inner_prod)
{
  
 
  assert( (init_base_end - init_base_first) == (ortho_base_end - ortho_base_first));
 
  typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
  _Distance init_base_first_size = (init_base_end - init_base_first);
  ortho_base_first[0] = init_base_first[0] / std::sqrt(inner_prod(init_base_first[0],init_base_first[0]));
  for(_Distance k = 1; k < init_base_first_size; ++k)
    {
      typename RandomAccessIterator::value_type sum = (inner_prod(ortho_base_first[0],init_base_first[k]) /  inner_prod(ortho_base_first[0],ortho_base_first[0])) * ortho_base_first[0];
      for(_Distance q = 1; q < k; ++q)
                {
                  sum += (inner_prod(ortho_base_first[q],init_base_first[k]) /  inner_prod(ortho_base_first[q],ortho_base_first[q])) * ortho_base_first[q];
                }
      ortho_base_first[k] = init_base_first[k] - sum;
      ortho_base_first[k] = ortho_base_first[k] /  std::sqrt(inner_prod(ortho_base_first[k],ortho_base_first[k]));
    }
}


  template<class RandomAccessIterator, class InnerProduct>
  void modified_gram_schmidt_normalization(RandomAccessIterator init_base_first,
                                                                                   RandomAccessIterator init_base_end,
                                                                                   RandomAccessIterator ortho_base_first,
                                                                                   RandomAccessIterator ortho_base_end,
                                                                                   InnerProduct inner_prod)
  {
  
   assert( (init_base_end - init_base_first) == (ortho_base_end - ortho_base_first));
   
   typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
   _Distance init_base_first_size = (init_base_end - init_base_first);
   
   std::copy(init_base_first,init_base_end,ortho_base_first);

  for(_Distance k = 0; k < init_base_first_size; ++k)
    {
      ortho_base_first[k] /=  std::sqrt(inner_prod(ortho_base_first[k],ortho_base_first[k]));
      for(_Distance q = (k + 1); q < init_base_first_size; ++q)
                {
                  ortho_base_first[q] -= inner_prod(ortho_base_first[q],ortho_base_first[k]) * ortho_base_first[k];
                } 
    }
}


  template<class RandomAccessIterator, class InnerProduct>
  void gram_schmidt_orthogonalization(RandomAccessIterator init_base_first,
                                                                      RandomAccessIterator init_base_end,
                                                                      RandomAccessIterator ortho_base_first,
                                                                      RandomAccessIterator ortho_base_end,
                                                                      InnerProduct inner_prod)
{
  assert( (init_base_end - init_base_first) == (ortho_base_end - ortho_base_first));
  typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
  _Distance init_base_first_size = (init_base_end - init_base_first);
  ortho_base_first[0] = init_base_first[0];
  for(_Distance k = 1; k < init_base_first_size; ++k)
    {
      typename RandomAccessIterator::value_type sum = (inner_prod(ortho_base_first[0],init_base_first[k]) /  inner_prod(ortho_base_first[0],ortho_base_first[0])) * ortho_base_first[0];
      for(_Distance q = 1; q < k; ++q)
                {
                  sum += (inner_prod(ortho_base_first[q],init_base_first[k]) /  inner_prod(ortho_base_first[q],ortho_base_first[q])) * ortho_base_first[q];
                }
      ortho_base_first[k] = init_base_first[k] - sum;
    }
}


  template<class RandomAccessIterator, class InnerProduct>
  void modified_gram_schmidt_orthogonalization(RandomAccessIterator init_base_first,
                                                                                       RandomAccessIterator init_base_end,
                                                                                       RandomAccessIterator ortho_base_first,
                                                                                       RandomAccessIterator ortho_base_end,
                                                                                       InnerProduct inner_prod)
  {
  
  assert( (init_base_end - init_base_first) == (ortho_base_end - ortho_base_first));
  typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
  _Distance init_base_first_size = (init_base_end - init_base_first);

  std::copy(init_base_first,init_base_end,ortho_base_first);

  for(_Distance k = 0; k < init_base_first_size; ++k)
    {
       typedef typename std::iterator_traits<RandomAccessIterator>::value_type BaseElement;

       BaseElement Bk = ortho_base_first[k] /  std::sqrt(inner_prod(ortho_base_first[k],ortho_base_first[k]));
      for(_Distance q = (k + 1); q < init_base_first_size; ++q)
                {
                  ortho_base_first[q] -= inner_prod(ortho_base_first[q],Bk) * Bk;
                } 
    }
}


 

  template<typename MatrixIterator1,
                   typename MatrixIterator2,
                   typename MatrixIterator3>
  void modified_gram_schmidt(MatrixIterator1 A_up,
                                                     MatrixIterator1 A_bot,
                                                     MatrixIterator2 Q_up,
                                                     MatrixIterator2 Q_bot,
                                                     MatrixIterator3 R_up,
                                                     MatrixIterator3 R_bot)
{
  assert((A_bot - A_up)[0] == (Q_bot - Q_up)[0]);
  assert((A_bot - A_up)[1] == (Q_bot - Q_up)[1]);
  assert((R_bot - R_up)[0] == (A_bot - A_up)[1]);
  assert((R_bot - R_up)[1] == (A_bot - A_up)[1]);
  assert((A_bot - A_up)[0] >= (A_bot - A_up)[1]);
  
    typedef typename std::iterator_traits<MatrixIterator1>::value_type value_type;
    typedef typename slip::lin_alg_traits<value_type>::value_type norm_type;
    typedef  typename MatrixIterator1::size_type size_type;
    typename MatrixIterator1::difference_type sizeA = A_bot - A_up;
    typename MatrixIterator2::difference_type sizeQ = Q_bot - Q_up;
    typename MatrixIterator3::difference_type sizeR = R_bot - R_up;

    size_type A_cols = sizeA[1];
    

   //init Q with 0
    std::fill(Q_up,Q_bot,value_type());
  //copy the first column of A in the first column of Q
  std::copy(A_up.col_begin(0),A_up.col_end(0),Q_up.col_begin(0));
  //init R with 0 and R[0][0] with 1
  std::fill(R_up,R_bot,value_type());
  *R_up = value_type(1);
  for(size_type k = 0; k < A_cols; ++k)
    {
      slip::DPoint2d<int> d(k,k);
      norm_type norm_col_k = std::sqrt(slip::L22_norm_cplx(A_up.col_begin(k),
                                                                                                                   A_up.col_end(k)));
      R_up[d] = norm_col_k;
      // divides the column k elements of A by norm_col_k and copy the result 
      //in the column k elements of Q
      std::transform(A_up.col_begin(k),A_up.col_end(k),
                                     Q_up.col_begin(k),std::bind2nd(std::divides<value_type>(),norm_col_k));
      for(size_type q = (k+1);  q < A_cols; ++q)
                {
                  slip::DPoint2d<int> d2(k,q);
                  R_up[d2] = slip::hermitian_inner_product(Q_up.col_begin(k),Q_up.col_end(k),
                                                                                                  A_up.col_begin(q),value_type());
                  slip::reduction(A_up.col_begin(q),A_up.col_end(q),
                                                  Q_up.col_begin(k),(-R_up[d2]));
                }
                                     
    } 
 }

  template<class RandomAccessIterator, class RandomAccessIterator2d, class InnerProduct>
  void gram_matrix(RandomAccessIterator init_base_first,
                                   RandomAccessIterator init_base_end,
                                   RandomAccessIterator2d matrix_upper_left,
                                   RandomAccessIterator2d matrix_bottom_right,
                                   InnerProduct inner_prod)
{
  assert((matrix_bottom_right - matrix_upper_left)[0] == (init_base_end - init_base_first));
  assert((matrix_bottom_right - matrix_upper_left)[1] == (init_base_end - init_base_first));
 
  typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
  typedef typename std::iterator_traits<RandomAccessIterator2d>::difference_type _Distance2d;

  _Distance init_base_first_size = (init_base_end - init_base_first);  
  _Distance2d matrix_size2d = (matrix_bottom_right - matrix_upper_left);
  
  for(_Distance i = 0 ; i < init_base_first_size; ++i)
    {
      for(_Distance j = 0 ; j < init_base_first_size; ++j)
                {
                  *(matrix_upper_left++) = inner_prod(init_base_first[i],init_base_first[j]);
                }
    }

}


}//slip::


#endif //SLIP_GRAM_SCHMIDT_HPP