Mathc matrices/h20a

Application

Installer ce fichier dans votre répertoire de travail.

	dot_au.h

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

/* ------------------------------------ */
/* ------- <u.v> = v^t A^t A u  ------- */
/* ------------------------------------ */
double dot_Au_R(
double **A,
double **u,
double **v
)
{
double **At       = i_mR(csize_R(A), rsize_R(A));
double **vt       = i_mR(R1, rsize_R(v));
double **T_vtAt   = i_mR(R1, rsize_R(v));
double **T_vtAtA  = i_mR(R1, rsize_R(v));
double **T_vtAtAu = i_mR(R1, C1);

double dotAuAv;

  transpose_mR(A,At);
  transpose_mR(v,vt);
  
        mul_mR(vt,At,T_vtAt);
        mul_mR(T_vtAt,A,T_vtAtA);
        mul_mR(T_vtAtA,u,T_vtAtAu);
        
  dotAuAv = T_vtAtAu[R1][C1];
  
  f_mR(At);
  f_mR(vt);
  f_mR(T_vtAt);
  f_mR(T_vtAtA);
  f_mR(T_vtAtAu);
  
  return (dotAuAv);
}
/* ------------------------------------ */
/* ------  ||u|| =  <u.u>^(1/2) ------- */
/* ------------------------------------ */
double norm_Au_R(
double **A,
double **u
)
{
    return ( pow( dot_Au_R(A,u,u),(1./2.)) );
}
/* ------------------------------------  */
/* ------ d(u,v) = ||u-v||^(1/2) ------- */
/* ------------------------------------  */
double dist_Au_R(
double **A,
double **u,
double **v
)
{
double **umnsv  =  i_mR(rsize_R(v),csize_R(v));

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

/*           d(u,v) = ||u-v||^(1/2)       */

  d = norm_Au_R(A,umnsv);
  
  f_mR(umnsv);
  
  return (d);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t A^t A u  ------- */
/* ------------------------------------ */
double **proj_Au_mR(
double **A,
double **u,
double **v,
double **projuv
)
{
           /* <u,v>     /  ||v||^2  */
double s = dot_Au_R(A,u,v) / dot_Au_R(A,v,v);

/* (<u,v>/||v||^2) v */
    smul_mR(s,v,projuv);
    
  return(projuv);   
}
/* ------------------------------------ */
/* ------- <u.v> = v^t A^t A u  ------- */
/* ------------------------------------ */
double **orth_Au_mR(
double **A,
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_Au_mR(A,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);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t A^t A u  ------- */
/* ------------------------------------ */
double **Normalize_Au_mR(
double **A,
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_Au_R(A,u),u,Nu);

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

  f_mR(u);  
  f_mR(Nu);  
  
 return(Orth);
}
/* ------------------------------------ */
/* ------- <u.v> = v^t A^t A u  ------- */
/* ------------------------------------ */
double **mul_Au_mR(
double **A,
double **U,
double **V,
double **UV
)
{
double **u     = i_mR(rsize_R(U),C1);
double **v     = i_mR(rsize_R(U),C1);

int c1;
int c2;

  for   ( c1=C1; c1<U[C_SIZE][C0]; c1++)
    for ( c2=C1; c2<U[C_SIZE][C0]; c2++) 
         {
		   c_c_mR(U,c1, u,C1);
		   c_c_mR(V,c2, v,C1);
		   		   
		   UV[c1][c2] = dot_Au_R(A,u,v);
		 }

  f_mR(u);
  f_mR(v); 	
  
 return(UV);
}
/* ------------------------------------ */
/* ---------- u.v = v^t u ------------- */
/* ------------------------------------ */
void QR_Au_mR(
double **A,
double **U,
double **Q,
double **R)
{
  orth_Au_mR(A,U,Q);
  Normalize_Au_mR(A,Q);  
  mul_Au_mR(A,Q,U,R);  
}
/* ------------------------------------ */
/* ------------------------------------ */
/* ------------------------------------ */
double **eigs_value_Au_mR(
double **A,
double **U,
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(U,T);
 
 for(i=0;i<LOOP_EIGSUV;i++)
     { 
	    QR_Au_mR(A,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 **xeigs_value_Au_mR(
double **A,
double **U,
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(U,T);
 
do
{ 
	    QR_Au_mR(A,T,Q,R);  
        mul_mR(R,Q,T);
        clrscrn();
        pE_mR(T,S5,P3,C6); 

} while(stop_w());
        getchar();
  
  for(rc=R1;rc<=r;rc++)
     { 
	   EigsValue[rc][C1] = T[rc][rc];
     }
     
  f_mR(T);
  f_mR(Q);
  f_mR(R); 
  
  return(EigsValue);    
}*/
/* ------------------------------------ */
/* ------------------------------------ */

Déclaration des fichiers h.