/* ------------------------------------ */
/*  Save as :   vgjsvdc1.h              */
/* ------------------------------------ */
/* ------------------------------------ */
double **svduvpow2_mR(
double **A,
double **SvdValue)
{
int rA = rsize_mR(A);
int cA = csize_mR(A);

double **A_T    =      i_mR(cA,rA);
double **S      =      i_mR(cA,cA);

  transpose_mR(A,A_T);
  mul_mR(A_T,A, S);   
  eigsuv_mR(S,SvdValue);
  
  f_mR(A_T);
  f_mR(S);
  
  return(SvdValue);    
}
/* ------------------------------------ */
/* ------------------------------------ */
void svdU_mR(
double **a,
double **EignVector
)
{
int i;	

int r = rsize_mR(a);
int c = csize_mR(a);
	
double **A         =      i_mR(r,r);
double **A_T       =      i_mR(c,r);	
	
double **EigsValue;
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 **sIde_A    =         i_mR(r,r);

/* -------------------------------------------- */ 
  if( r>=c)
     {
	  EigsValue   = i_mR(c,C1);
      svduvpow2_mR(a,EigsValue);
     }
  else
     {
	  EigsValue   = i_mR(r,C1);
      svduvpow2_mR(transpose_mR(a,A_T),EigsValue);
     }
/* -------------------------------------------- */ 

  mul_mR(a,A_T,A);
  
  for(i = R1; i <= rsize_mR(EigsValue); i++)
  {  
   smul_mR(c_coef_R(EigsValue,i,C1),Ide,sIde);
   sub_mR(sIde,A,sIde_A);
   Ab = c_A_b_Ab_mR(sIde_A,b,i_AbR0_mR(r,r,C1));
   
   gj_freevariable_eignvector_mR(Ab,eignV);
   
   c_c_mR(eignV,C2,EignVector,i);  
     
  }    
   
  Normalizeuv_mR(EignVector);  

  f_mR(EigsValue);
  f_mR(eignV);
  
  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(sIde_A);
}

/* ------------------------------------ */
/* ------------------------------------ */

void svdV_mR(
double **a,
double **EignVector
)
{
int i;	

int r = rsize_mR(a);
int c = csize_mR(a);
	
double **A         =      i_mR(c,c);
double **A_T       =      i_mR(c,r);	
	
double **T_EigsValue;
double **EigsValue =      i_mR(c,C1);
double **eignV     =      i_mR(c,c);

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

/* -------------------------------------------- */ 
  if( r>=c)
     {
	  T_EigsValue   = i_mR(c,C1);
      svduvpow2_mR(a,T_EigsValue);
     }
  else
     {
	  T_EigsValue   = i_mR(r,C1);
      svduvpow2_mR(transpose_mR(a,A_T),T_EigsValue);
     }
/* -------------------------------------------- */ 
  
  mul_mR(A_T,a,A);
  
  c_mR(T_EigsValue,EigsValue);
  
  for(i = R1; i <= rsize_mR(EigsValue); i++)
  {  
   smul_mR(c_coef_R(EigsValue,i,C1),Ide,sIde);
   sub_mR(sIde,A,sIde_A);
   Ab = c_A_b_Ab_mR(sIde_A,b,i_AbR0_mR(c,c,C1));
   
   gj_freevariable_eignvector_mR(Ab,eignV);
   
   c_c_mR(eignV,C2,EignVector,i);  
  }    
   
  Normalizeuv_mR(EignVector);  

  f_mR(T_EigsValue);
  f_mR(EigsValue);
  f_mR(eignV);
  
  f_mR(Ab);
  f_mR(b);
  
  f_mR(Ide);
  f_mR(sIde);
  f_mR(sIde_A);
}
/* ------------------------------------ */
/* ------------------------------------ */