\documentclass[10pt]{article}
\usepackage[latin1]{inputenc}
%\usepackage{epsfig}
\usepackage{hyperref}
\usepackage{color}
%\usepackage{srcltx}\SRCOKfalse
\usepackage{amsfonts}
%\usepackage{keyboard}
\usepackage{FF}
\usepackage[frenchb]{babel}
\usepackage{mathptmx,times}
\definecolor{orange}{rgb}{1,0.5,0}
\definecolor{darkblue}{rgb}{0,0,0.4}
\definecolor{darkgreen}{rgb}{0,0.4,0}
\definecolor{rose}{rgb}{1,0.25,0.25}
\definecolor{cyan}{rgb}{0,0.6,0.6}
\definecolor{magenta}{rgb}{0.6,0,0.6}
\definecolor{lightgreen}{rgb}{0.5,1,0.5}
%%   for the frames/page numbers etc.
\def\textRed{\color{red}}
\def\textBlack{\color{black}}
\def\textWhite{\color{white}}
\def\textBlue{\color{blue}}
\def\textRose{\color{rose}}
\def\textYellow{\color{yellow}}
\def\textGreen{\color{green}}
\def\textCyan{\color{cyan}}
\def\textlg{\color{lightgreen}}
\def\Blue#1{{\color{blue}{#1}}}
\def\Black#1{{\color{black}{#1}}}
\def\White#1{{\color{white}{#1}}}
\def\Red#1{{\color{red}{#1}}}
\def\Green#1{{\color{darkgreen}{#1}}}
\def\Yellow#1{{\color{yellow}{#1}}}
\def\Orange#1{{\color{orange}{#1}}}
\def\Rose#1{{\color{rose}{#1}}}
\def\Cyan#1{{\color{cyan}{#1}}}
\def\Magenta#1{{\color{magenta}{#1}}}
\def\LightGreen#1{{\color{lightgreen}{#1}}}
%\usepackage{french}
\def\R{\mathrm{I\!R}}
\def\Rtwo{\mathrm{I\!R}^2}
\def\Th{\mathcal{T}_h}
\def\Vh{V_{h}}\def\Voh{V_{0h}}
\def\IR{\hbox{${\rm I\!R}$}}   
\def\ftpsrc{ftp\string://ann.jussieu.inria.fr/~hecht/ftp/cpp}
\def\ftpurl#1{~\url{\ftpsrc/#1}}
\newtheorem{remarque}{\sc Remarque}
\newenvironment{rem}{\begin{remarque}\rm}{\end{remarque}}

\topmargin -1.54cm      %marge en haut a 2cm %% on monte de 1.65 cm car dvi2ps marche en letter 11in
\oddsidemargin 0cm   %marge a 2 cm
\evensidemargin 0cm  %marge a 2 cm
\newcommand{\indente}{\hbox to \parindent {\hss}}
\parindent 0cm
\headsep 0.5cm
\topskip .5cm 
\footskip 1.5cm
\headheight 1.0cm
\textwidth  16.5cm
%  \parindent 0cm
\textheight 24cm
\author{Frédéric Hecht \\ Université Paris 6/ INRIA \\~\\ 
\url{mailto\string:Frederic.hecht@upmc.fr} \\
\url{http\string://www.ann.jussieu.fr/~hecht}}
\newtheorem{listing}{\bf Listing}
\newenvironment{listg}[1] %
{\par\medskip\begin{listing}\vbox{\hfill (#1)\\
\hbox{\rule[1ex]{\textwidth}{.05cm}}}
 }% 
{\end{listing}\hbox{\rule[1ex]{\textwidth}{.05cm}}
}
\newtheorem{Algo}{Algorithme}
\newenvironment{algo}{\begin{Algo}}{\end{Algo}}

\title{\Red{Des petites classes C++ pour tableaux avec des opérations 
vectorielles}}
\begin{document}
    \maketitle
    \tableofcontents
    
\section{Introduction}
\def\K{\mathbb{K}}
Dans ce chapitre nous allons nous intéresser à la construction de 
classes pour modéliser des tableaux à 1,2 où 3 indices,
soit $\K^n$,$\K^{n,m}$, et
$\K^{n,m,k}$ où 
$\K^{n,m}$ est l'espace matrice $n\times m$ à coefficient dans $\K$ et ou
$\K^{n,m,k}$ est l'espace de tenseurs d'ordre 3  à coefficient dans $\K$.

Donc nous allons pleinement utiliser les patrons \ttCC{<<template>>}  pour
construire ces  trois classes, où $\K$ sera le type paramétré.
 Nous n'utilisons pas les classes  \ttCC{valarray} du \texttt{C++} qui
 ne fonctionne pas encore très bien,  de plus les classes définies dans 
 le package {\tt Blitz++} \url{http\string://www.oonumerics.org/blitz} sont trops
 lourde pour être utiliser pour faire l'enseignement.

Nous commencerons sur une version didactique, nous verrons
 la version complète  accessible sur  le web 
 
où dans le fichier \og tar compressé \fg \ftpurl{RNM.tar.gz}. 


\def\K{\mathbb{K}} Dans ce chapitre, nous allons nous int\'eresser \`a la
construction de classes pour mod\'eliser des tableaux \`a 1,2 ou 3 indices, soit
$\K^n$,$\K^{n, m}$, et $\K^{n,m,k}$ o\`u $\K^{n,m}$ est l'espace des matrices
$n\times m$ \`a coefficient dans $\K$ et ou $\K^{n,m,k}$ est l'espace des tenseurs
d'ordre 3 \`a coefficient dans $\K$.

Donc nous allons pleinement utiliser les patrons \og template \fg pour construire
ces trois classes.

Nous commencerons sur une version didactique, nous verrons la version compl\`ete
qui est dans \ftpurl{source/RNM} o\`u dans le fichier \og tar compress\'e \fg
\ftpurl{RNM.tar.gz}.

\section{Version simple}
Mais avant toute chose, me para\^{\i}t clair qu'un vecteur sera un classe 
qui contient au minimum la taille n, et un pointeur sur les valeurs. 
Que faut'il dans cette classe minimale (not\'e \texttt{A}).

\bFF
  typedef double K; // d\'efinition du corps
  @class A { @public:  // version 1  -------
       
    @int n;            //  la taille du vecteur 
    @K *v;             //  pointeur sur les  n valeurs  
    A() { @cerr <<" Pas de constructeur pas d\'efaut " << @endl; exit(1);} 
    A(int i) : n(i),v(@new @K[i]) { assert(v);} // constructeur
    ~A() {@delete [] v;}  // destructeur 
     K & @operator[](@int i) @const  {assert(i>=0 && i <n);@return v[i];}
 };
\eFF

Cette classe ne fonctionne pas car le constructeur par copie par d\'efaut
fait une copie bit \`a bit et donc le pointeur \texttt{v} est perdu, il faut donc
\'ecrire :


\bFF
  @class A { @public:  // version 2  -------       
    @int n;            //  la taille du vecteur 
    @K *v;             //  pointeur sur les  n valeurs  
    A() { @cerr <<" Pas de constructeur pas d\'efaut " << endl; exit(1);} 
    A(const A& a) :n(a.n),v(new K[a.n]) // constructeur par copie
         { operator=(a);} 
    A(int i) : n(i),v(@new @K[i]) { assert(v);} // constructeur
    A& @operator=(A &a) {  //  copie 
          assert(n==a.n);         
          @for(@int i=0;i<n;i++) v[i]=a.v[i];
          @return *this;}  
    ~A() {@delete [] v;}  // destructeur 
    K & @operator[](@int i) @const  {assert(i>=0 && i <n);@return v[i];}
  };
\eFF

Maintenant nous voulons ajouter les op\'erations vectorielles $+, -,*,..$.

\bFF
  @class A { @public:  // version 3  -------       
    @int n;            //  la taille du vecteur 
    @K *v;             //  pointeur sur les  n valeurs  
    A() { @cerr <<" Pas de constructeur pas d\'efaut " << endl; exit(1);} 
    A(const A& a) :n(a.n),v(new K[a.n]) { operator=(a);}
    A(int i) : n(i),v(@new @K[i]) { assert(v);} // constructeur
    A& @operator=(A &a) {assert(n==a.n); 
                        @for(@int i=0;i<n;i++) v[i]=a.v[i];
                        @return *@this}  
    ~A() {@delete [] v;}  // destructeur 
    K & @operator[](@int i) @const  {assert(i>=0 && i <n);@return v[i];}
    A operator+(const &a) @const ; // addition 
    A operator*(K &a) @const ;     // espace vectoriel \`a droite 
	
  @private:             //  constructeur priv\'e pour faire des optimisations	
    A(int i, K* p) : n(i),v(p){assert(v);}	
    @friend A @operator*(const K& a,@const A& ) ;  // multiplication \`a gauche 
	
  };
\eFF

Il
 faut faire attention dans les param\`etres d'entr\'ee et de sortie des op\'erateurs. 
Il est clair que l'on ne veut pas travailler sur des copies, mais la sortie est
forc\'ement un objet et non une r\'ef\'erence sur un objet car il faut allouer de la
m\'emoire dans l'op\'erateur et si l'on retourne une r\'ef\'erence aucun destructeur ne
sera appel\'e et donc cette m\'emoire ne sera jamais lib\'er\'ee. 


\bFF
//version avec avec une copie du tableau au niveau du return 
 A A::@operator+(@const A & a) @const {
  A b(n); assert(n == a.n);
  @for (@int i=0;i<n;i++) b.v[i]= v[i]+a.v[i];
  @return b;  // ici l'op\'erateur  A(const A& a) est appeler
 }
\eFF
 Pour des raisons
optimisation nous ajoutons  un nouveau constructeur \verb;A(int,K*); qui
\'evitera de faire une copie du tableau.
\bFF
//  --- version optimisée sans copie----
 A A::@operator+(@const A & a) @const {
  K  *b(@new K[n]); assert(n == a.n);
  @for (@int i=0;i<n;i++) b[i]= v[i]+a.[i];
  @return A(n,b);  // ici l'op\'erateur  A(n,K*)  ne fait pas de copie 
 }
\eFF

Pour la multiplication par un scalaire \`a droite on a:

\bFF
//  --- version optimisée ----
 A A::@operator*(@const K & a) @const {
  K  *b(@new K[n]); assert(n == a.n);
  for (@int i=0;i<n;i++) b[i]= v[i]*a;
  @return A(n,b);  // ici l'op\'erateur  A(n,K*)  ne fait pas de copie 
 }
\eFF

Pour la version \`a gauche, il faut d\'efinir \texttt{operator*} ext\'erieurement \`a la
classe car le terme de gauche n'est pas un vecteur.

\bFF
 A @operator*(const K & a,const T & c) {
  K  *b(@new K[n]); assert(n == a.n);
  @for (@int i=0;i<n;i++) b[i]= c[i]*a;
  @return A(n,b);  // attention c'est op\'erateur est priv\'e donc 
                  // cette fonction doit est ami ( @friend) de la classe 
 }
\eFF




Maintenant regardons ce qui est ex\'ecut\'e dans le cas d'une expression vectorielle.

\bFF
  int n=100000; 
  A a(n),b(n),c(n),d(n);
  ... //  initialisation des tableaux a,b,c 
  d = (a+2.*b)+c*2.0;
\eFF
voil\`a le pseudo code g\'en\'er\'e avec les 3 functions suivantes:
\texttt{add(a,b,ab)},\texttt{mulg(s,b,sb)},\texttt{muld(a,s,as)} , \texttt{copy(a,b)} o\`u
le dernier argument retourne le r\'esultat.
\bFF
  A a(n),b(n),c(n),d(n);
  A t1(n),t2(n),t3(n),t4(n);
  muld(2.,b),t1);     //  t1 = 2.*b 
  add(a,t1,t2);       //  t2 = a+2.*b
  mulg(c,2.,t3);      //  t3 = c*2.
  add(t2,t3,t4);         //  t4 = (a+2.*b)+c*2.0;
  copy(t4,d);         //  d = (a+2.*b)+c*2.0;
\eFF

Nous voyons que quatre tableaux interm\'ediaires sont cr\'e\'es ce qui est excessif. 
D'o\`u l'id\'ee de ne pas utiliser toutes ces possibilit\'es car le code g\'en\'er\'e sera
trop lent.

\begin{remarque} Il est toujours possible de cr\'eer des classes interm\'ediaires
pour des op\'erations pr\'ed\'efinies afin obtenir se code g\'en\'erer :
\bFF
  A a(n),b(n),c(n),d(n);
  for (int i;i<n;i++)
    d[i] = (a[i]+2.*b[i])+c[i]*2.0;
\eFF
 Ce cas me para\^{\i}t trop compliquer, mais nous pouvons optimiser raisonnablement
toutes les  combinaisons lin\'eaires \`a 2 termes.
\end{remarque}
 
Mais, il est toujours possible d'\'eviter ces probl\`emes en utilisant les
op\'erateurs \texttt{+=, -=, *=, /=} ce que donnerai de ce cas 
\bFF 
A a(n),b(n),c(n),d(n); 
@for (int i;i<n;i++) 
   d[i] = (a[i]+2.*b[i])+c[i]*2.0;
\eFF
 D'o\`u l'id\'ee de d\'ecouper les classes vecteurs en deux types de classe les
 classes de terminer par un \texttt{\_} sont des classes sans allocation (cf. 
 new), et les autres font appel \`a l'allocateur \texttt{new} et au d\'eallocateur
 \texttt{delete}.  De plus comme en fortran 90, il est souvent utile de voir une
 matrice comme un vecteur ou d'extraire une ligne ou une colonne ou m\^eme une
 sous-matrice.  Pour pouvoir faire tous cela comme en fortran 90, nous allons
 consid\`ere un tableau comme un nombre d'\'el\'ement $n$, un incr\'ement $s$ et un
 pointeur $v$ et du tableau suivant de m\^eme type afin de extraire des
 sous-tableau.



\section{Les classes tableaux}


  Nous d\'efinissons des classes tableaux \`a un, deux ou trois indices avec ou sans
  allocation.  Ces tableaux est d\'efini par un pointeur sur les valeurs et par la
  forme de l'indice (class \ttCC{ShapeOfArray}) qui donne la taille, le pas
  entre de valeur et le tableau suivant de m\^eme type (ces deux donn\'ees
  suppl\'ementaire permettent extraire des colonnes ou des lignes de matrice,
  \ldots).

Nous voulons faire les op\'erations classiques sur \ttCC{A, C,D} tableaux de type
\ttCC{KN<R>} suivante par exemples:

\bFF   
  A = B; A += B ; A -= B ; A = 1.0; A = 2*.C ; A /=C; A *=C;
  A = A+B; A = 2.0*C+ B; C = 3.*A-5*B; 
  R c = A[i]; 
  A[j] = c;
\eFF

Pour des raisons \'evidentes nous ne sommes pas all\'es plus loin que des
combinaisons lin\'eaires \`a plus de deux termes. 
Toutes ces op\'erations sont faites
sans aucune allocation et avec une seule boucle.

De plus nous avons d\'efini, les tableaux 1,2 ou 3 indices, il est possible
extraire une partie d'un tableau, une ligne ou une colonne.



\bFF
  tyepdef double R; 
  KNM<R> A(10,20); // un matrice
  . . . 
  KN_<R> L1(A(1,'.'); // la ligne 1 de la matrice A;
  KN<R> cL1(A(1,'.');  // copie de la ligne 1 de la matrice A;
  KN_<R> C2(A('.',2); // la colonne 2 de la matrice A;
  KN<R> cC2(A('.',2); // copie de la colonne 2 de la matrice A;
  KNM_<R> pA(FromTo(2,5),FromTo(3,7)); //  partie de la matrice A(2:5,3:7) 
  // vue comme un matrice 4x5
  
  KNM B(n,n);
  B(SubArray(n,0,n+1)) // le vecteur diagonal de B;
  KNM_ Bt(B.t());     //  la matrice transpose sans copie 
  
\eFF

Pour l'utilisation,  utiliser l'ordre \ttCC{\#include "RNM.hpp"}, et les 
flags de compilation {\tt -DCHECK\_KN} ou en  d\'efinisant  la variable 
du preprocesseur \texttt{cpp} du \texttt{C++}  avec l'ordre \ttCC{\#defined CHECK\_KN}, 
avant la ligne include.

Les d\'efinitions des classes sont faites dans 4 fichiers 
\texttt{RNM.hpp}, \texttt{RNM\_tpl.hpp}, \texttt{RNM\_op.hpp}, 
\texttt{RNM\_op.hpp}.
  
Pour plus de d\'etails voici un exemple d'utilisation assez complet.

\section{Exemple d'utilisation}
\bFF
@namespace std

#define CHECK_KN
#include "RNM.hpp"
#include "assert.h"


@using @namespace std;  
//  definition des 6 types de base des tableaux a 1,2 et 3 parametres
@typedef @double R;
@typedef KN<R> Rn;
@typedef KN_<R> Rn_;
@typedef KNM<R> Rnm;
@typedef KNM_<R> Rnm_;
@typedef KNMK<R> Rnmk;
@typedef KNMK_<R> Rnmk_;
R f(@int i){@return i;}    
R g(@int i){@return -i;}
@int main()
{
  
  @const @int n= 8;	
  cout << "Hello World, @this is RNM use!" << endl << endl;
  Rn a(n,f),b(n),c(n);
  b =a;
  c=5;
  b *= c; 
  
  cout << " a = " << (KN_<const_R>) a << endl;
  cout << " b = " <<  b << endl;

  //  les operations vectorielles 
  c = a + b;
  c = 5. *b + a;
  c = a + 5. *b;
  c = a - 5. *b;
  c = 10.*a - 5. *b;
  c = 10.*a + 5. *b;
  c += a + b;
  c += 5. *b + a;
  c += a + 5. *b;
  c += a - 5. *b;
  c += 10.*a - 5. *b;
  c += 10.*a + 5. *b;
  c -= a + b;
  c -= 5. *b + a;
  c -= a + 5. *b;
  c -= a - 5. *b;
  c -= 10.*a - 5. *b;
  c -= 10.*a + 5. *b;
  
  
  cout <<" c = " << c << endl;
  Rn u(20,f),v(20,g);   //  2 tableaux u,v de 20 
                        //initialiser avec $u_i=f(i),v_j=g(i)$
  Rnm A(n+2,n);           //  Une matrice $n+2\times n$
  
  @for (@int i=0;i<A.N();i++)        // ligne 
	 @for (@int j=0;j<A.M();j++) // colonne
	   A(i,j) = 10*i+j; 
	   
  cout << "A=" << A << endl;
  cout << "Ai3=A('.', 3 ) = " << A('.', 3 ) << endl;   // la colonne 3
  cout << "A1j=A( 1 ,'.') = " << A( 1 ,'.') << endl;   // la ligne  1
  Rn CopyAi3(A('.', 3 )); //  une copie de la colonne 3  
  cout << "CopyAi3 = " << CopyAi3;
  
  Rn_ Ai3(A('.', 3 )); //  la colonne 3 de la matrice 
  CopyAi3[3]=100;
  cout << CopyAi3 << endl;
  cout << Ai3 << endl;
  
  assert( & A(0,3) == & Ai3(0)); //  verification des adresses
   

  Rnm S(A(SubArray(3),SubArray(3))); //  sous matrice 3x3
  
  Rn_ Sii(S,SubArray(3,0,3+1));  // la diagonal de la matrice  sans copy
  
  cout << "S=  A(SubArray(3),SubArray(3)  = " << S <<endl;
  cout << "Sii     = " << Sii <<endl;
  b = 1;
 //  Rn Ab(n+2) = A*b; error
  Rn Ab(n+2);
  Ab = A*b;
  cout << " Ab = A*b  =" << Ab << endl;        
  
  Rn_ u10(u,SubArray(10,5));    //  la partie [5,5+10[ du tableau u     
  cout << "u10 " << u10 << endl;
  v(SubArray(10,5)) += u10;
  cout << " v = " << v << endl;
  cout << " u(SubArray(10))    "  << u(SubArray(10)) << endl;
  cout << " u(SubArray(10,5))  "  << u(SubArray(10,5)) << endl;
  cout << " u(SubArray(8,5,2))"  << u(SubArray(8,5,2)) 
       << endl;
  

  cout << " A(5,'.')[1] " << A(5,'.')[1] << " " << " A(5,1) = " 
       << A(5,1) << endl;
  cout << " A('.',5)(1) = "<< A('.',5)(1) <<  endl;
  cout << " A(SubArray(3,2),SubArray(2,4)) = " << endl;
  cout <<  A(SubArray(3,2),SubArray(2,4)) << endl;	
  A(SubArray(3,2),SubArray(2,4)) = -1;
  A(SubArray(3,2),SubArray(2,0)) = -2;
  cout << A << endl;
  
  Rnmk B(3,4,5);
  @for (@int i=0;i<B.N();i++)        // ligne 
	 @for (@int j=0;j<B.M();j++) // colonne
	     @for (@int k=0;k<B.K();k++) // colonne
	   B(i,j,k) = 100*i+10*j+k;       
  cout << " B      = " << B << endl;
  cout << " B(1  ,2  ,'.') " << B(1  ,2  ,'.') << endl;
  cout << " B(1  ,'.',3  ) " << B(1  ,'.',3  ) << endl;
  cout << " B('.',2  ,3  ) " << B('.',2  ,3  ) << endl;
  cout << " B(1  ,'.','.') " << B(1,'.'  ,'.') << endl;
  cout << " B('.',2  ,'.') " << B('.',2  ,'.') << endl;
  cout << " B('.','.',3  ) " << B('.','.',3  ) << endl;
  
  cout << " B(1:2,1:3,0:3)   = " 
       << B(SubArray(2,1),SubArray(3,1),SubArray(4,0)) << endl;
       
 // copie du sous tableaux
        
  Rnmk Bsub(B(FromTo(1,2),FromTo(1,3),FromTo(0,3)));       
  B(SubArray(2,1),SubArray(3,1),SubArray(4,0))  = -1;
  B(SubArray(2,1),SubArray(3,1),SubArray(4,0)) += -1;
  cout << " B      = " << B << endl;
  cout << Bsub << endl;
   
  @return 0;
}

\eFF
  
\section{Un resolution de syst\'eme lin\'eaire avec le gradient conjugu\'e}


   L'algorithme du gradient conjugu\'{e} pr\'{e}sent\'{e} dans cette section est utilis\'{e}

   pour r\'{e}soudre le syst\`{e}me

   lin\'{e}aire $ A x = b$,  o{\`u} $A$ est une matrice sym\'{e}trique positive

   $n\times n$.

   Cet algorithme est bas\'{e} sur la minimisation de la fonctionnelle

quadratique $E : \R^n \rightarrow \R$

   suivante :

   $$ E(x) = {1\over2} (A x,x) - (b,x), $$

   o{\`u} $(.,.)$ est le produit scalaire classique de $\R^n$.



\begin{Algo} Le gradient conjugu\'{e}: \label{algo gc}

   soit  $x \in \IR^n$ donn\'{e}.

   \begin{itemize}

     \item $  g^{0} = A x - b $

     \item $  h^{0} = - g^{0}$

     \item  ~pour~ $i=0$ ~\`{a}~ $ n$

    \begin{itemize}

     \item[-] $ \displaystyle \rho = -{(g^{i},h^i) \over (h,A h)}$

     \item[-] $ x^{i+1} = x^{i} + \rho h^{i}$

     \item[-] $ g^{i+1} = g^{i} + \rho A h^{i} $

     \item[-] $ \displaystyle \gamma = { ( g^{i+1},g^{i+1} ) \over (  g^{i},g^{i} )}$

     \item[-]  $h^{i+1} =  -g^{i+1} + gamma* h^{i}$

     \item[] si $( g^{i+1},g^{i+1} ) < \varepsilon $ stop

    \end{itemize}



  \end{itemize}
 \end{Algo}


Voila comment \'ecrire un gradient conjugue avec ces classes.

\subsection{Gradient conjugu\'e pr\'econditionn\'e}
\begin{listg}{GC.hpp}
\inputFF{GC.FF}
\end{listg}

\subsection{Teste du gradient conjugu\'e}
Pour finir voil\`a, un petit programme pour le test\'e sur cas
diff\'erent. Le troisi\`eme cas \'etant la r\'esolution de l'\'equation 
au diff\'erencielle 1d $ -u'' = 1 $ sur $[0,1]$  avec comme conditions
aux limites $u(0)=u(1)=0$, par la m\'ethode de l'\'el\'ement fini. La solution
exact est $ f(x) = x(1-x)/2$, nous verifions donc l'erreur sur le resultat.

\begin{listg}{GradConjugue.cpp}
\inputFF{GradConjugue.FF}
\end{listg}
\subsection{Sortie du teste}
\bFF
10x10   :
          10  -1  -1  -1  -1  -1  -1  -1  -1  -1
          -1  11  -1  -1  -1  -1  -1  -1  -1  -1
          -1  -1  18  -1  -1  -1  -1  -1  -1  -1
          -1  -1  -1  37  -1  -1  -1  -1  -1  -1
          -1  -1  -1  -1  74  -1  -1  -1  -1  -1
          -1  -1  -1  -1  -1 135  -1  -1  -1  -1
          -1  -1  -1  -1  -1  -1 226  -1  -1  -1
          -1  -1  -1  -1  -1  -1  -1 353  -1  -1
          -1  -1  -1  -1  -1  -1  -1  -1 522  -1
          -1  -1  -1  -1  -1  -1  -1  -1  -1 739
        GradienConjugue preconditionne par la diagonale 
6  ro = 0.990712 ||g||^2 = 1.4253e-24
 solution : A*x= 10     :
1.60635e-15       1       2       3       4       5       6       7       8       9

GradienConjugue preconditionnee par la identity  
9  ro = 0.0889083 ||g||^2 = 2.28121e-15
10      :
6.50655e-11       1       2       3       4       5       6       7       8       9

 GMRES 5 4.11771e-08 0
 solution : A*x= 10     :
-3.60463e-10      1       2       3       4       5       6       7       8       9

GradienConjugue laplacien 1D  par la identity  
48  ro = 0.00505051 ||g||^2 = 1.55006e-32
 Temps cpu = 0s
Fin err=5.55112e-17
\eFF


\begin{thebibliography}{77}
    
 
\bibitem{Buzzi}
G. Buzzi-Ferrari:  Scientific C++ Addison-Wesley 1993

\bibitem{Koenig}
A. Koenig (ed.) : Draft Proposed International Standard for Information Systems
- Programming Language C++. ATT report X3J16/95-087 (ark@research.att.com)

\bibitem{Pironneau}
B. Lucquin, O. Pironneau:  Introduction au calcul scientifique. Masson (1996).

\bibitem{Shapiro}
J. Shapiro:  A C++ Toolkit, Prentice Hall 1991.

\bibitem{Stroustrup}
B. Stroustrup: The C++ Programming Language Addison Wesley 1991

\end{thebibliography}

\end{document}