C DDOT2     SOURCE    PV        10/12/09    21:15:01     6811           
      function ddot2(n,x,y)
*
*   produit scalaire compensé de deux vecteurs
*   algorithme de Ogita,Rump & Oishi
*   http://perso.ens-lyon.fr/nicolas.louvet/TheseLouvet07.pdf
*
      implicit real*8 (a-h,o-z)
      real*8 x(n),y(n)
*  dans la constante, l'exposant est la moitie du nombre de bits de la mantisse
*  soit 26 pour le 64 bits IEE
*  a ajuster eventuellement
      parameter (c=2.D0**26+1.D0)
*
      if (n.le.0) then
       ddot2=0.d0
       return
      endif
*  code non compense
**    res=0.D0
**    do i=1,n
**     res=res+x(i)*y(i)
**    enddo
**    ddot2=res
**    return
*
* twoprod x(1),y(1)
      z=c*x(1)
      xx=z-(z-x(1))
      xy=x(1)-xx
      z=c*y(1)
      yx=z-(z-y(1))
      yy=y(1)-yx
*
      s1=x(1)*y(1)
      c1=xy*yy-(((s1-xx*yx)-xx*yy)-xy*yx)
* boucle
      do i=2,n
* twoprod x(i),y(i)
      z=c*x(i)
      xx=z-(z-x(i))
      xy=x(i)-xx
      z=c*y(i)
      yx=z-(z-y(i))
      yy=y(i)-yx
      s2=x(i)*y(i)
      c2=xy*yy-(((s2-xx*yx)-xx*yy)-xy*yx)
* twosum s1 s2
      xx=s1+s2
      z=xx-s1
      xy=(s1-(xx-z))+(s2-z)
      s1=xx
      c1=c1+(xy+c2)
      enddo
      ddot2=s1+c1

      end