\LargeCaract�ristiques arithm�tiques des processeurs graphiques
Marc Daumas
Guillaume Da Gra�a et David Defour% \HAL{
[12pt] DALI-LP2A (Universit\'e de Perpignan)
52 avenue Paul Alduy --- 66860 Perpignan --- France
[6pt] LIRMM (CNRS, Universit� de Montpellier 2)
161 rue Ada --- 34392 Montpellier Cedex 5 --- France
#1
daumas_dagraca_defour_sympa2006.tex
% RECHERCHER LES TODO

\newcommand{\HAL}[1]{#1}
\newcommand{\NON}[1]{}


%====================================
% Title                 : Example d'article pour RenPar'2000
% Last Update           : 6 octobre 2001
% Modification          : 
%===========================================================
%%%
%  JM  Bruel nous a  permis  de r\'eutiliser  le style d\'evelopp\'e pour
%   ETR'97. On a change des trucs.
%%%% Si vous avez besoin d'ajouter des fichiers de
%%%% style faites le ou apparaissent les trois
%%%% points (...) et pas ailleurs. Pour ce qui de Latex 2e ca marche
%%%% pas mal aussi.
%%%%
% gillemets a la francaise
\def\ggf{\leavevmode\hbox{\raise.4ex\hbox{%
	$\scriptscriptstyle<\!<\;$}}}
\def\gdf{\leavevmode\hbox{\raise.4ex\hbox{%
	$\scriptscriptstyle\;>\!>\;$}}}
\def\leftnote#1{\leavevmode\vadjust{\setbox1=\vtop{\hsize 20mm
	\parindent=0pt\small\baselineskip=9pt
	\rightskip=4mm plus 4mm#1}
	\hbox{\kern-2cm\smash{\box1}}}}
% encore quelques petits symboles particuliers
	\font\myl=manfnt

	\def\panneau{{\myl\char"7F}}
	\def\boxone{{\myl\char"1C}}
	\def\boxtwo{{\myl\char"1D}}
	\def\ortf{{\myl\char"1E}}
	\def\fleurone{{\myl\char"26}}
	\def\fleurtwo{{\myl\char"27}}
	\def\diams{{\myl\char"23}}
	\def\cible{{\myl\char"24}}
	\def\carre{{\myl\char"25}}
	\def\fleche{{\myl\char"79}}
	\def\panneaubis{{\myl\char"7E}}
\def\panneauinverse{{\myl\char"00}}

% Conserver ces commandes (Debut)

\documentclass[twoside\HAL{,a4paper}]{article}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage[french]{babel}
\usepackage{listings,color}
\NON{\usepackage{perpi2006}}
\HAL{\usepackage{a4wide}}
\usepackage{euler,palatino}
\lstset{numbers=left,
  numberstyle=\tiny,
  stepnumber=1,
  numbersep=5pt,
  lineskip=-1pt,
  extendedchars=true,
  basicstyle=\footnotesize,
  breaklines,
  showstringspaces=false,
  frame=single,
  language={C},
  firstnumber=last,
  escapeinside={(*@}{@*)}}
% Conserver ces commandes (Fin)

\usepackage{epsfig,shadow,url}

\usepackage{amsmath}
\usepackage{amssymb}
\def\figdir{\HAL{.}\NON{../Fig}}
\def\figext{\HAL{pstex_t}\NON{pdftex_t}}

%===========================================================
%                               Title
%===========================================================
\newcommand{\filet}{\noindent\rule[0mm]{\textwidth}{0.2mm}}
\pagestyle{plain}

%\toappear{1} % Conserver cette ligne pour la version finale

\begin{document}

\parindent=0pt

\title{\Large\bf Caract�ristiques arithm�tiques des processeurs graphiques}
\NON{\shorttitle{Caract�ristiques arithm�tiques des processeurs graphiques}}

\author{Marc Daumas, Guillaume Da Gra�a et David Defour%
\HAL{
  \\[12pt]
  DALI-LP2A (Universit\'e de Perpignan)\\
  52 avenue Paul Alduy --- 66860 Perpignan --- France\\
  \{marc.daumas, guillaume.dagrac, david.defour\}@univ-perp.fr\\[6pt]
  LIRMM  (CNRS, Universit� de Montpellier 2)\\
  161 rue Ada --- 34392 Montpellier Cedex 5 --- France\\
  marc.daumas@lirmm.fr
}}%

\NON{
  \address{%
    DALI-LP2A (Universit\'e de Perpignan)\\
    52 avenue Paul Alduy --- 66860 Perpignan --- France\\
    \{marc.daumas, guillaume.dagrac, david.defour\}@univ-perp.fr\\[6pt]
    LIRMM  (CNRS, Universit� de Montpellier 2)\\
    161 rue Ada --- 34392 Montpellier Cedex 5 --- France\\
    marc.daumas@lirmm.fr
  }
}

\HAL{\date{}}

\maketitle

%===========================================================         %
%R\'esum\'e
%===========================================================  
\HAL{\newcommand{\Resume}[1]{\begin{abstract}#1\end{abstract}}}
\Resume{%
  Les unit�s graphiques ({\em Graphic Processing Units} --- GPU) sont
  d�sormais des processeurs puissants et flexibles. Les derni�res
  g�n�rations de GPU contiennent des unit�s programmables de traitement
  des sommets ({\em vertex shader}) et des pixels ({\em pixel shader})
  supportant des op�rations en virgule flottante sur 8, 16 ou 32 bits.
  La repr�sentation flottante sur 32 bits correspond � la simple pr�cision de la
  norme IEEE sur l'arithm�tique en virgule flottante (IEEE-754). Les
  GPU sont bien adapt�s aux applications avec un fort parall�lisme de
  donn�es.  Cependant ils ne sont que peu utilis�s en dehors des
  calculs graphiques ({\em General Purpose computation on GPU} ---
  GPGPU).  Une des raisons de cet �tat de faits est la pauvret� des
  documentations techniques fournies par les fabricants (ATI et
  Nvidia), particuli�rement en ce qui concerne l'implantation des
  diff�rents op�rateurs arithm�tiques embarqu�s dans les diff�rentes
  unit�s de traitement. Or ces informations sont essentielles pour
  estimer et contr�ler les erreurs d'arrondi ou pour mettre en
  {\oe}uvre des techniques de r�duction ou de compensation afin de
  travailler en pr�cision double, quadruple ou arbitrairement �tendue.
  Nous proposons dans cet article un ensemble de programmes qui
  permettent de d�couvrir les caract�ristiques  principales des GPU en ce qui concerne l'arithm�tique � virgule flottante.
  Nous donnons les r�sultats obtenus sur deux cartes graphiques
  r�centes: la Nvidia 7800GTX et l'ATI RX1800XL.
%
}%

\NON{\MotsCles{GPGPU, GPU, arithm�tiques virgule flottante, pr�cision.}}
	% 5 mots cl\'es seulement!!!

%=========================================================
\section{Introduction et pr�sentation des unit�s graphiques (GPU)}
%=========================================================

Les unit�s graphiques (GPU) sont des unit�s sp�cialis�es dans le
calcul intensif et r�gulier qui d�veloppent une puissance de calcul
bien sup�rieure � celle disponible sur les processeurs g�n�ralistes
\cite{Pharr2005}.  Avec l'arriv�e des GPU de derni�re g�n�ration, il devient 
int�ressant d'utiliser ces processeurs GPU pour faire des calculs
g�n�ralistes (GPGPU) \cite{Man05}~\footnote{Voir aussi le site
  \url{http://www.gpgpu.org/}.}.  Les GPU deviennent alors des
processeurs sp�cialis�s pour des applications r�guli�res et � fort parall�lisme de
donn�es \cite{LegRob03}.

% %=========================================================
% \subsection{Plan de ce document}
% %=========================================================

Apr�s une pr�sentation des unit�s graphiques (section \ref{sec:pipeline_graphique}) et des diff�rentes
implantations de l'arithm�tique � virgule flottante (section \ref{sec:virgule_flottante}), cette introduction
se termine par un �tat de l'art des travaux li�s aux tests des propri�t�s flottantes d'un syst�me (section \ref{sec:etat_art}). Nous pr�sentons ensuite nos algorithmes, leurs implantations et les
r�sultats obtenus (section \ref{sec:resultat}).

% Cet article se termine par des premi�res remarques
% et de nouvelles questions issues de nos tests exp�rimentaux.

%=========================================================
\subsection{Mod�le de fonctionnement des cartes graphiques: le pipeline graphique}
%=========================================================
\label{sec:pipeline_graphique}

\begin{figure}
  \begin{center}
    \input \NON{../Fig/}graphics_pipeline.\figext
  \end{center}
  \caption{Vue d'ensemble du pipeline graphique}
  \label{fig:graphics_pipeline}
\end{figure}

Les GPU traitent principalement des objets g�om�triques et des pixels.
Les images sont cr��es en appliquant des transformations g�om�triques
aux sommets et en d�coupant les objets en fragments ou pixels.  Les
calculs sont r�alis�s par diff�rents �tages de ce que l'on appelle le
pipeline graphique, comme pr�sent� � la
figure~\ref{fig:graphics_pipeline}. La machine h�te envoie des sommets
pour positionner dans l'espace des objets g�om�triques primitifs (polygones, lignes,
points). Ces objets primitifs subissent des transformations (rotations, translation, illumination\ldots) avant d'�tre assembl�s pour cr�er un objet plus complexe. Ces op�rations sont r�alis�es dans l'unit� de traitement des sommets ({\em vertex shader}).

Quand un objet a atteint sa position, sa forme et son �clairage
finals, il est d�coup� en fragments ou pixels. Une interpolation est
effectu�e pour obtenir les propri�t�s de chaque pixel. Les pixels sont
ensuite trait�s par l'unit� de traitement des pixels ({\em pixel
  shader}) qui n'effectue pas des transformations g�om�triques mais des
t�ches d'affichage comme par exemple appliquer une texture ou calculer
la couleur d'un pixel.

Le pipeline effectivement implant� diff�re l�g�rement de celui
repr�sent� dans la figure~\ref{fig:graphics_pipeline}. Selon les
cartes et les circuits, les constructeurs d�placent, partagent,
dupliquent ou ajoutent certaines ressources. La figure montre les
diff�rentes �tapes sur l'exemple d'un triangle.  Les vertex shaders traitent 3 sommets alors que les pixel shaders  traitent 17 pixels. Pour une figure donn�e, le
nombre de pixels est presque toujours plus important que le nombre de
sommets et les architectures modernes contiennent plus d'unit�s de
traitement des pixels que d'unit�s de traitement des sommets. Le ratio
actuel est par exemple de 24 pour 8.

% � l'aide, sur les processeurs de derni�re g�n�ration comme le Nvidia
% 7800GTX, d'une unit� MAD, d'une unit� de calcul des fonctions
% sp�ciales (sinus, cosinus, exponentielle, logarithme). Ces unit�s
% ont �galement acc�s � une m�moire appel�e texture. Le sch�ma d'un
% des vertex processor de la Nvidia 7800GTX est donn�e dans le
% graphique XXXX.


% , avec l'aide, sur les processeurs de derni�re g�n�ration comme le
% Nvidia 7800GTX, de 2 MAD cascad�s. Il a �galement acc�s au cache de
% texture. Le sch�ma d'un des vertex processor de la Nvidia 7800GTX
% est donn�e dans le graphique XXXX.

\begin{figure}
  \begin{center}
    \input \NON{../Fig/}vertexshader.\figext
  \end{center}
  \caption{D�tails d'un vertex shader pr�sent dans la Nvidia 7800GTX}
  \label{fig:vertexshader}
\end{figure}

La figure \ref{fig:vertexshader} reprise de \cite{Pharr2005} pr�sente
l'un des 8 vertex shaders de la Nvidia 7800 GTX. Dans la
classification de Flynn \cite{Fly66,Fly72}, les 8 shaders fonctionnement en
mode MIMD ({\em Multiple Instruction Multiple Data}) entre eux.  Chaque vertex
shader est capable d'initier � chaque cycle une op�ration {\em Multiply
  and Accumulate} (MAD) sur 4 tripl�s dans l'unit� vectorielle et une
op�ration {\em special} dans l'unit� scalaire.  Les op�rations  {\em special}
implant�es sont les fonctions exponentielles (exp, log), trigonom�triques (sin, cos) et deux fonctions inverses ($1/x$ et $1/\sqrt{x}$). Avec l'arriv�e de la
version 3 du support mat�riel Direct~3D, les shaders sont capables d'acc�der � la m�moire de texture
gr�ce � une unit� d�di�e.

\begin{figure}
  \begin{center}
    \input \NON{../Fig/}pixelshader.\figext
  \end{center}
  \caption{D�tails d'un pixel shader pr�sent dans la Nvidia 7800GTX}
  \label{fig:pixelshader}
\end{figure}

La figure~\ref{fig:pixelshader} d�crit l'un des 24 pixel shaders de la
Nvidia 7800 GTX.  Les 24 shaders fonctionnent en mode SIMD ({\em Single
  Instruction Multiple Data}) entre eux. La premi�re unit� flottante ex�cute
selon le programme 4 MAD ou un acc�s � la texture via
l'unit� de traitement des textures.  Le r�sultat est envoy� � la
deuxi�me unit� flottante qui ex�cute 4 MAD.  Dans le cas
de la Nvidia 7800 GTX, chaque pixel shader dispose d'une m�moire d�di�e appel�e cache
de texture de niveau 1  et d'une unit� capable de premiers traitements
sur les textures.

%=========================================================
\subsection{Implantation sur GPU de l'arithm�tique � virgule flottante en regard des normes ANSI et ISO}
%=========================================================
\label{sec:virgule_flottante}
Les GPU travaillent sur diff�rents formats de repr�sentation des
nombres. L'arithm�tique � virgule flottante offre un spectre de
nombres repr�sentables plus �tendu que l'arithm�tique enti�re ou �
virgule fixe sans que le programmeur ait besoin de g�rer manuellement
d�calages et recadrages. L'implantation de l'arithm�tique � virgule
flottante \cite{Dau05} est normalis�e par deux normes am�ricaines en cours de
r�vision\footnote{Voir le site \url{http://grouper.ieee.org/groups/754/}.}
IEEE-ANSI 754 \cite{Ste.81,Ste.87} et 854 \cite{CodKar.84} et par une
norme internationale IEC-ISO 60559 \cite{IEC89}. Des implantations
tr�s diff�rentes de l'arithm�tique � virgule flottante existent
\cite{BolDau04a} mais elles ne sont pas utilis�es par les GPU.

\begin{table}
  \caption{Format de repr�sentation des nombres � virgule flottante sur GPU et CPU}
  \label{fpformat}
  \begin{center}
    \begin{tabular}{|l|c|c|c|c|c|}\hline
      R�f�rence  & \multicolumn{4}{|c|}{Nombre de bits}                              & Valeurs non   \\ \cline{2-5}
                 & \parbox{2cm}{\centering Total} & \parbox{2cm}{\centering Signe} & \parbox{2cm}{\centering Exposant} & \parbox{2cm}{\centering Fraction} & num�riques    \\ \hline \hline
      Nvidia     & 16 & 1 & 5  & 10 & NaN, Inf        \\ \cline{2-5}
                 & 32 & 1 & 8  & 23 &                 \\ \hline
      ATI        & 16 & 1 & 5  & 10 & Absentes        \\ \cline{2-5}
                 & 24 & 1 & 7  & 16 &                 \\ \cline{2-6}
                 & 32 & 1 & 8  & 23 & Non document�es \\ \hline \hline
      ANSI-ISO   & 32 & 1 & 8  & 23 & NaN, Inf        \\ \cline{2-5}
                 & 64 & 1 & 11 & 52 &                 \\ \hline
    \end{tabular}
  \end{center}
\end{table}

Les normes pr�c�dentes repr�sentent un nombre par trois champs, une
mantisse $m$, un exposant $e$ et un signe $s$, la base de l'exposant
est fix�e � 2 ou � 10. Les GPU utilisent tous la base 2. La mantisse
est un nombre en virgule fixe avec 1 bit avant la virgule, l'exposant
est un entier biais�. Les normes pr�c�dentes imposent que le bit de
la mantisse � gauche de la virgule soit �gal � 1 ce qui est toujours
possible sauf si le nombre � repr�senter est trop petit. On parle selon les cas de repr�sentation normalis�e ou d�normalis�e.
Pour �conomiser de la m�moire, on ne stocke pas ce premier bit de la mantisse �gal � 1 pour les nombres normalis�s et on appelle
fraction $f$ les bits de la mantisse � droite de la virgule.  Ainsi,
un nombre flottant normalis� $x$, s'�crit
$$
x ~ = ~ (-1)^s ~ \cdot ~ m ~ \cdot ~ 2^e ~ = ~ (-1)^s ~ \cdot ~ 1,f ~ \cdot ~ 2^e.
$$

Les normes d�finissent deux formats pr�sent�s dans le
tableau~\ref{fpformat}. Le format simple pr�cision qui est cod� sur 32
bits et le format double pr�cision qui est cod� sur 64 bits. Pour
traiter les d�passements de capacit� ou les situations
exceptionnelles (division par z�ro), les normes d�finissent les valeurs non num�riques que
sont les infinis ou les non-nombres ({\em Not a Number} --- NaN).

Le r�sultat d'une op�ration arithm�tique ($+$, $\times$, $/$,
$\sqrt{~}$) n'est presque jamais repr�sentable bien que les op�randes
sont repr�sentables. Il faut l'arrondir.  Les normes imposent que
l'on utilise un des quatre modes d'arrondi pr�d�finis: par d�faut, par
exc�s, tronqu�, ou au plus pr�s. Des travaux anciens \cite{Coo78}
montrent que l'on peut arrondir pr�cis�ment les quatre op�rations
usuelles pr�c�dentes, en utilisant uniquement trois bits
suppl�mentaires par rapport au format du r�sultat. C'est-�-dire comme si les op�rations interm�diaires
avaient �t� infiniment pr�cises. On parle des  bits {\em guard}, {\em
  round} et {\em sticky}.

Une r�gle d'arrondi au plus proche ne d�finit pas un r�sultat unique
quand le r�sultat interm�diaire exact se situe � �gales distances
entre deux nombres repr�sentables.  Pour des raisons statistiques, l'arrondi normalis� favorise le
nombre dont la mantisse se termine par un 0 et on parle d'arrondi pair
\cite{ReiKnu75}.

Le tableau~\ref{fpformat} reprend les formats implant�s sur les cartes
graphiques que nous avons �tudi�es. Les donn�es concernant le support pour les valeurs non num�riques sont issues de \cite{Cebenoyan2005}. Le format sur 32 bits est li� � la
version 3.0 de Direct~3D.  Ce format s'inspire du format simple
pr�cision des norme ANSI-IEEE et ISO-IEC, mais les constructeurs ne garantissent
pas la compatibilit� de leur unit� avec ces normes. La documentation disponible sur les cartes graphiques ne permet pas de d�terminer en quoi ces implantations diff�rent de la norme. Il est donc n�cessaire de mieux
comprendre le fonctionnement de l'arithm�tique flottante sur GPU pour
ma�triser le comportement num�rique des applications scientifiques ex�cut�es sur GPU.


%=========================================================
\subsection{�tat de l'art et travaux ant�rieurs sur CPU et sur GPU}
%=========================================================
\label{sec:etat_art}
% %=========================================================
% \subsubsection{Sur les processeurs g�n�ralistes}
% %=========================================================

Les premiers logiciels r�alis�s pour tester l'arithm�tique � virgule
flottante sur CPU avaient pour unique but de d�couvrir les
caract�ristiques des fonctionnalit�s implant�es
\cite{Cod88,GenMar74,Sch81}. Suite � l'adoption  massive de la norme
IEEE-754, certains logiciels sont apparus pour v�rifier la bonne implantation de cette
norme par une s�rie de tests bien choisis
\cite{Kar85,Par99,VerCuyVer01a} alors que d'autres proposaient
d'implanter et de tester des fonctions �l�mentaires, sp�ciales ou
complexes \cite{Cod93a,Cod93b}.
%qui ne seront normalis�es que plus tard par les normes dites LIA \cite{IEC94,IEC95,IEC01}. 
Certains logiciels tels UCBTest\footnote{Ce programme est disponible sur le d�p�t Netlib
  \url{http://www.netlib.org}.} testent � la fois la conformit� des
fonctionnalit�s normalis�es, la qualit� des autres fonctionnalit�s
usuelles et le bon fonctionnement de la cha�ne de  compilation.

% Avec l'arriv�e de la norme IEEE-754 en 1985, Karpinski  a
% adapt� un logiciel de William Kahan �crit dans les ann�es 1980, pour
% explorer les caract�ristiques flottantes d'un syst�me. Ce syst�me
% teste l'addition, la multiplication, la division. Il est bas� sur la
% g�n�ration de cas difficiles � arrondir pour la multiplication, la
% division et la racine carr�e. De nombreux param�tres sont test�s tel
% que la pr�sence de bit de garde, des modes d'arrondi ou des
% d�normalis�s. 


% Un autre ensemble de programme permet de tester certains cas
% difficiles de la norme IEEE-754. Il s'agit de UCBTEST. Cet ensemble
% est bas� sur une partie de paranoia compl�t� par d'autre programme
% testant certaines caract�ristiques particuli�res du syst�me tel que la
% pr�sence du bug de la division du Pentium.I

% %=========================================================
% \subsubsection{Sur les GPU}
% %=========================================================

Comme nous l'avons �voqu�  � la section pr�c�dente, les GPU ne respectent pas la norme IEEE-754 contrairement � la
majorit� des CPU.  Certaines  hypoth�ses sur le comportement arithm�tique des op�rateurs flottants doivent �tre v�rifi�es avant
de pouvoir porter un programme pour CPU sur un GPU comme, par exemple, �muler des op�rateurs arithm�tiques en pr�cision multiple \cite{Graca2006}.

% Un m�me calcul peut donc produire des r�sultats diff�rents sur deux
% machines diff�rentes (un GPU fabriqu� par Nvidia et un par ATI). Dans
% ces conditions,

{\em Paranoia} \cite{Kar85} est un outil qui teste certaines
propri�t�s de l'arithm�tique flottante des CPU. Un sous-ensemble de
ces tests a �t� adapt� pour pouvoir �tre ex�cut� sur certains GPU
\cite{Hillesland2004}. \NON{Nous publions en annexe de notre rapport
  de recherche\footnote{Bient�t disponible sur HAL-CCSD
    \url{http://}.} les r�sultats complets du test de la carte
  graphique Nvidia 7800GTX.} \HAL{Les r�sultats complets du test de la
  carte graphique Nvidia 7800GTX sont publi�s en annexe.} En plus des
binaires du logiciel, les auteurs fournissent un tableau de synth�se
sur l'ATI R300 \cite{Hillesland2004}. � ce jour, ce programme ne
fonctionne pas sur notre ATI RX1800XL.
% Il permet de construire le
% tableau~\ref{paranoia}.
%
% \begin{table}
%   \begin{center}
%     \begin{tabular}{|l|c|c|c|c|} \hline
%       Operation&			Exact &	Chopped&	R300&	NV35 \\ 
%                &			rounding &	&	&	 \\ \hline \hline
%       Addition&				[-0.5, 0.5]&	(-1, 0]&	[-1.0, 0.0]&	[-1.0, 0.0]\\ \hline
%       Substraction&		[-0.5, 0.5]&  (-1, 1)&	[-1.0, 1.0]&	[-0.75, 0.75]\\ \hline
%       Multiplication&	[-0.5, 0.5]&	(-1, 0]&	[-0.989, 0.125]&	[-0.782, 0.625]\\ \hline
%       Division&				[-0.5, 0.5]&	(-1, 0]&	[-2.869, 0.094]&	[-1.199, 1.375]\\ \hline
%     \end{tabular}
%     \caption{R�sultat de l'ex�cution de Paranoia sur diff�rentes
%       cartes graphiques. Les r�sultats obtenus correspondent �
%       l'erreur flottante exprim� en ulp}
%     \label{paranoia}
%   \end{center}
% \end{table}
Ses sources ne sont pas disponibles ce qui nous
emp�che d'aller plus loin que les seuls messages g�n�r�s par
son ex�cution.  Par exemple, nous ne pouvons pas d�terminer
quels op�rateurs ont �t� test�s parmi tous les op�rateurs disponibles
dans les vertex et pixel shaders et quels vecteurs de test ont �t�
utilis�s au del� de {\em Paranoia}.

Certaines caract�ristiques se d�gagent de ces premiers tests:
\begin{itemize}
\item L'addition et la multiplication sont tronqu�es sur les deux GPU.
\item La soustraction semble b�n�ficier d'un bit de garde sur Nvidia
  mais pas sur ATI.
\item La multiplication implante un arrondi fid�le sur les deux GPU \cite{Pri92}.
\item L'erreur de division semble indiquer que la division est
  implant�e par la multiplication du dividende par une approximation
  de l'inverse du diviseur \cite{BriMulRai04}.
\end{itemize}

En conclusion, les travaux initi�s montrent que de nombreuses
questions restent sans r�ponse et que les r�ponses ne pourront �tre
apport�es que par un logiciel dont les sources sont disponibles. C'est
le cas pour les logiciels cit�s pr�c�demment et ciblant les CPU et pour les implantations des algorithmes d�crits dans
la suite qui sont disponibles sur demande aupr�s des auteurs pour �valuation.

% \begin{center}
% \url{http://gpgpu.univ-perp.fr}
% \end{center}

%=========================================================
\section{D�couverte des caract�ristiques des unit�s arithm�tiques des GPU}
%=========================================================
\label{sec:resultat}
Nous avons d�fini et implant� des algorithmes pour mieux comprendre le
fonctionnement de l'addition, de la multiplication et le stockage des
nombres � virgule flottante dans les registres et en m�moire. Bien que
ces algorithmes peuvent �tre adapt�s � d'autres formats de donn�es, nous
avons cibl� le format simple pr�cision sur 32 bits. Les algorithmes
ont �t� �crits dans une version pr�liminaire en OpenGL et utilisent les
{\em Frame Buffer Object} pour le stockage dans les textures. Ces
premiers programmes fonctionnent parfaitement sur la carte Nvidia 7800
GTX avec le driver ForceWare 81.98.  En revanche ils ne fonctionnent
pas sur la carte ATI RX1800XL et le driver Catalyst 6.3 et nous avons
r��crit les programmes en DirectX.  Nos algorithmes sont ainsi
disponibles pour OpenGL et DirectX.


%Il est n�cessaire de bien comprendre le fonctionnement de ces
%op�rateurs pour pouvoir d�velopper des algorithmes num�riques et des
%op�rateurs � pr�cision multiple \cite{Graca2006}. D'apr�s les
%r�sultats pr�liminaires \cite{Hillesland2004}, l'additionneur de la
%Nvidia comporte un bit de garde. Nous avons donc d�velopp� un algorithme
%qui nous permet de d�terminer le nombre exact de bits de garde de
%l'additionneur. Cet algorithme nous donne �galement la fa�on dont est
%r�alis� la troncature de l'addition.


% Nous avons pour cela
% cr�� des algorithmes nous permettant de d�duire des propri�t�s au fur
% et � mesure de nos r�sultats.

%=========================================================
\subsection{Stockage des nombres � virgule flottante dans les registres et en m�moire}
%=========================================================

% \paragraph{Valeurs sp�ciales}

Le GPU est souvent utilis� comme coprocesseur. Des donn�es sont
g�n�r�es dans le CPU pour �tre transf�r�es dans la m�moire de texture avant d'�tre trait�es par le GPU.
Pour savoir si des conversions interviennent lors des transferts du
CPU vers le GPU, nous avons envoy� des nombres d�normalis�s, des NaN
et des nombres infinis du CPU vers le GPU pour ensuite les r�cup�rer
sur le CPU. Nous avons aussi effectu� des op�rations sur GPU g�n�rant
ces valeurs sp�ciales et r�cup�r� le r�sultat sur CPU.

Les r�sultats obtenus montrent que lors d'un transfert sans aucune
op�ration, les nombres d�normalis�s sont remplac�s par 0 et les
infinis ne sont pas modifi�s sur Nvidia ou ATI. Les NaN ne sont pas
modifi�s sur Nvidia mais ATI transforme les sNaN ({\em signaling NaN})
en qNaN ({\em quiet NaN}).

Sur certaines architectures, les registres internes stockent les nombres avec une
pr�cision sup�rieure � celle des donn�es en m�moire \cite{Int97a} et avec une
dynamique plus grande pour l'exposant \cite{Mar2K}. La conversion vers le format final est r�alis�e lors de l'�criture en m�moire. Nous avons test� si certaines unit�s
des GPU se comportent de mani�re similaire. Pour tester l'exposant
maximum, nous avons utilis� des vecteurs pour calculer
$$(\text{MAX\_FLOAT} + \text{MAX\_FLOAT}) - \text{MAX\_FLOAT}$$
o� $\text{MAX\_FLOAT}$ est le plus grand nombre repr�sentable au format
simple pr�cision. Nous avons v�rifi� le nombre de bits de la mantisse en utilisant des vecteurs pour calculer
$$(1,5 - 2^{-i}) + 2^{-i}$$
pour $i$ variant de 1 � 64.



Nous avons voulu savoir si les diff�rents MAD, o� $\text{MAD}(x, y,
z)$ effectue l'op�ration $x \times y + z$, conservent dans l'accumulateur plus de bits de
mantisse du produit $x \times y$ que la pr�cision de travail pour les
utiliser quand les chiffres de poids forts du produit sont compens�s
par le troisi�me op�rande $z$. C'est le cas des architectures
normalis�es r�centes qui implantent un FMA en arithm�tique � virgule
flottante \cite{Mar2K,Fre05}.  Nous avons g�n�r� des vecteurs de tests
al�atoires et compar� le reste de la multiplication sur GPU avec la
valeur exacte calcul�e sur CPU,
$$\text{MAD} (x, y, - x \otimes y)$$
o� $x \otimes y$ est la multiplication arrondie au plus pr�s au format
simple pr�cision sur CPU.

% vecteur $Z[i] = X[i] * Y[i]$.
% et calcul� $z = x * y$ sur pour d�terminer la pr�cision de l'accumulateur.

% \begin{lstlisting}[caption={Test de l'accumulateur du MAD},firstnumber=1,label=algo_mad]
% // Cr�ation des donn�es sur le CPU
% for(i=0; i<N; i++) {
% 	dataX[i] = random_float();
% 	dataY[i] = random_float();
% 	dataZ[i] = - dataX[i]*dataY[i]
% }

% // Algorithme ex�cut� dans le shader
% return (dataX * dataY) + dataZ
% \end{lstlisting}

% \paragraph{Pr�cision des registres temporaires}


L'ex�cution des tests montre qu'aucun des registres temporaires des
vertex et pixel shaders des cartes ATI et Nvidia n'utilise une plage
d'exposant �tendue pour �viter les d�bordements ou un nombre de bits de
mantisse sup�rieur pour augmenter la pr�cision des calculs. Par
ailleurs, aucun des MAD ne conserve le produit sur une pr�cision
�tendue au del� des 24 bits de pr�cision des nombres flottants simple
pr�cision.

%=========================================================
\subsection{La multiplication}
%=========================================================

Les r�sultats de l'ex�cution de {\em Paranoia} nous laissent penser que
la multiplication est tronqu�e � la fois sur ATI et sur Nvidia. Nous
avons test� la fa�on dont cette troncature �tait effectu�e �
l'int�rieur des pixel shaders. Pour cela nous avons souhait�
d�terminer si tous les produits partiels �tait g�n�r�s ou si pour une
question de rapidit�, les bits de poids faibles de certains produits
partiels �taient ignor�s comme pr�sent� � la
figure~\ref{fig:produits_partiels} \cite{SchSwa93b}. Dans ce cas, une constante est
ajout�e � la somme des produits partiels pour r�duire  le biais
statistique introduit par la troncature d'une quantit� toujours
positive.

\begin{figure}
  \begin{center}
    \input \NON{../Fig/}Mul_1.\figext
    % \rotatebox{270}{\resizebox{!}{6cm}{\includegraphics{Mul_1.pdf}}}
  \end{center}
  \caption{Produits partiels tronqu�s pour acc�l�rer le multiplieur et r�duire sa taille}
  \label{fig:produits_partiels}
\end{figure}

Nous avons multipli� deux vecteurs de 24 nombres: un premier vecteur compos� uniquement de
la valeur $2^{23}+3$ et un deuxi�me vecteur d�fini par la r�currence suivante pour $i > 0$
$$
    y_1      = 1                  ~~~~~~~~~~ 
    y_{2i}   = y_{2i - 1} + 1     ~~~~~~~~~~
    y_{2i+1} = 4 y_{2i - 1} + 2.
$$
Cette suite peut aussi s'interpr�ter � l'aide des op�rateurs binaires \emph{ou logique} ($|$) et 
\emph{d�calage vers la gauche} ($<<$) en :
$$
    y_1      = 1                            ~~~~~~~~~~ 
    y_{2i}   = y_{2i - 1}        ~ | ~ 01_2  ~~~~~~~~~~
    y_{2i+1} = (y_{2i - 1} << 2) ~ | ~ 10_2.
$$
Ce deuxi�me vecteur alterne une valeur qui cr�e une cha�ne de 1 sans
propagation dans les bits de poids faibles pour v�rifier la
troncature, avec une valeur g�n�rant une propagation de retenue �
partir de la colonne la plus � droite jusqu'au premier bit du r�sultat
attendu.

Les r�sultats des tests nous indiquent que le multiplieur des pixel
shaders traite les produits partiels  jusqu'� la 9�me colonne
pour ceux d'ATI et jusqu'� la 6�me colonne pour Nvidia. Ces r�sultats
sont confirm�s par des tests al�atoires qui montrent qu'un
biais est ajout� pour compenser les produits partiels manquants.

Un deuxi�me test compos� de 2 vecteurs de $2^{23}$ nombres
nous donne plus de pr�cision. Il calcule les produits
$$(2^{23} + 1) \times (2i + 1) ~~~~~ \text{pour} ~~~~~0 < i < 2^{23}.$$
Ce test d�termine la valeur de la constante de
r�duction du biais statistique. 

% Sur ATI, la constante semble �tre de $1$ et sur Nvidia elle vaut $41$.

% TODO j'ai remplac� 2^-12 abh�rant avec un op�rateur qui ne garde que 9 bits par 2^-9...


%----------------------------------------------------------
% \paragraph{V�rification de l'�galit� d'arrondi lors de changement de signe}

% Avant de s'attaquer plus en profondeur dans la compr�hension du
% fonctionnement del'arrondi pour l'op�rateur de multiplication, il a
% fallu d'abord s'int�resser � l'�galit� des r�sultats obtenus lors d'un
% simple changement de signe effectu� sur les param�tres. On a donc
% chercher � savoir si, pour deux nombres A et B quelconques, les
% r�sultats des op�rations A*B, (-A)*(-B) sont �gaux, ainsi que (-A)*B
% et A*(-B) et enfin A*B et -(A*(-B)), comme le montre l'algorithme.

Nous avons test� sur deux  milliards de valeurs al�atoires sur ATI et
Nvidia que $A \times B$, $(-A) \times (-B)$ et $-(A \times (-B))$
sont �gaux, ainsi que $(-A) \times B$ et $A \times (-B)$. Selon toute
probabilit�, la notation signe-valeur absolue est utilis�e et le
multiplieur calcule de fa�on s�par�e le signe du r�sultat et sa valeur
absolue.





%=========================================================
\subsection{L'addition}
%=========================================================

%-------------------------------------------------------------
%\paragraph{Recherche du nombre de bit(s) de garde}

Lors de l'utilisation de l'outil {\em Paranoia}, nous avons remarqu�
que la soustraction sur les processeurs Nvidia semblait disposer d'un
bit de garde. En l'absence d'informations suppl�mentaires, nous avons
cherch� � clarifier le comportement de la soustraction des pixel et
vertex shaders de l'ATI et de la Nvidia. Les pixel shaders des deux
architectures test�es comportent deux MAD cascad�s. Nous avons d�termin�
le comportement du premier additionneur en calculant
$$
1,5 - 2^{-i}
$$
pour $i$ variant de 1 � 64. Nous avons remarqu� que le r�sultat obtenu
�tait �gal � $1,5$ d�s que $i \geq 26$ sur ATI et Nvidia dans les
vertex et pixel shaders. Ce r�sultat laisse penser que la premi�re
soustraction b�n�ficie de deux bits de garde.

De fa�on similaire nous avons test� le deuxi�me additionneur des pixel
shaders en calculant les deux soustractions suivantes
$$
(1,5 - 2^{-i})-1.5
$$
pour $i$ variant de 1 � 64. Le r�sultat retourn� �tait �gal � $0$
lorsque $i \geq 25$ sur ATI et $i \geq 26$ sur Nvidia. Les deux
additionneurs cascad�s des pixel shaders de la Nvidia se comportent donc
de fa�on similaire. En revanche, nous pouvons penser que dans le cas
o� une seule addition est lanc�e dans les pixel shaders d'ATI, alors
elle est effectu�e par le deuxi�me additionneur qui dispose de deux bits de
garde alors que lorsque deux additions sont lanc�es, les deux additionneurs
sont utilis�s et le premier ne dispose que d'un seul bit de garde.

L'utilisation d'un bit de garde permet de ne pas commettre d'erreur
d'arrondi quand on calcule la diff�rence de deux nombres proches mais
dont les exposants diff�rent d'une unit�.  En revanche, la pr�sence
d'un ou de plusieurs bits de garde suppl�mentaires dans une arithm�tique tronqu�e peut
faire que la propri�t� suivante n'est plus v�rifi�e \cite{Pic76,PicVig93}
$$e(x) - e(y) > t ~~~ \Longrightarrow ~~~ x \oplus y = x$$
% Th�or�me 6.1
o� $x \oplus y$ est le r�sultat de l'addition sur GPU, $e(\cdot)$ est
la valeur de l'exposant de la variable consid�r�e et $t$ est le nombre
de bits de mantisse du format de travail. Cette propri�t� est n�cessaire pour que la diff�rence
$$x + y - x \oplus y$$
est toujours repr�sentable en machine sauf d�passement de capacit�. Il faut par la suite modifier   certains algorithmes de calcul � pr�cision multiple 
\cite{Pri92,Graca2006}.
Nous continuons nos travaux pour comprendre
pourquoi les additionneurs des GPU test�s disposent de deux bits de
garde et non pas d'un seul.

% Pour mieux comprendre ces r�sultats, n

Nous avons lanc� les m�mes
additions que pr�c�demment mais en utilisant des flottants sur 16 bits
au lieu de 32. Les r�sultats obtenus sur Nvidia montrent que les
additions sont effectu�es � pr�cision maximale (26 bits de pr�cision)
pour ensuite �tre arrondies dans le format de destination souhait�
(11 bits).



%=========================================================
\section{Conclusion et perspectives}
%=========================================================

Nous avons vu des algorithmes adapt�s aux caract�ristiques des
op�rateurs flottants des processeurs graphiques.  Gr�ce � ces
algorithmes, nous avons pu tester les propri�t�s des additionneurs et
multiplieurs des vertex et pixel shaders de la Nvidia 7800 GTX ainsi
que de l'ATI RX1800XL. Nous avons, entre autre, montr� que les
registres temporaires stockent les nombres flottants uniquement sur
32 bits. Nous avons �galement d�termin� le comportement des
multiplieurs qui utilisent un biais pour compenser la troncature. Et
enfin nous avons montr� que les additionneurs disposent de 2 bits de
garde ce qui explique certaines surprises lors du calcul
de l'erreur d'arrondi. Ces r�sultats constituent une
premi�re �tape pour la construction d'algorithmes num�riques plus
pr�cis et plus rapides sur GPU. Cependant, de nombreux tests
compl�mentaires restent � �crire pour caract�riser pr�cis�ment le
comportement de ces op�rateurs et de ceux que nous n'avons pas encore
test�s (fonctions trigonom�triques, logarithmiques et inverses).
 
%===========================================================
% \section*{Bibliographie}
%===========================================================

\bibliographystyle{plain}
%\bibliography{alias,perso,groupe,saao,these,livre,arith,alternate}
\bibliography{alternate,daumas_ref}

\HAL{
  \appendix
  
  \section{Messages g�n�r�s par Paranoia pour la Nvidia 7800GTX}
  
  {\small
    \input resultat_nvidia_paranoia.txt
  }
}

\end{document}


daumas_dagraca_defour_sympa2006.bbl
\begin{thebibliography}{10}

\bibitem{IEC89}
Binary floating-point arithmetic for microprocessor systems.
\newblock ISO / IEC Standard 60559, International Electrotechnical Commission,
  1989.

\bibitem{BolDau04a}
Sylvie Boldo and Marc Daumas.
\newblock Properties of two's complement floating point notations.
\newblock {\em International Journal on Software Tools for Technology
  Transfer}, 5(2-3):237--246, 2004.

\bibitem{BriMulRai04}
Nicolas Brisebarre, Jean-Michel Muller, and Saurabh~Kumar Raina.
\newblock Accelerating correctly rounded floating point division when the
  divisor is known in advance.
\newblock {\em IEEE Transactions on Computers}, 53(8):1069--1072, 2004.

\bibitem{Cebenoyan2005}
Cem Cebenoyan.
\newblock Floating point specials on the {GPU}.
\newblock Technical report, Nvidia, february 2005.

\bibitem{Cod88}
William~J. Cody.
\newblock {MACHAR}: a subroutine to dynamically determine machine parameters.
\newblock {\em ACM Transactions on Mathematical Software}, 14(4):303--311,
  1988.

\bibitem{Cod93a}
William~J. Cody.
\newblock Algorithm 714: {CELEFUNT}: a portable test package for complex
  elementary functions.
\newblock {\em ACM Transactions on Mathematical Software}, 19(1):1--21, 1993.

\bibitem{Cod93b}
William~J. Cody.
\newblock Algorithm 715: {SPECFUN} -- a portable {F}ortran package of special
  function routines and test drivers.
\newblock {\em ACM Transactions on Mathematical Software}, 19(1):22--30, 1993.

\bibitem{CodKar.84}
William~J. Cody, Richard Karpinski, et~al.
\newblock A proposed radix and word-length independent standard for floating
  point arithmetic.
\newblock {\em IEEE Micro}, 4(4):86--100, 1984.

\bibitem{Coo78}
Jerome~T. Coonen.
\newblock Specification for a proposed standard for floating point arithmetic.
\newblock Memorandum ERL M78/72, University of California, Berkeley, 1978.

\bibitem{Dau05}
Marc Daumas.
\newblock {\em Informatique r�partie}, chapter Implantations de la norme IEEE
  754 de l'arithm�tique � virgule flottante.
\newblock Herm�s, 2005.

\bibitem{DauMat94}
Marc Daumas and David~W. Matula.
\newblock Rounding of floating point intervals.
\newblock {\em Interval Computations}, (4):28--45, 1994.

\bibitem{Fly66}
Michael Flynn.
\newblock Very high-speed computing systems.
\newblock {\em Proceedings of the {IEEE}}, 54(12):1901--1909, 1966.

\bibitem{Fly72}
Michael Flynn.
\newblock Some computer organizations and their effectiveness.
\newblock {\em IEEE Transactions on Computers}, 21(9):948--960, 1972.

\bibitem{Fre05}
Freescale Semiconductor.
\newblock {\em DSP56800: Family Manual}, 2005.

\bibitem{GenMar74}
W.~Morven Gentleman and Scott~B. Marovitch.
\newblock More on algorithms that reveal properties of floating point
  arithmetic units.
\newblock {\em Communications of the ACM}, 17(5), 1974.

\bibitem{Graca2006}
Guillaume~Da Gra�a and David Defour.
\newblock Implementation of float-float operators on graphics hardware.
\newblock Technical Report ccsd-00021443, HAL-CCSD, 2006.

\bibitem{Hillesland2004}
Karl Hillesland and Anselmo Lastra.
\newblock {GPU} floating-point paranoia.
\newblock In {\em ACM Workshop on General Purpose Computing on Graphics
  Processors}, page~C8, August 2004.

\bibitem{Int97a}
Intel.
\newblock {\em Pentium {II} Processor: Developer's Manual}, 1997.

\bibitem{Kar85}
Richard Karpinski.
\newblock {PARANOIA}: a floating-point benchmark.
\newblock {\em Byte}, 10(2):223--235, 1985.

\bibitem{LegRob03}
Arnaud Legrand and Yves Robert.
\newblock {\em Algorithmique parall�le}.
\newblock Dunod, 2003.

\bibitem{Man05}
Dinesh Manocha.
\newblock General purpose computations using graphics processors.
\newblock {\em IEEE Computer}, 38(8):85--88, 2005.

\bibitem{Mar2K}
Peter Markstein.
\newblock {\em {IA}-64 and elementary functions: speed and precision}.
\newblock Prentice Hall, 2000.

\bibitem{Par99}
Michael Parks.
\newblock Number theoretic test generation for directed rounding.
\newblock In Israel Koren and Peter Kornerup, editors, {\em Proceedings of the
  14th Symposium on Computer Arithmetic}, pages 241--248, Adelaide, Australia,
  1999.

\bibitem{Pharr2005}
Matt Pharr, editor.
\newblock {\em GPUGems 2 : Programming Techniques for High-Performance Graphics
  and General-Purpose Computation}.
\newblock Addison-Wesley, 2005.

\bibitem{Pic76}
Mich\`ele Pichat.
\newblock {\em Contributions \`a l'etude des erreurs d'arrondi en
  arithm\'etique \`a virgule flottante}.
\newblock PhD thesis, Universit\'e Scientifique et M\'edicale de Grenoble,
  Grenoble, France, 1976.

\bibitem{PicVig93}
Mich\`ele Pichat and Jean Vignes.
\newblock {\em Ing\'enierie du contr\^ole de la pr\'ecision des calculs sur
  ordinateur}.
\newblock Editions Technip, 1993.

\bibitem{Pri92}
Douglas~M. Priest.
\newblock {\em On properties of floating point arithmetics: numerical stability
  and the cost of accurate computations}.
\newblock PhD thesis, University of California at Berkeley, Berkeley,
  California, 1992.

\bibitem{ReiKnu75}
John~F. Reiser and Donald~E. Knuth.
\newblock Evading the drift in floating point addition.
\newblock {\em Information Processing Letters}, 3(3):84--87, 1975.

\bibitem{Sch81}
N.~L. Schryer.
\newblock A test of computer's floating-point arithmetic unit.
\newblock Technical report~89, AT\&T Bell Laboratories, 1981.

\bibitem{SchSwa93b}
Michael~J. Schulte and Earle~E. Swartzlander.
\newblock Truncated multiplication with correction constant.
\newblock In {\em Proceedings of the 6th IEEE Workshop on VLSI Signal
  Processing}, pages 388--396. IEEE Computer Society Press, 1993.

\bibitem{Ste.81}
David Stevenson et~al.
\newblock A proposed standard for binary floating point arithmetic.
\newblock {\em IEEE Computer}, 14(3):51--62, 1981.

\bibitem{Ste.87}
David Stevenson et~al.
\newblock An {A}merican national standard: {IEEE} standard for binary floating
  point arithmetic.
\newblock {\em ACM SIGPLAN Notices}, 22(2):9--25, 1987.

\bibitem{VerCuyVer01a}
Brigitte Verdonk, Annie Cuyt, and Dennis Verschaeren.
\newblock A precision and range independent tool for testing floating-point
  arithmetic {I}: basic operations, square root and remainder.
\newblock {\em ACM Transactions on Mathematical Software}, 27(1):92--118, 2001.

\end{thebibliography}
