Mathc matrices/h05b

Un livre de Wikilivres.
Aller à la navigation Aller à la recherche


Sommaire


Installer et compiler ces fichiers dans votre répertoire de travail. Installer ce fichier dans votre répertoire de travail.

Crystal Clear mimetype source h.png vdotu.h
'
/* ------------------------------------ */
/*  Save as :   vdotu.h                */
/* ------------------------------------ */

/* ------------------------------------ */
/* ------- u.v = v^t u ---------------- */
/* ------------------------------------ */
double dot_uv_R(
double **u,
double **v
)
{
double **vt  = i_mR(R1, rsize_R(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 norm_uv_R(
double **u
)
{
    return (pow(dot_uv_R(u,u),(1./2.)));
}
/* ------------------------------------  */
/* -- d(u,v) = ||u-v||                -- */
/* ------------------------------------  */
double dist_uv_R(
double **u,
double **v
)
{
double **umnsv  =  i_mR(rsize_R(v),csize_R(v));

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

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


    smul_mR(s,v,projuv);
    
  return(projuv); 
}
/* ------------------------------------ */
/* ------------------------------------ */
double **orth_uv_mR(
double **U,
double **Orth
)
{
double **Vect_U      = i_mR(rsize_R(U),C1);
double **Vect_T      = i_mR(rsize_R(U),C1);
double **Vect_Orth   = i_mR(rsize_R(U),C1);
double **Vect_projuv = i_mR(rsize_R(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); 
		                
		    proj_uv_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 **Normalize_uv_mR(
double **Orth
)
{
double **u  = i_mR(rsize_R(Orth),C1);
double **Nu = i_mR(rsize_R(Orth),C1);	

int c_Orth   = csize_R(Orth);
int c;

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

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

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

int c_Orth   = csize_R(Orth);
int c;

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

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

  f_mR(u);  
  f_mR(Nu);  
  
 return(Norm);
}
/* ------------------------------------ */
/* ------------------------------------ */
void QR_uv_mR(
double **U,
double **Q,
double **R)
{
double **Q_T = i_mR(csize_R(Q),rsize_R(Q));

  orth_uv_mR(U,Q);
  Normalize_uv_mR(Q);  
  transpose_mR(Q,Q_T);
  mul_mR(Q_T,U, R);  
}
/* ------------------------------------ */
/* ------------------------------------ */
void r_Q_uv_mR(
double **U,
double **Q)
{
  orth_uv_mR(U,Q);
  Normalize_uv_mR(Q);  
}
/* ------------------------------------ */
double **eigs_value_uv_mR(
double **A,
double **EigsValue
)
{
int r = rsize_R(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<LOOP_EIGSUV;i++)
     { 
	   QR_uv_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 **svd_sqrt_uv_mR(
double **A,
double **sqrtA
)
{
int r;
int c;

  dif_sizes_mR(A,sqrtA,"svd_sqrt_mR();","(A or sqrtA)");

  for   ( r=R1; r<A[R_SIZE][C0]; r++)
    for ( c=C1; c<A[C_SIZE][C0]; c++)

            sqrtA[r][c] = sqrt(A[r][c]);
            
 return(sqrtA);
}
/* ------------------------------------ */
double **svd_inv_D_uv_mR(
double **S,
double **invS
)
{
int r;
int c;

  dif_sizes_mR(S,invS,    "svd_inv_D_uv_mR();","(S,invS)");

  for   ( r=R1; r<S[R_SIZE][C0]; r++)
    for ( c=C1; c<S[C_SIZE][C0]; c++)

         if(r==c)
            {
			 if(fabs(S[r][c]) < ERROR_E) 
			             invS[r][c] = 0;
             else        invS[r][c] = (1./S[r][c]);
		    }
         else
            invS[r][c] = 0;
         
 return(invS);
}
/* ------------------------------------ */
/* ------------------------------------ */
/* ------------------------------------ */
double **svd_value_uv_mR(
double **A,
double **SvdValue)
{
int rA = rsize_R(A);
int cA = csize_R(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);   
  eigs_value_uv_mR(S,SvdValue_2);
  svd_sqrt_uv_mR(SvdValue_2,SvdValue);
  
  f_mR(A_T);
  f_mR(S);
  f_mR(SvdValue_2);
  
  return(SvdValue);    
}
/* ------------------------------------ */
/* ------------------------------------ */


Déclaration des fichiers h.