Esope source files

   1 : C DDOT2     SOURCE    PV        10/12/09    21:15:01     6811           
   2 :       function ddot2(n,x,y)
   3 : *
   4 : *   produit scalaire compensé de deux vecteurs
   5 : *   algorithme de Ogita,Rump & Oishi
   6 : *   http://perso.ens-lyon.fr/nicolas.louvet/TheseLouvet07.pdf
   7 : *
   8 :       implicit real*8 (a-h,o-z)
   9 :       real*8 x(n),y(n)
  10 : *  dans la constante, l'exposant est la moitie du nombre de bits de la mantisse
  11 : *  soit 26 pour le 64 bits IEE
  12 : *  a ajuster eventuellement
  13 :       parameter (c=2.D0**26+1.D0)
  14 : *
  15 :       if (n.le.0) then
  16 :        ddot2=0.d0
  17 :        return
  18 :       endif
  19 : *  code non compense
  20 : **    res=0.D0
  21 : **    do i=1,n
  22 : **     res=res+x(i)*y(i)
  23 : **    enddo
  24 : **    ddot2=res
  25 : **    return
  26 : *
  27 : * twoprod x(1),y(1)
  28 :       z=c*x(1)
  29 :       xx=z-(z-x(1))
  30 :       xy=x(1)-xx
  31 :       z=c*y(1)
  32 :       yx=z-(z-y(1))
  33 :       yy=y(1)-yx
  34 : *
  35 :       s1=x(1)*y(1)
  36 :       c1=xy*yy-(((s1-xx*yx)-xx*yy)-xy*yx)
  37 : * boucle
  38 :       do i=2,n
  39 : * twoprod x(i),y(i)
  40 :       z=c*x(i)
  41 :       xx=z-(z-x(i))
  42 :       xy=x(i)-xx
  43 :       z=c*y(i)
  44 :       yx=z-(z-y(i))
  45 :       yy=y(i)-yx
  46 :       s2=x(i)*y(i)
  47 :       c2=xy*yy-(((s2-xx*yx)-xx*yy)-xy*yx)
  48 : * twosum s1 s2
  49 :       xx=s1+s2
  50 :       z=xx-s1
  51 :       xy=(s1-(xx-z))+(s2-z)
  52 :       s1=xx
  53 :       c1=c1+(xy+c2)
  54 :       enddo
  55 :       ddot2=s1+c1
  56 : 
  57 :       end
  58 : 
  59 :