linear_gramschmidt.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: gramschmidt.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_LINEAR_GRAMSCHMIDT_HPP
#define realroot_SOLVE_SBDSLV_LINEAR_GRAMSCHMIDT_HPP
//--------------------------------------------------------------------
#include <assert.h>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace linear
{

  template<typename real_t>
  bool gramschmidt( real_t * M, unsigned m, unsigned n )
  {
    unsigned i,j,k;
    assert(m<=n);
    for ( i = 0; i < m; i ++ )
      {
    real_t s = 0.0;
    for ( k = 0; k < n; k ++ )
      s += M[i*n+k]*M[i*n+k];
    if ( s < numerics::epsilon<real_t>::result ) return false;
    for ( j = i+1; j < m; j ++ )
      {
        real_t d = 0.0;
        for ( k = 0; k < n; k ++ )
          d += M[i*n+k]*M[j*n+k];
        d /= s;
        for ( k = 0; k < n; k ++ )
          M[j*n+k] -= d*M[i*n+k];
      };
      };
    return true;
  };
  
  /* on suppose que M admet un noyau de dimension 1 
   * ( et puis c tout !)
   */
  
  template<typename real_t> 
  void library( real_t * M, unsigned n, real_t * v  )
  {

    //int m = n;
    int goods[n];
    int c = 0;
    int i,j,k,p;
    //int already =0;
    real_t sc[ n ];
    
    c = 0;
    for ( i = 0; i < n; i ++ )
      {
    real_t s = 0.0;
    for ( k = 0; k < n; s += M[i*n+k]*M[i*n+k], k ++ ) ;
    if ( s < numerics::epsilon<real_t>::result )   continue; 
    goods[c++] = i;
    for ( j = i+1; j < n; j ++ )
      {
        real_t d = 0.0;
        for ( k = 0; k < n; k ++ )
          d += M[i*n+k]*M[j*n+k];
        d /= s;
        for ( k = 0; k < n; k ++ )
          M[j*n+k] -= d*M[i*n+k];
      };
    sc[i] = s;
      };
    if ( c == n ) return;
    k = 0;
    while( k < c && goods[k] == k ) k ++;
    while ( k < c )
      {
    std::copy( M + goods[k]*n, M + (goods[k]+1)*n, 
           M + k );
    sc[k] = sc[goods[k]];
    k++;
      };
    //    vctops::print(goods,c);
    //    std::cout << std::endl;
    //    vctops::print(M,n);
    //    std::cout << std::endl;
    //    vctops::print(M+n,n);
    //    std::cout << std::endl;
    //    std::cout << "c = " << c << std::endl;
    //    std::cout << "n = " << n << std::endl;
    for ( i = c; i < n; i ++ )
      for ( k = 0; k < n; k ++ )
    M[i*n+k] = ((real_t)rand())/((real_t)RAND_MAX); 
    //    vctops::print(M+c*n,n);
    //    std::cout << std::endl;
    for ( i = c; i  < n; i ++ )
      {
    for ( p = i-1; p >= 0; p-- )
      {
        //      std::cout << sc[p] << std::endl;
        real_t d = 0.0;
        for ( k = 0; k < n; k ++ )
          d += M[p*n+k]*M[i*n+k];
        d /= sc[p];
        for ( k = 0; k < n; k ++ )
          {
        //      std::cout << M[i*n+k] << std::endl;
        M[i*n+k] -= d*M[p*n+k];
        //      std::cout << M[i*n+k] << std::endl;
          };
      };
    real_t s = 0.0;
    for ( k = 0; k < n; k ++ )
      s += M[i*n+k]*M[i*n+k];
    sc[i] = s;
      };
  };

  
};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //