Mathc matrices/h12c

Bibliothèque

Installer ce fichier dans votre répertoire de travail.

	vdsvdcn.h

/* ------------------------------------ */
/*  Save as :   vdsvdcn.h               */
/* ------------------------------------ */
/* ------------------------------------ */
void svd_U_Cn_mR(
double **A,
double **U
)
{
int i;	

int r = rsize_R(A);
int c = csize_R(A);

double **A_T       =      i_mR(c,r);	
double **AA_T      =      i_mR(r,r);
	
double **EigsValue =      i_mR(r,C1);
double **eignV     =      i_mR(r,r);

double **Ab;
double **b         =   m0_mR(i_mR(r,C1));

double **Ide       =  eye_mR(i_mR(r,r));
double **sIde      =         i_mR(r,r);
double **AA_TmnssIde =      i_mR(r,r);

/* ------------------------------------ 
   eigs_value_AA_T_mR(A,  eignValue); if(c > r)
   eigs_value_AA_T_mR(A_T,eignValue); if(r > c)
   ------------------------------------ */
  eigs_value_AA_T_mR(A,EigsValue);         
                          
  transpose_mR(A,A_T) ;                                             
  mul_mR(A,A_T,AA_T);
  
  for(i = R1; i <= rsize_R(EigsValue); i++)
  {  
   smul_mR(EigsValue[i][C1],Ide,sIde);
   MmnsD_mR(AA_T,sIde,AA_TmnssIde);
   Ab = c_A_b_Ab_mR(AA_TmnssIde,b,i_Abr_Ac_bc_mR(r,r,C1));
   
   GJ_PP_FreeV_mR(Ab,eignV);
   
   c_c_mR(eignV,C2,U,i);  
     
  }    
   
  Normalize_mR(U);  

  f_mR(EigsValue);
  f_mR(eignV);
  
  f_mR(A_T);  
  f_mR(AA_T);

  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(AA_TmnssIde);
}

/* ------------------------------------ */
/* ------------------------------------ */
void svd_V_Cn_mR(
double **A,
double **V
)
{
int i;
int j;	

int r = rsize_R(A);
int c = csize_R(A);

double **A_T       =         i_mR(c,r);	
double **A_TA      =         i_mR(c,c);

double **EigsValue =         i_mR(c,C1);
double **eignV     =         i_mR(c,c);

double **Ab;
double **Ab0;
double **b         =   m0_mR(i_mR(c,C1));

double **Ide       =  eye_mR(i_mR(c,c));
double **sIde      =         i_mR(c,c);
double **A_TAmnssIde =         i_mR(c,c);

/* ------------------------------------ 
   eigs_value_AA_T_mR(A,  eignValue); if(c > r)
   eigs_value_AA_T_mR(A_T,eignValue); if(r > c)
   ------------------------------------ */
   eigs_value_AA_T_mR(A,EigsValue);
   
  transpose_mR(A,A_T);
  
  mul_mR(A_T,A,A_TA);
    
  for(i = R1; i <= rsize_R(EigsValue); i++)
  { 
   if(EigsValue[i][C1])
     {   
      smul_mR(EigsValue[i][C1],Ide,sIde);
      MmnsD_mR(A_TA,sIde,A_TAmnssIde);
      
      Ab = i_Abr_Ac_bc_mR(rsize_R(A_TAmnssIde),
                     csize_R(A_TAmnssIde),
                     csize_R(b));
      
      c_A_b_Ab_mR(A_TAmnssIde,b, Ab); 
               
      GJ_PP_FreeV_mR(Ab,eignV);
      c_c_mR(eignV,C2,V,i);  
     } 
   else
    {
      smul_mR(EigsValue[i][C1],Ide,sIde);
      MmnsD_mR(A_TA,sIde,A_TAmnssIde);
           
      Ab0 = i_Abr_Ac_bc_mR(rsize_R(A_TAmnssIde),   
                      csize_R(A_TAmnssIde),       
                      (C1+c-r)                          
                      );
                           
      c_A_b_Ab_mR(A_TAmnssIde,b, Ab0);
      
      GJ_PP_FreeV_Value0_mR(Ab0,eignV);		     
      		
	 for(j=C2; i <= rsize_R(EigsValue); i++,j++)
	       c_c_mR(eignV,j,V,i);
    }
  }    
   
  Normalize_mR(V);  

  f_mR(EigsValue);
  f_mR(eignV);
   
  f_mR(A_TA);
  
  f_mR(Ab);
  f_mR(Ab0);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(A_TAmnssIde);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **pseudo_Cn_mR(
double **A,
double **Pinv,
double   factor_E
)
{
int r = rsize_R(A);
int c = csize_R(A);

double **a   = smul_mR(factor_E,A,i_mR(r,c));

double **U      = i_mR(r,r);
double **U_T    = i_mR(r,r);
double **V      = i_mR(c,c);
double **U_TA   = i_mR(r,c);
double **S      = i_mR(r,c);
double **S_T    = i_mR(c,r);
double **invS_T = i_mR(c,r); 
double **VS_T   = i_mR(c,r);
	
 svd_U_Cn_mR(a,U);   
 svd_V_Cn_mR(a,V);
 
/*  S = U_T * B * V   */
  transpose_mR(U,U_T); 
  mul_mR(U_T,A,U_TA);
  mul_mR(U_TA,V,S); 

/*  invS_T =        */   
  transpose_mR(S,S_T);
  inv_svd_mR(S_T,invS_T); 

/* PseudoInverse = V * invS_T * U_T */   
  mul_mR(V,invS_T,VS_T);
  mul_mR(VS_T,U_T,Pinv);  

 f_mR(a);   
 f_mR(U);
 f_mR(V);
 f_mR(U_T);  
 f_mR(U_TA);
 f_mR(S);
 f_mR(S_T);
 f_mR(invS_T);
 f_mR(VS_T); 
    
 return(Pinv);	
}
/* ------------------------------------ */
/* ------------------------------------ */

Déclaration des fichiers h.