Mathc matrices/h12c

Bibliothèque

Installer ce fichier dans votre répertoire de travail.

	vdsvdcn.h

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

  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(U);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_V_Cn_mR(
double **A,
double **V
)
{
int i;
int j;	

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

double **Ab;
double **Ab0;
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);
  
        mul_mR(tA,A,tAA);
  
        eye_mR(Ide);
    
 for(i=R1; i<= rsize_R(EValue); i++)
  { 
   if(EValue[i][C1])
    {   
      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,V,i);  
    } 
  else
   {
/* smul_mR(EValue[i][C1],Ide,sIde);          EValue[i][C1] = 0
    sub_mR(tAA,sIde,tAAsIde); */           

     Ab0 = c_A_b_Ab_mR(tAA,b, 
     i_Abr_Ac_bc_mR(c,c,(C1+c-r)));
      
    GJ_PP_FreeV_Value0_mR(Ab0,EV);		     
      		
   for(j=C2; i<=rsize_R(EValue); i++,j++)
	 
	     c_c_mR(EV,j,V,i);
    }
  }    
   
  Normalize_mR(V);  

  f_mR(EValue);
  f_mR(EV);
   
  f_mR(tAA);
  f_mR(AtA);
    
  f_mR(Ab);
  f_mR(Ab0);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(tAAsIde);

  return(V);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **Pinv_Cn_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 **U           = i_mR(r,r);
double **tU          = i_mR(r,r);
double **V           = i_mR(c,c);
double **tUA         = i_mR(r,c);
double **EValue      = i_mR(r,c);
double **tEValue     = i_mR(c,r);
double **invtEValue  = i_mR(c,r); 
double **VinvtEValue = i_mR(c,r);
	
 svd_U_Cn_mR(FA,U);   
 svd_V_Cn_mR(FA,V);
 
/*  EValue = tU  A  V   */
  transpose_mR(U,tU); 
        mul_mR(tU,A,tUA);
        mul_mR(tUA,V,EValue); 

  transpose_mR(EValue,tEValue);
    inv_svd_mR(tEValue,invtEValue); 

/* PseudoInverse = V invtEValue tU */   
  mul_mR(V,invtEValue,VinvtEValue);
  mul_mR(VinvtEValue,tU,Pinv);  

 f_mR(FA);   
 f_mR(U);
 f_mR(V);
 f_mR(tU);  
 f_mR(tUA);
 f_mR(EValue);
 f_mR(tEValue);
 f_mR(invtEValue);
 f_mR(VinvtEValue); 
    
 return(Pinv);	
}
/* ------------------------------------ */
/* ------------------------------------ */

Déclaration des fichiers h.