ORTHONORMALISATION DE GRAM-SCHMIDT

#ifndef _UTIL_H_
	#include "util.h"
#endif // _UTIL_H_

#ifndef _MATH_H_
	#include "math.h"
#endif // _MATH_H_

#ifndef _LIST_H_
	#include "list.h"
#endif // _LIST_H_

#ifndef _WINUTIL_H_
	#include "winutil.h"
#endif // _WINUTIL_H_

#include "conio.h"


//-------------------------------------------------
// OUTILS
//-------------------------------------------------
void printFloat(double x)
{
printf("%lf",x);
} // printFloat()

//-------------------------------------------------
void printVect(P_VECTND u)
{
int i;
printf("(");
printFloat(u->coord[0]);
for(i=1;i<u->dim;i++)
  {
  printf(" ,");
  printFloat(u->coord[i]);
  }
printf(" )\n");
} // printVect()

//-------------------------------------------------
// GRAM-SCHMIDT
//-------------------------------------------------
typedef double (*DOT_PRODUCT)(P_VECTND,P_VECTND);

//-------------------------------------------------
double dotProduct(DOT_PRODUCT phi,P_VECTND u,P_VECTND v)
{
return phi(u,v);
} // dotProduct()

//-------------------------------------------------
double norm(DOT_PRODUCT phi,P_VECTND u)
{
return sqrt(phi(u,u));
} // norm()

//-------------------------------------------------
PP_VECTND GramSchmidt(PP_VECTND vect,int nb,DOT_PRODUCT phi)
{
PP_VECTND res;
int       dim,i;
double    n;

AssertPointer(vect);
AssertPointer(phi);
Assert(nb >= 1);

dim = vect[0]->dim;
res = Malloc(P_VECTND,nb);
for(i=0;i<nb;i++)
  {
  res[i] = CreateVectND(dim);
  }

// le premier vecteur est juste normalise
n = norm(phi,vect[0]);
Assert(!CmpDouble0(n));
MulByScalarVectND(res[0],vect[0],1./n);
for(i=1;i<nb;i++)
  {
  int j;

  Assert(vect[i]->dim == dim);
  CopyVectND(res[i],vect[i]);

  // on hote les composantes suivants les vecteurs deja "orthonormalises"
  for(j=0;j<i;j++)
    {
    double k;
    k = dotProduct(phi,res[i],res[j]); // la composante est le produit scalaire ...
    LinCombinVectND(res[i],res[i],res[j],1.,-k);
    }
  
  // on normalise
  n = norm(phi,res[i]);
  Assert(!CmpDouble0(n));
  MulByScalarVectND(res[i],res[i],1./n);
  }

return res;
} //GramSchmidt()

//-------------------------------------------------
// MAIN
//-------------------------------------------------
double myDotProduct(P_VECTND u,P_VECTND v)
{
double x1,y1,z1;
double x2,y2,z2;

x1 = u->coord[0];
y1 = u->coord[1];
z1 = u->coord[2];

x2 = v->coord[0];
y2 = v->coord[1];
z2 = v->coord[2];

return x1*x2+y1*y2+z1*z2+0.5*(x1*y2+x2*y1+x1*z2+x2*z1+y1*z2+y2*z1);
} // myDotProduct()

//-------------------------------------------------
int main(int argc,char **argv)
{
InitializationLibUtil();
  {
  P_VECTND  vect[3];
  PP_VECTND res;

  vect[0] = CreateVectND(3);vect[0]->coord[0] = 1;vect[0]->coord[1] = -1;vect[0]->coord[2] = 0;
  vect[1] = CreateVectND(3);vect[1]->coord[0] = 1;vect[1]->coord[1] =  1;vect[1]->coord[2] = 0;
  vect[2] = CreateVectND(3);vect[2]->coord[0] = 1;vect[2]->coord[1] = -1;vect[2]->coord[2] = 1;

  res = GramSchmidt(vect,3,myDotProduct);

  printVect(vect[0]);
  printVect(vect[1]);
  printVect(vect[2]);
  printf("\nPar le procede d'orthonormalisation de Gram-Schmidt\n");
  printf("la famille de vecteurs ci-dessus devient :\n\n");
  printVect(res[0]);
  printVect(res[1]);
  printVect(res[2]);

  DeleteVectND(vect[0]);
  DeleteVectND(vect[1]);
  DeleteVectND(vect[2]);

  DeleteVectND(res[0]);
  DeleteVectND(res[1]);
  DeleteVectND(res[2]);

  Free(res);
  }
CloseLibUtil();
CheckingClosingLibUtil();

printf("\nfin du programme ...\n");
getch();
return 0 ;
} // main()