strgy_cvmatrix.h

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: cvmatrix.h,v 1.1 2005/07/11 10:03:57 jppavone Exp $
 **********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_STRGY_CVMATRIX_HPP
#define realroot_SOLVE_SBDSLV_STRGY_CVMATRIX_HPP
//--------------------------------------------------------------------
#include <realroot/strgy_empty.h>
#include <realroot/linear_householder.hpp>
#include <realroot/linear_gramschmidt.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace realroot
{
  namespace strgy
  {
    template< class system, class _previous_ = strgy::empty<system> >
    struct cvmatrix : _previous_
    {
      typedef typename system::creal_t creal_t;
      typedef typename system::sz_t    sz_t;
      typedef typename system::interval_t interval_t;
      
      void symeigen2( creal_t * base, creal_t * m, bool normalise = false )
      {
    numerics::rounding< creal_t > rnd( numerics::rnd_nr() );
    creal_t b, c, delta;
    creal_t vp[2];
    // valeurs propres
    b = -(m[0]+m[3]);
    c =   m[0]*m[3] - m[1]*m[1];
    delta = b*b-4*c;
    assert(delta>=0);
    delta = std::sqrt(delta);
    vp[0] = -b - delta;
    vp[1] = -b + delta;
    // vecteurs propres
    base[0] = -m[1];
    base[1] =  m[0]-vp[0];
    creal_t sc; 
    if ( normalise ) sc = std::sqrt( base[0] * base[0] + base[1] * base[1] ) * vp[0];
    else sc = std::max(base[0],base[1]);
    base[0] /= sc;
    base[1] /= sc;
    base[2] =-(m[3]-vp[1]);
    base[3] =  m[1];
    if ( normalise ) sc = std::sqrt( base[2] * base[2] + base[3] * base[3] ) * vp[1];
    else sc = std::max(base[2],base[3]);
    base[2] /= sc;
    base[3] /= sc;
      };
      
      void symeigen( system * slv, creal_t * base, creal_t * m, unsigned n, const creal_t& eps = 1e-12 )
      {
    if ( n == 2 ) { symeigen2( base, m ); return; };
    creal_t vp[n];
    std::copy(m,m+n*n,base);
    linear::symeigen(base,vp,n);
    /* normalisation (vecteur propre / valeur propre) */
    unsigned i,j;
    for (  i = 0; i < n; i ++ )
      for ( j = i+1; j < n; j ++ )
        std::swap(base[i*n+j],base[j*n+i]);

    for (  i = 0; i < n ; i ++ )
      {
        //      creal_t s = 0.0;
        if ( std::abs(vp[i]) < 1e-3 ) 
          for ( unsigned k = 0; k < n; base[i*n+k] /= vp[i], k ++ ) ;
      };
      };
      
      void applycvmatrix( system * slv, sz_t * eqs, sz_t n  )
      {
    slv->restrict( eqs, n ); /* restriction des équations */
    slv->elevate ( eqs, n ); /* élévation des équations   */
    creal_t mean[n];
    creal_t cvmatrix[ n * n ];

    unsigned nsmp = slv->m_elenv.sz();

    for ( sz_t ie = 0; ie < n; ie ++ )
    {
      creal_t * data = slv->elldata(eqs[ie]);
      mean[ie] = 0;
      for ( unsigned k = 0; k < nsmp; k ++ )
        mean[ie] += (data[k]-mean[ie])/(k+1);
    };
    
    /* calcul de la partie supérieure de cvmatrix */
    for ( sz_t ie = 0; ie < n; ie ++ )
      {
        creal_t * data0 = slv->elldata(eqs[ie]); 
        creal_t      m0 = mean[ie];
        for ( sz_t id = ie; id < n; id ++ ) 
          {
        creal_t * data1 = slv->elldata(eqs[id]);
        creal_t m1 = mean[id];
        creal_t cm = 0.0;
        for ( int k = 0; k < slv->m_elenv.sz(); k ++ )
          cm += ((data0[k]-m0)*(data1[k]-m1)-cm)/(k+1);
        cvmatrix[ ie*n + id ] = cm;
          };
      };
    
    /* matrice ring_sparseetrique */
    for ( sz_t ie = 0; ie < n; ie ++ )
      for ( sz_t id = 0; id < ie; id ++ )
        cvmatrix[ie*n+id] = cvmatrix[id*n+ie];
    /* calcul de la base des vecteurs propres */
    creal_t base[n*n];
    symeigen( slv, base, cvmatrix, n ); 
    slv->transform( base, eqs, n ); /* transformation du systeme */
    slv->project();
      };

      bool process( system * slv )
      {
    if ( ! _previous_::process(slv) ) return false;
    applycvmatrix( slv, slv->m_eqssel, slv->neqs() );
    return true;
      };
    };
  };
};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //