CoreFlows/Models/inc/utilitaire_algebre.h

   1 #ifndef utilitaire_algebre_h
   2 #define utilitaire_algebre_h
   3
   4 /*      rpoly.cpp -- Jenkins-Traub real polynomial root finder.
   5  *
   6  *      (C) 2000, C. Bond.  All rights reserved.
   7  *
   8  *      Translation of TOMS493 from FORTRAN to C. This
   9  *      implementation of Jenkins-Traub partially adapts
  10  *      the original code to a C environment by restruction
  11  *      many of the 'goto' controls to better fit a block
  12  *      structured form. It also eliminates the global memory
  13  *      allocation in favor of local, dynamic memory management.
  14  *
  15  *      The calling conventions are slightly modified to return
  16  *      the number of roots found as the function value.
  17  *
  18  *      INPUT:
  19  *      op - double precision vector of coefficients in order of
  20  *              decreasing powers.
  21  *      degree - integer degree of polynomial
  22  *
  23  *      OUTPUT:
  24  *      zeror,zeroi - output double precision vectors of the
  25  *              real and imaginary parts of the zeros.
  26  *
  27  *      RETURN:
  28  *      returnval:   -1 if leading coefficient is zero, otherwise
  29  *                  number of roots found.
  30  */
  31
  32 /* A roots_polynoms class */
  33
  34 /*! \class roots_polynoms utilitaire_algebre.hxx "utilitaire_algebre.hxx"
  35  *  \brief  Computes the roots of a given polynomial
  36  *  \details  Translation of TOMS493 from FORTRAN to C. This
  37  *      implementation of Jenkins-Traub partially adapts
  38  *      the original code to a C environment by restruction
  39  *      many of the 'goto' controls to better fit a block
  40  *      structured form. It also eliminates the global memory
  41  *      allocation in favor of local, dynamic memory management.
  42  *      The calling conventions are slightly modified to return
  43  *      the number of roots found as the function value.
  44  */
  45 #include <math.h>
  46 #include <vector>
  47 #include <complex>
  48 #include <iostream>
  49
  50 using namespace std;
  51
  52 class roots_polynoms
  53 {
  54 public:
  55
  56         roots_polynoms();
  57
  58         /** \fn quad
  59          * \brief calcule les racines d'un polynome de degré 2
  60          * \Details Calculate the zeros of the quadratic a*x² + b*x + c.
  61          *  The quadratic formula, modified to avoid overflow, is used
  62          *  to find the larger zero if the zeros are real and both
  63          *  are complex. The smaller real zero is found directly from
  64          *  the product of the zeros c/a
  65          * @param a double correspond au coefficient de x²
  66          * @param b double correspond au coefficient de x
  67          * @param c double correspond au coefficient de 1
  68          * @param *sr double correspond à la pratie réelle de la premiere racine
  69          * @param *si double correspond à la pratie imaginaire de la premiere racine
  70          * @param *lr double correspond à la pratie réelle de la deuxième racine
  71          * @param *li double correspond à la pratie imaginaire de la deuxième racine
  72          * @return calcule les racine du polynome dans sr+i*si et lr+i*li
  73          * */
  74         void quad(double a,double b,double c,double *sr,double *si, double *lr,double *li);
  75
  76
  77         /** \fn fxshfr
  78          * \Details  Computes up to L2 fixed shift k-polynomials,
  79          *  testing for convergence in the linear or quadratic
  80          *  case. Initiates one of the variable shift
  81          *  iterations and returns with the number of zeros found
  82          * @param l2 entier
  83          * @param *nz entier
  84          * @return */
  85         void fxshfr(int l2, int *nz);
  86
  87         /** \fn quadit
  88          * \brief
  89          * \Details /*  Variable-shift k-polynomial iteration for a
  90          *  quadratic factor converges only if the zeros are
  91          *  equimodular or nearly so.
  92          * @param *uu, *vv double coefficients of starting quadratic
  93          * @param *nz entier number of zeros found
  94          * @return */
  95         void quadit(double *uu,double *vv,int *nz);
  96
  97         /** \fn realit
  98          * \brief Variable-shift H polynomial iteration for a real zero
  99          * @param  *sss starting iterate (double)
 100          * @param *nz number of zeros found (int)
 101          * @param *iflag flag to indicate a pair of zeros near real axis (int)
 102          * @return */
 103         void realit(double *sss, int *nz, int *iflag);
 104
 105         /** \fn calcsc
 106          * \Details This routine calculates scalar quantities used to
 107          *  compute the next k polynomial and new estimates of
 108          *  the quadratic coefficients.
 109          * @param *type integer variable set here
 110          * indicating how the calculations are normalized to avoid overflow
 111          * @return */
 112         void calcsc(int *type);
 113
 114         /** \fn nextk
 115          * \brief Computes the next k polynomials using scalars computed in calcsc().
 116          * @param *type integer variable set here
 117          * indicating how the calculations are normalized to avoid overflow
 118          * @return */
 119         void nextk(int *type);
 120
 121         /** \fn newest
 122          * \brief Compute new estimates of the quadratic coefficients
 123          *  using the scalars computed in calcsc ().
 124          * @param type
 125          * @param *uu, *vv double coefficients of starting quadratic
 126          * @return */
 127         void newest(int type,double *uu,double *vv);
 128
 129         /** \fn quadsd
 130          * \brief Divides p by the quadratic 1,u,v placing the quotient
 131          *  in q and the remainder in a,b.
 132          * @return */
 133         void quadsd(int n,double *u,double *v,double *p,double *q,
 134                         double *a,double *b);
 135
 136         /** \fn rpoly
 137          * \brief
 138          * @param *op
 139          * @param degree
 140          * @param *zeror
 141          * @param *zeroi
 142          * @return */
 143         int rpoly(double *op, int degree, double *zeror, double *zeroi);
 144
 145 private:
 146         double infin;
 147         double smalno;
 148         double eta;
 149         double base;
 150         double *p,*qp,*k,*qk,*svk;
 151         double sr,si,u,v,a,b,c,d,a1,a2;
 152         double a3,a6,a7,e,f,g,h,szr,szi,lzr,lzi;
 153         double are,mre;
 154         int n,nmi;
 155 };
 156
 157 class Polynoms
 158 {
 159 public:
 160
 161         void abs_par_interp_directe(int n,  vector< complex< double > > vp,   double * A, int sizeA, double epsilon, double *result,vector< complex< double > > y);
 162
 163         /** \fn polynome_caracteristique
 164          * \brief Calcule le polynome caracteristique
 165          * @return */
 166         vector< double > polynome_caracteristique(double alpha_v, double alpha_l, double Uv, double Ul, double rhov, double rhol, double   cv_2, double cl_2, double dPiv, double dPil);//c= inverse vitesse du son!
 167         /** \fn polynome_caracteristique
 168          * \brief
 169          * @param alpha1
 170          * @param alpha2
 171          * @param u1
 172          * @param u2
 173          * @param rho1
 174          * @param rho2
 175          * @param invc1sq
 176          * @param invc2sq
 177          * @param dpi1
 178          * @param dpi2
 179          * @param g2press
 180          * @param g2alpha
 181          * @param g2
 182          * @param epsilon
 183          * @return */
 184         vector< double > polynome_caracteristique(double alpha1, double alpha2, double u1, double u2, double rho1, double rho2, double invc1sq, double invc2sq, double dpi1, double dpi2, double g2press, double g2alpha, double g2, double epsilon);
 185
 186         /** \fn module
 187          * \brief calcule le module d'un nombre complexe
 188          * @param z est un nombre complexe
 189          * @return  calcule le module de z*/
 190         double module(complex< double > z);
 191         /** \fn modulec calcule le module² de z */
 192         double modulec(complex< double > z);
 193         /** \fn abs_generalise
 194          * \brief calcule la valeur absolue d'un nombre complexe
 195          * \Details calcule la valeur absolue d'un nombre complexe en prenant en compte
 196          * que la partie réelle c-a-d si z= a+ib abs_generalize() renvoie |a|+ib
 197          * @param z
 198          * @return si z = a+ib elle renvoie |a|+ib */
 199         complex< double > abs_generalise(complex< double > z);
 200
 201         int new_tri_selectif(vector< complex< double > > &L, int n, double epsilon);
 202
 203         /** \fn tri_selectif
 204          * \brief
 205          * @param L
 206          * @param n
 207          * @param epsilon
 208          * @return */
 209         template< typename T >
 210         int tri_selectif(T& L, int n, double epsilon);
 211
 212         /** \fn matrixProduct
 213          * \brief calcule le produit de 2 matrices .
 214          * @param *Matrix1 la 1ere matrice
 215          * @param R1,C1 la taille de la matrice1 Nbr de ligne et le nbr de colonne
 216          * @param Matrix2 la 2ieme matrice
 217          * @param R2,C2 la taille de la matrice2 Nbr de ligne et le nbr de colonne
 218          * @param out le resultat de Matrix1*Matrix2
 219          * @return */
 220         void matrixProduct
 221         (
 222                         const double *Matrix1,
 223                         const int &R1,
 224                         const int &C1,
 225                         const double *Matrix2,
 226                         const int &R2,
 227                         const int &C2,
 228                         double *out
 229         );
 230
 231         /** \fn matrixProdVec
 232          * \brief Calcule le produit Matrice vecteur
 233          * @param *Matrix correspond à la matrice du produit
 234          * @param R1,C1 la taille de la matrice( Nbr de ligne et le nbr de colonne)
 235          * @param Vector correspond au vecteur du produit
 236          * @param out le resultat de Matrix*Vector
 237          * @return   */
 238         void matrixProdVec
 239         (       const double *Matrix,
 240                         const int &R1,
 241                         const int &C1,
 242                         const double *Vector,
 243                         double *out
 244         );
 245
 246         bool belongTo(complex< double > a , vector< complex <double > > v, double epsilon);
 247
 248 private:
 249
 250 /** \fn add
 251  * \brief Calcule la somme de 2vecteurs ( U=U+V)
 252  * @param *U,*V les 2 vecteurs
 253  * @param N taille des 2 vecteurs
 254  * @return la somme de U et V dans U*/
 255 void add
 256 (
 257                 double *U,              //vecteur auquel on ajoute
 258                 int N,                  //taille du vecteur
 259                 const double *V //ajout
 260 );
 261 /** \fn tensor
 262  * \brief Calcule le tenseur de 2 vecteurs
 263  * @param *a le 1er vecteur (*double)
 264  * @param N la taille du premier vecteur (int)
 265  * @param *b le 2ieme vecteur (*double)
 266  * @param M la taille du 2ieme vecteur (int)
 267  * @param *ab le resultat
 268  * @return */
 269
 270 void tensor
 271 (       const double *a,        //vecteur gauche
 272                 int N,                  //taille du vecteur gauche
 273                 const double *b,        //vecteur droit
 274                 int M,                  //taille du vecteur droit
 275                 double *ab              //conteneur de sortie
 276 );
 277
 278 /** \fn shift_diagonal
 279  * \brief Calcule la trace d'une matrice carrée
 280  * @param *Matrix correspond  à la matrice
 281  * @param N taille de la matrice
 282  * @param shift résultat
 283  * @return renvoie la trace de la matrice Matrix dans shift */
 284 void shift_diagonal( double *Matrix, const int N, double shift);
 285 //qques outils d'algebre des polynomes
 286
 287 /** \fn remplir_trinome
 288  * \brief remplie un trinome
 289  */
 290 void remplir_trinome(double coeff, double u, double alpha, vector< double >& trinome);
 291
 292 /** \fn additionne
 293  * \brief Additionne 2 vecteurs de tailles differentes
 294  * @param P,Q les 2 vecteurs à additionner
 295  * @param n,m les 2 tailles respectives des vecteurs P et Q
 296  * @return un vecteur qui correspond a P+Q*/
 297 vector< double > additionne(const vector< double >& P, const vector< double >& Q, int n, int m);
 298
 299 /** \fn multiplie
 300  * \brief Calcule le produit scalaire de 2 vecteurs de tailles differentes
 301  * @param P,Q les 2 vecteurs
 302  * @param n,m les 2 tailles respectives des 2 vecteurs
 303  * @return un vecteur correspond à (P,Q) */
 304 vector< double > multiplie(const vector< double >& P, const vector< double >& Q, int n, int m);
 305
 306 //Pour le tri des valeurs propres
 307
 308 void ordre_croissant_abs(vector< complex< double > > &L, int n);
 309
 310 //Calcul des coefficients du polynome d'interpolation x->y dans la base de Newton
 311 template<class T>
 312
 313 /** \fn dif_div
 314  * \brief
 315  * @param n
 316  * @param x
 317  * @param y
 318  * @param epsilon
 319  * @return */
 320 vector< complex< double > > dif_div(int n, const vector< complex< double > >& x, const vector< T >& y, double epsilon);//n=nb valeurs à interpoler
 321
 322 //attention, Les vp complexes conjugu�es doivent se suivre dans la liste x
 323 void appliquer_dif_div(int n, const vector< complex< double > >& dif_div, const vector< complex< double > >& x, const double* A, const int sizeA, double epsilon, double *p);//p=result
 324
 325 double avr(double a, double b);
 326
 327 };
 328 #endif