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 **tA      = i_mR(c,r);	
double **AtA     = i_mR(r,r);
	
double **EValue  = i_mR(r,C1);
double **EV      = i_mR(r,r);

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

double **Ide     = i_mR(r,r);
double **sIde    = i_mR(r,r);
double **AtAsIde = i_mR(r,r);

  transpose_mR(A,tA) ;                                             
        mul_mR(A,tA,AtA);
       eigs_mR(AtA,EValue);
       
       eye_mR(Ide);
  
 for(i=R1; i<=rsize_R(EValue); i++)
  {  
    smul_mR(EValue[i][C1],Ide,sIde);
     sub_mR(AtA,sIde,AtAsIde);
           
   Ab = c_A_b_Ab_mR(AtAsIde,b,
         i_Abr_Ac_bc_mR(r,r,C1)); 
         
     GJ_PP_FreeV_mR(Ab,EV);
       
         c_c_mR(EV,C2,V,i);  
  }    
   
  Normalize_mR(V);  

  f_mR(EValue);
  f_mR(EV);
  
  f_mR(tA);  
  f_mR(AtA);

  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(AtAsIde);

  return(V);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_U_Rn_mR(
double **A,
double **U
)
{
int i;

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

double **tA      = i_mR(c,r);	
double **tAA     = i_mR(c,c);
double **AtA     = i_mR(r,r);

double **EValue  = i_mR(r,C1);
double **EV      = i_mR(c,r);

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

double **Ide     = i_mR(c,c);
double **sIde    = i_mR(c,c);
double **tAAsIde = i_mR(c,c);

  transpose_mR(A,tA) ;                                             
        mul_mR(A,tA,AtA);
       eigs_mR(      AtA,EValue);

       eye_mR(Ide);
  
  mul_mR(tA,A,tAA);
    
 for(i=R1; i<=rsize_R(EValue); i++)
  {  
    smul_mR(EValue[i][C1],Ide,sIde);
     sub_mR(tAA,sIde,tAAsIde);
      
  Ab = c_A_b_Ab_mR(tAAsIde,b, 
       i_Abr_Ac_bc_mR(c,c,C1));  
         
 GJ_PP_FreeV_mR(Ab,EV);
 
   c_c_mR(   EV,C2,U,i);  
  }    
   
  Normalize_mR(U);  

  f_mR(EValue);
  f_mR(EV);

  f_mR(tA);    
  f_mR(tAA);
  f_mR(AtA);
    
  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(tAAsIde);

  return(U);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **Pinv_Rn_mR(
double **A,
double **Pinv,
double   FACTOR
)
{
int r = rsize_R(A);
int c = csize_R(A);

double **FA         = smul_mR(FACTOR,A, 
                      i_mR(r,c));
                  
double **tFA        = transpose_mR(FA, 
                      i_mR(c,r));

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

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

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

double **invEValue  = i_mR(c,c);
double **VinvEValue = i_mR(c,c); 

  svd_U_Rn_mR(tFA,U);    
  svd_V_Rn_mR(tFA,V); 
  
/*  EValue = tU A V  */ 
 transpose_mR(U,tU);
       mul_mR(tU,A,tUA);
       mul_mR(tUA,V,EValue);	

   inv_svd_mR(EValue,invEValue);

/* PseudoInverse = V invEValue tU */  
   mul_mR(V,invEValue,VinvEValue);
   mul_mR(VinvEValue,tU,Pinv);    

  f_mR(U);
  f_mR(tU); 
  f_mR(V);  

  f_mR(FA);  
  f_mR(tFA);
  f_mR(tUA);
  
  f_mR(EValue); 

  f_mR(invEValue);
  f_mR(VinvEValue);  
            
 return(Pinv);
}
/* ------------------------------------ */
/* ------------------------------------ */