La QR décomposition


Installer et compiler ces fichiers dans votre répertoire de travail.


c00a.c
/* ------------------------------------ */
/*  Save as :   c00a.c                  */
/* ------------------------------------ */
#include "w_a.h"
/* ------------------------------------ */
void fun(int r)
{
double **A    = i_mZ(r,r);
double **b    = i_mZ(r,C1);
double **Ab   = i_Abr_Ac_bc_mZ(r,r,C1);

double **Q    = i_mZ(r,r);
double **R    = i_mZ(r,r);
double **Q_T  = i_mZ(r,r);
double **invR = i_mZ(r,r);

double **invR_Q_T = i_mZ(r,r);
double **x        = i_mZ(r,C1); // x invR * Q_T * b

   do{
   r_mZ(A,99.);
   printf(".");
  }while(!det_Z(A).r||!det_Z(A).i);

  clrscrn();
  printf(" A :");
  p_mZ(A,S5,P0,S3,P0,C6);
  printf(" b :");
  r_mZ(b,99.);
  p_mZ(b,S5,P0,S3,P0,C6);
  printf(" Ab :");
  c_A_b_Ab_mZ(A,b,Ab);
  p_mZ(Ab,S5,P0,S3,P0,C6);
  stop();
    
  clrscrn();
  QR_mZ(A,Q,R);    
  printf(" Q :");
  p_mZ(Q,S12,P5,S10,P5,C6); 
  printf(" R :");
  p_mZ(R,S12,P5,S10,P5,C6);
  stop();

  clrscrn();
  ctranspose_mZ(Q,Q_T);   
  printf(" Q_T :");
  pE_mZ(Q_T,S9,P5,S7,P5,C6);
  inv_mZ(R,invR); 
  printf(" invR :");
  pE_mZ(invR,S9,P5,S7,P5,C6);
  stop();
  
  clrscrn();
  printf(" Copy/Past into the octave window.\n\n");
  p_Octave_mZ(Ab,"Ab",P0,P0);
  printf(" rref(Ab,.00000000001)\n\n");

  printf(" gj_mZ(Ab) :");
  gj_mZ(Ab);
  p_mZ(Ab,S9,P5,S7,P5,C6);

  mul_mZ(invR,Q_T,invR_Q_T);
  mul_mZ(invR_Q_T,b,x);
  printf(" x = invR * Q_T * b :");
  p_mZ(x,S9,P5,S7,P5,C6);

  f_mZ(A);
  f_mZ(b);
  f_mZ(Ab);
  f_mZ(Q);
  f_mZ(Q_T);
  f_mZ(R);
  f_mZ(invR);  
  f_mZ(invR_Q_T); 
  f_mZ(x); 
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));
do
{
  fun(R3);


} while(stop_w());

  return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */


On obtient le même résultat qu'avec la méthode de gauss jordan


Exemple de sortie écran :
 A :
  +76.000 +13.00i   -47.000 -20.00i   -86.000  +6.00i 
  -88.000 +82.00i   -10.000  -2.00i   -95.000 -31.00i 
  +53.000 -12.00i   +34.000 +36.00i   +23.000 +83.00i 

 b :
 -578.000 +58.00i 
  -72.000+141.00i 
 +663.000+890.00i 

 Ab :
  +76.000 +13.00i   -47.000 -20.00i   -86.000  +6.00i  -578.000 +58.00i 
  -88.000 +82.00i   -10.000  -2.00i   -95.000 -31.00i   -72.000+141.00i 
  +53.000 -12.00i   +34.000 +36.00i   +23.000 +83.00i  +663.000+890.00i 

 Press return to continue.  


 Q :
   +0.497  +0.09i    -0.578  -0.39i    -0.510  +0.04i 
   -0.576  +0.54i    -0.118  +0.19i    -0.460  +0.35i 
   +0.347  -0.08i    +0.531  +0.43i    -0.615  -0.17i 

 R :
 +152.859  +0.00i   -11.422 +15.72i    -2.728+109.69i 
   +0.000  +0.00i   +69.191  +0.00i  +100.425 +19.11i 
   +0.000  -0.00i    -0.000  +0.00i   +48.443  +0.00i 

 Press return to continue. 


 Q_T :
+4.972e-01-8.50e-02i -5.757e-01-5.36e-01i +3.467e-01+7.85e-02i 
-5.779e-01+3.88e-01i -1.177e-01-1.90e-01i +5.308e-01-4.29e-01i 
-5.098e-01-3.51e-02i -4.598e-01-3.45e-01i -6.148e-01+1.74e-01i 

 invR :
+6.542e-03-4.49e-19i +1.080e-03-1.49e-03i -2.457e-03-1.22e-02i 
+0.000e+00+0.00e+00i +1.445e-02-4.15e-18i -2.996e-02-5.70e-03i 
+0.000e+00-0.00e+00i -0.000e+00+0.00e+00i +2.064e-02-3.24e-18i 

 Press return to continue. 


 Copy/Past into the octave window.

 Ab=[
+76+13*i,-47-20*i,-86+6*i,-578+58*i;
-88+82*i,-10-2*i,-95-31*i,-72+141*i;
+53-12*i,+34+36*i,+23+83*i,+663+890*i]

 rref(Ab,.00000000001)

 Press return to continue. 


 gj_mZ(Ab) :
   +1.000  +0.00i    -0.000  -0.00i    +0.000  +0.00i    -4.373  +4.33i 
   +0.000  +0.00i    +1.000  -0.00i    +0.000  +0.00i   +18.379 +14.42i 
   +0.000  +0.00i    +0.000  +0.00i    +1.000  -0.00i    -3.797  -9.93i 

 x = invR * Q_T * b :
   -4.373  +4.33i 
  +18.379 +14.42i 
   -3.797  -9.93i 


 Press   return to continue
 Press X return to stop