Mathc complexes/h05d

Bibliothèque

Installer ce fichier dans votre répertoire de travail.

	wddotu.h

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

/* ------------------------------------ */
/* ------- u.v = conj(v)^t u ---------- */
/* ------------------------------------ */
nb_Z dot_Z(
double **u,
double **v
)
{
double **vt  = i_mZ(R1, rsize_Z(v));
double **vtu = i_mZ(R1, C1);

nb_Z   udotv;
  
  ctranspose_mZ(v,vt);
         mul_mZ(vt,u,vtu);

  udotv = i_Z( vtu[R1][C1],vtu[R1][C1+C1]);
  
  f_mZ(vt);
  f_mZ(vtu);
  
  return (udotv);
}

/* ------------------------------------ */
/* ------  ||u|| =  (u.u)^(1/2) ------- */
/* ------------------------------------ */
double norm_Z(
double **u
)
{
double x = Re_Z(dot_Z(u,u));
	
    return (pow(x,(1./2.)));
}

/* ------------------------------------  */
/* -- d(u,v) = ||u-v||                -- */
/* ------------------------------------  */
double dist_Z(
double **u,
double **v
)
{
double **umnsv  =  i_mZ(rsize_Z(v),csize_Z(v));

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

  d = norm_Z(umnsv);
  
  f_mZ(umnsv);
  
  return (d);
}

/* ------------------------------------ */
/* ------- (<u,v>/||v||^2).v ---------- */
/* ------------------------------------ */
double **proj_mZ(
double **u,
double **v,
double **projuv
)
{               /* z = <u,v> /  ||v||^2 */
nb_Z  z = div_Z(dot_Z(u,v),dot_Z(v,v));            
                  
  return(zmul_mZ(z,v,projuv)); /* (<u,v>/||v||^2).v */
}

/* ------------------------------------ */
/* ------- u.v = conj(v)^t u ---------- */
/* ------------------------------------ */
double **orth_mZ(
double **U,
double **Orth
)
{
double **Vect_U      = i_mZ(rsize_Z(U),C1);
double **Vect_T      = i_mZ(rsize_Z(U),C1);
double **Vect_Orth   = i_mZ(rsize_Z(U),C1);
double **Vect_projuv = i_mZ(rsize_Z(U),C1);

int c_U;	
int c_Orth;

  c_c_mZ(U,C1, Orth,C1);

  for( c_U = C2; c_U <= csize_Z(U); c_U++)
      {
		c_c_mZ(U,c_U, Vect_U,C1);
		c_c_mZ(U,c_U, Vect_T,C1);

		for( c_Orth = C1; c_Orth < c_U; c_Orth++)
          {		
		    c_c_mZ(Orth,c_Orth, Vect_Orth,C1); 
		 		                
		    proj_mZ(Vect_U,Vect_Orth,Vect_projuv);  
		    sub_mZ(Vect_T,Vect_projuv,Vect_Orth);
		    
		    c_mZ(Vect_Orth,Vect_T);    
	      }
	      c_c_mZ(Vect_T,C1, Orth,c_Orth);		
	  }
	  
  f_mZ(Vect_U);  
  f_mZ(Vect_T);   
  f_mZ(Vect_Orth);  
  f_mZ(Vect_projuv);  
         
 return(Orth);
}

/* ------------------------------------ */
/* ------- u.v = conj(v)^t u ---------- */
/* ------------------------------------ */
double **Normalize_mZ(
double **Orth
)
{
double **u  = i_mZ(rsize_Z(Orth),C1);
double **Nu = i_mZ(rsize_Z(Orth),C1);	

int c_Orth   = csize_Z(Orth);
int c;

		for( c = C1; c <= c_Orth; c++)
          {		
		    c_c_mZ(Orth,c, u,C1); 		    		           
		    smul_mZ(1./norm_Z(u),u,Nu);
            c_c_mZ(Nu,C1, Orth,c); 
	      }

  f_mZ(u);  
  f_mZ(Nu);  
  
 return(Orth);
}
/* ------------------------------------ */
/* ------- u.v = conj(v)^t u ---------- */
/* ------------------------------------ */
double **Normalize2_mZ(
double **Orth,
double **Norm
)
{
double **u  = i_mZ(rsize_Z(Orth),C1);
double **Nu = i_mZ(rsize_Z(Orth),C1);	

int c_Orth   = csize_Z(Orth);
int c;

		for( c = C1; c <= c_Orth; c++)
          {		
		    c_c_mZ(Orth,c, u,C1); 		    		           
		    smul_mZ(1./norm_Z(u),u,Nu);
            c_c_mZ(Nu,C1, Norm,c); 
	      }

  f_mZ(u);  
  f_mZ(Nu);  
  
 return(Norm);
}
/* ------------------------------------ */
/* ------- u.v = conj(v)^t u ---------- */
/* ------------------------------------ */
void QR_mZ(
double **U,
double **Q,
double **R)
{
double **Q_T = i_mZ(csize_Z(Q),rsize_Z(Q));

  orth_mZ(U,Q);
  Normalize_mZ(Q);  
  ctranspose_mZ(Q,Q_T);
  mul_mZ(Q_T,U, R);
  
  f_mZ(Q_T);    
}
/* ------------------------------------ */
double **eigs_mZ(
double **A,
double **EigsValue
)
{
int r = rsize_Z(A);

double **T      =      i_mZ(r,r);
double **Q      =      i_mZ(r,r);
double **R      =      i_mZ(r,r);
int i;
int re, ce;

 c_mZ(A,T);
 
 for(i=0;i<1000;i++)
     { 
	   QR_mZ(T,Q,R);  
        mul_mZ(R,Q,T);
     }
     
  for(re=R1,ce=C1;re<=r;re++,ce+=C2)
     { 
	   EigsValue[re][C1] = T[re][ce];
	   EigsValue[re][C2] = T[re][ce+C1];
     }
     
  f_mZ(T);
  f_mZ(Q);
  f_mZ(R); 
  
  return(EigsValue);    
}
/* ------------------------------------ */
double **svd_mZ(
double **A,
double **SvdValue)
{
int rA = rsize_Z(A);
int cA = csize_Z(A);

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

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

  ctranspose_mZ(A,A_T);
  mul_mZ(A_T,A, S);   
  eigs_mZ(S,SvdValue_2);
  sqrt_svd_mZ(SvdValue_2,SvdValue);
  
  f_mZ(A_T);
  f_mZ(S);
  f_mZ(SvdValue_2);
  
  return(SvdValue);    
}
/* ------------------------------------ */
/* ------------------------------------ */

Déclaration des fichiers h.