Télécharger ddot2.eso

Retour à la liste

Numérotation des lignes :

  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.  
  60.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales