/* ------------------------------------ */
/*  Save as :   vbdotu.h                */
/* ------------------------------------ */

/* ------------------------------------ */
/* ------- u.v = v^t u ---------------- */
/* ------------------------------------ */
double dotuv_R(
double **u,
double **v
)
{
double **vt  = i_mR(R1, rsize_mR(v));
double **vtu = i_mR(R1, C1);

double udotv;
  
  transpose_mR(v,vt);
        mul_mR(vt,u,vtu);

  udotv = vtu[R1][C1];
  
  f_mR(vt);
  f_mR(vtu);
  
  return (udotv);
}
/* ------------------------------------ */
/* ------  ||u|| =  (u.u)^(1/2) ------- */
/* ------------------------------------ */
double normuv_R(
double **u
)
{
    return (pow(dotuv_R(u,u),(1./2.)));
}
/* ------------------------------------  */
/* -- d(u,v) = ||u-v||                -- */
/* ------------------------------------  */
double distuv_R(
double **u,
double **v
)
{
double **umnsv  =  i_mR(rsize_mR(v),csize_mR(v));

double d;
  
  sub_mR(u,v,umnsv);

  d = normuv_R(umnsv);
  
  f_mR(umnsv);
  
  return (d);
}
/* ------------------------------------ */
/*        (<u,v>/||v||^2)  * v          */
/* ------------------------------------ */
double **projuv_mR(
double **u,
double **v,
double **projuv
)
{
           /* <u,v>     /  ||v||^2  */
double s = dotuv_R(u,v) / dotuv_R(v,v);


    smul_mR(s,v,projuv);
    
  return(projuv); 
}
/* ------------------------------------ */
/* ------------------------------------ */
double **orthuv_mR(
double **U,
double **Orth
)
{
double **Vect_U      = i_mR(rsize_mR(U),C1);
double **Vect_T      = i_mR(rsize_mR(U),C1);
double **Vect_Orth   = i_mR(rsize_mR(U),C1);
double **Vect_projuv = i_mR(rsize_mR(U),C1);

int c_U;	
int c_Orth;

  c_c_mR(U,C1, Orth,C1);

  for( c_U = C2; c_U < U[C_SIZE][C0]; c_U++)
      {
		c_c_mR(U,c_U, Vect_U,C1);
		c_c_mR(U,c_U, Vect_T,C1);

		for( c_Orth = C1; c_Orth < c_U; c_Orth++)
          {		
		    c_c_mR(Orth,c_Orth, Vect_Orth,C1); 
		                
		    projuv_mR(Vect_U,Vect_Orth,Vect_projuv);  
		    sub_mR(Vect_T,Vect_projuv,Vect_Orth);
		    
		    c_mR(Vect_Orth,Vect_T);    
	      }
	      c_c_mR(Vect_T,C1, Orth,c_Orth);		
	  }
	  
  f_mR(Vect_U);  
  f_mR(Vect_T);   
  f_mR(Vect_Orth);  
  f_mR(Vect_projuv);  
         
 return(Orth);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **Normalizeuv_mR(
double **Orth
)
{
double **u  = i_mR(rsize_mR(Orth),C1);
double **Nu = i_mR(rsize_mR(Orth),C1);	

int c_Orth   = csize_mR(Orth);
int c;

		for( c = C1; c <= c_Orth; c++)
          {		
		    c_c_mR(Orth,c, u,C1); 
		    		           
		    smul_mR(1./normuv_R(u),u,Nu);

            c_c_mR(Nu,C1, Orth,c); 
	      }

  f_mR(u);  
  f_mR(Nu);  
  
 return(Orth);
}
/* ------------------------------------ */
/* ------------------------------------ */
double **normalizeuv_mR(
double **Orth,
double **Norm
)
{
double **u  = i_mR(rsize_mR(Orth),C1);
double **Nu = i_mR(rsize_mR(Orth),C1);	

int c_Orth   = csize_mR(Orth);
int c;

		for( c = C1; c <= c_Orth; c++)
          {		
		    c_c_mR(Orth,c, u,C1); 
		    		           
		    smul_mR(1./normuv_R(u),u,Nu);

            c_c_mR(Nu,C1, Norm,c); 
	      }

  f_mR(u);  
  f_mR(Nu);  
  
 return(Norm);
}
/* ------------------------------------ */
/* The number of rows of U must be superior
   or equal of the number of columns of U */
/* ------------------------------------ */
void QRuv_mR(
double **U,
double **Q,
double **R)
{
double **Q_T = i_mR(csize_mR(Q),rsize_mR(Q));

  orthuv_mR(U,Q);
  Normalizeuv_mR(Q);  
  transpose_mR(Q,Q_T);
  mul_mR(Q_T,U, R);  
}
/* ------------------------------------ */
/* ------------------------------------ */
void r_Quv_mR(
double **U,
double **Q)
{
  orthuv_mR(U,Q);
  Normalizeuv_mR(Q);  
}
/* ------------------------------------ */
double **eigsuv_mR(
double **A,
double **EigsValue
)
{
int r = rsize_mR(A);

double **T      =      i_mR(r,r);
double **Q      =      i_mR(r,r);
double **R      =      i_mR(r,r);
int i;
int rc;

 c_mR(A,T);
 
 for(i=0;i<1000;i++)
     { 
	   QRuv_mR(T,Q,R);  
        mul_mR(R,Q,T);
     }
     
  for(rc=R1;rc<=r;rc++)
     { 
	   EigsValue[rc][C1] = T[rc][rc];
     }
     
  f_mR(T);
  f_mR(Q);
  f_mR(R); 
  
  return(EigsValue);    
}
/* ------------------------------------ */
double **svduv_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);

double **SvdValue_2 =  i_mR(cA,C1);

  transpose_mR(A,A_T);
  mul_mR(A_T,A, S);   
  eigsuv_mR(S,SvdValue_2);
  sqrt_mR(SvdValue_2,SvdValue);
  
  f_mR(A_T);
  f_mR(S);
  f_mR(SvdValue_2);
  
  return(SvdValue);    
}
/* ------------------------------------ */
/* ------------------------------------ */