Mathc matrices/h12b

Bibliothèque

Installer ce fichier dans votre répertoire de travail.

	vdsvdrn.h

/* ------------------------------------ */
/*  Save as :   vdsvdrn.h              */
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_V_Rn_mR(
double **A,
double **V
)
{
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,V,i);  
     
  }    
   
  Normalize_mR(V);  

  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);

  return(V);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_U_Rn_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 **A_TA      =         i_mR(c,c);

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

double **Ab;
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++)
  {  
      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,U,i);  
  }    
   
  Normalize_mR(U);  

  f_mR(EigsValue);
  f_mR(eignV);
   
  f_mR(A_TA);
  
  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(A_TAmnssIde);

  return(U);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **pseudo_Rn_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 **a_T = transpose_mR(a,
                  i_mR(c,r));

double **U      = i_mR(r,c);
double **U_T    = i_mR(c,r);
double **V      = i_mR(c,c);

double **U_TA   = i_mR(c,c);

double **S      = i_mR(c,c);
double **S_T    = i_mR(c,c);

double **invS_T = i_mR(c,c);
double **VS_T   = i_mR(c,c); 

  svd_U_Rn_mR(a_T,U);    
  svd_V_Rn_mR(a_T,V); 
  
/*  S = U_T * A * 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(U);
  f_mR(U_T); 
  f_mR(V);  

  f_mR(a);  
  f_mR(a_T);
  f_mR(U_TA);
  
  f_mR(S); 
  f_mR(S_T);

  f_mR(invS_T);
  f_mR(VS_T);  
            
 return(Pinv);
}
/* ------------------------------------ */
/* ------------------------------------ */