Fonctions matricielles dans C. Matrices Non symétriques conjugués


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


c00a.c
/* ------------------------------------ */
/*  Save as :   c00a.c                  */
/* ------------------------------------ */
#include "w_a.h"
/* ------------------------------------ */
#define FACTOR_E    +1.E-0         
#define RCA          RC3  
/* ------------------------------------ */
/* ------------------------------------ */
double **EV_cos2A_mZ(
double **EV,
double **EV_cos2A
)
{
int  r;
int  c;

nb_Z cos2A;

  for (   r=R1; r<EV[R_SIZE][C0]; r++)
    for ( c=C1; c<EV[C_SIZE][C0]; c+=C2)

        if(((r-R1)*C2)==c-C1)
          {
           cos2A = cos_Z(i_Z(2.*EV[r][c],2.*EV[r][c+C1]));
           
           EV_cos2A[r][c   ] = cos2A.r;
           EV_cos2A[r][c+C1] = cos2A.i;
          }
              
 return(EV_cos2A);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **EV_cosAP2_mns_sinAP2_mZ(
double **EV,
double **EV_cosAP2_mns_sinAP2
)
{
int  r;
int  c;

nb_Z cosAP2_mns_sinAP2;

nb_Z cos;
nb_Z cosP2;

nb_Z sin;
nb_Z sinP2;

  for (   r=R1; r<EV[R_SIZE][C0]; r++)
    for ( c=C1; c<EV[C_SIZE][C0]; c+=C2)

        if(((r-R1)*C2)==c-C1)
          {
           cos   = cos_Z(i_Z(EV[r][c],EV[r][c+C1]));
           cosP2 = mul_Z(cos,cos);			  
			  
           sin   = sin_Z(i_Z(EV[r][c],EV[r][c+C1]));
           sinP2 = mul_Z(sin,sin);
           
           cosAP2_mns_sinAP2 = sub_Z(cosP2,sinP2);
           
           EV_cosAP2_mns_sinAP2[r][c   ] = cosAP2_mns_sinAP2.r;
           EV_cosAP2_mns_sinAP2[r][c+C1] = cosAP2_mns_sinAP2.i;
          }
              
 return(EV_cosAP2_mns_sinAP2);
}
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{                               
double **A                            = i_mZ(RCA,RCA);

double **V                            = i_mZ(RCA,RCA);
double **invV                         = i_mZ(RCA,RCA);
double **T                            = i_mZ(RCA,RCA);

double **EV                    = i_mZ(RCA,RCA);
double **EV_cos2a              = i_mZ(RCA,RCA);
double **EV_cosAP2_mns_sinAP2  = i_mZ(RCA,RCA);

double **cos2a                        = i_mZ(RCA,RCA);
double **cosAP2_mns_sinAP2            = i_mZ(RCA,RCA);

   do
  {
       rE_mZ(A,999,+1.E-4);  
   eigs_V_mZ(A,V,FACTOR_E);
  }while(!det_Z(V).r);
  
  clrscrn();
  printf(" A :");
  p_mZ(A, S9,P4, S8,P4, C3);

  printf(" V :");
  eigs_V_mZ(A,V,FACTOR_E);
  pE_mZ(V, S12,P4, S12,P4, C3);
  
  printf(" inv(V) ... Some time the matrix is not invertible :");
  inv_mZ(V,invV);
  pE_mZ(invV, S12,P4, S12,P4, C3);
  
  printf(" EV = invV * A * V");
  mul_mZ(invV,A,T);
  mul_mZ(T,V,EV); 
  pE_mZ(clean_eyes_mZ(EV), S12,P4, S12,P4, C3);
  stop();

  clrscrn();  
  printf(" cos(2*A)        :");
  EV_cos2A_mZ(EV,EV_cos2a);
   
  mul_mZ(V,EV_cos2a,T);
  mul_mZ(T,invV,cos2a); 
  pE_mZ(cos2a, S12,P4, S8,P4, C3);

  printf(" cos(A)**2-sin(A)**2 :");
  EV_cosAP2_mns_sinAP2_mZ(EV,EV_cosAP2_mns_sinAP2);
   
  mul_mZ(V,EV_cosAP2_mns_sinAP2,T);
  mul_mZ(T,invV,cosAP2_mns_sinAP2); 
  pE_mZ(cosAP2_mns_sinAP2, S12,P4, S8,P4, C3);  
  
  f_mZ(A);
  
  f_mZ(V);  
  f_mZ(invV);  
  f_mZ(T);  
  
  f_mZ(EV);
  
  f_mZ(EV_cos2a);   
  f_mZ(cos2a); 
  
  f_mZ(EV_cosAP2_mns_sinAP2);   
  f_mZ(cosAP2_mns_sinAP2);     
}
/* ------------------------------------ */
int main(void)
{
time_t t;

  srand(time(&t));

do
{
    fun();
    
} while(stop_w());

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


Exemple de sortie écran :
 A :
  +0.0079 -0.0683i   +0.0472 -0.0940i   +0.0943 +0.0500i 
  -0.0662 -0.0341i   +0.0673 +0.0078i   -0.0565 -0.0033i 
  -0.0678 -0.0344i   +0.0182 -0.0402i   +0.0885 +0.0365i 

 V :
 +6.0580e-01 +4.5086e-01i  +3.1697e-01 -5.2921e-01i  -2.8260e-01 +1.8510e-01i 
 +4.3666e-01 -2.2596e-01i  +2.6245e-01 +6.3163e-01i  -5.1110e-02 -8.0316e-01i 
 +4.3359e-01 +6.0173e-18i  +3.8940e-01 +7.4493e-19i  +4.8806e-01 -1.1231e-18i 

 inv(V) ... Some time the matrix is not invertible :
 +6.7605e-01 -6.2461e-01i  +6.0221e-01 +9.9447e-02i  +5.3979e-02 +3.8337e-01i 
 -1.3294e-01 +4.6378e-01i  -1.1964e-01 -6.2436e-01i  +1.1139e+00 +5.6689e-02i 
 -4.9454e-01 +1.8487e-01i  -4.3955e-01 +4.0981e-01i  +1.1123e+00 -3.8582e-01i 

 EigsValue = invV * A * V
 +2.6923e-02 -1.3203e-01i  +0.0000e+00 +0.0000e+00i  +0.0000e+00 +0.0000e+00i 
 +0.0000e+00 +0.0000e+00i  +6.4033e-02 +1.0307e-01i  +0.0000e+00 +0.0000e+00i 
 +0.0000e+00 +0.0000e+00i  +0.0000e+00 +0.0000e+00i  +7.2744e-02 +4.9646e-03i 

 Press return to continue. 


 cos(2*A)        :
 +1.0313e+00+6.2536e-03i  -3.1624e-03+2.5763e-02i  -1.5574e-02-1.3925e-02i 
 +6.7737e-03-7.2145e-03i  +1.0060e+00-1.5730e-02i  +2.6438e-02+1.9016e-02i 
 +2.0519e-02-1.7653e-03i  +3.6390e-03+1.4921e-03i  +9.9862e-01-4.0989e-03i 

 cos(A)**2-sin(A)**2 :
 +1.0313e+00+6.2536e-03i  -3.1624e-03+2.5763e-02i  -1.5574e-02-1.3925e-02i 
 +6.7737e-03-7.2145e-03i  +1.0060e+00-1.5730e-02i  +2.6438e-02+1.9016e-02i 
 +2.0519e-02-1.7653e-03i  +3.6390e-03+1.4921e-03i  +9.9862e-01-4.0989e-03i 


 Press   return to continue
 Press X return to stop