Fast and Accurate Computation ofSum and Dot ProductTakeshi Ogita (JST/Waseda Univ.)
joint work with
Siegfried M. Rump (TUHH) and Shin’ichi Oishi (Waseda Univ.)
Dagstuhl Seminar, Sep. 29, 2005
Paper available in SISC 26:6 (2005), 1955–1988.

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Contents
1. Outline
2. Error-free “sum” and “product”
3. Error-free vector summation
4. Accurate summation
5. Accurate dot product6. Numerical examples
7. Applications
8. Conclusion

Fast and Accurate Computation of Sum and Dot Product – 1

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Outline
The quality of a result obtained by ﬂoating point arithmetic on
computers sometimes happens to be poor because of the roundingerrors.
Purpose To compute summation
nPn
P

xi and dot product xT y =i=1xi yi accurately by ﬂoating point arithmetic.

i=1=⇒ one of the most basic task in scientiﬁc computing
Fast and Accurate Computation of Sum and Dot Product – 2

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Background
There are a number of algorithms to compute a summation or adot product. With respect to the result of algorithm we can grossly
distinguish three diﬀerent approaches:i) heuristic/numerical evidence for improvement of result accuracy,ii) result accuracy as if computed in higher precision,iii) result with prescribed accuracy.
Our methods belong to the second class.
Fast and Accurate Computation of Sum and Dot Product – 3

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Disadvantages of previous methods• sorting of input data is necessary, eitheri) by absolute value or,ii) by exponent,• the inner loop contains branches,• besides working precision, some extra (higher) precision isnecessary,• access to mantissa and/or exponent is necessary.
Fast and Accurate Computation of Sum and Dot Product – 4

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Our approach• Only using ordinary (one working precision) ﬂoating point
addition, subtraction and multiplication• Based on Knuth’s algorithm and Veltkamp-Dekker’s algorithm• No branch, no access to mantissa nor exponent=⇒ Extra higher precision or special hardware is NOT necessary=⇒ From these properties, our algorithms are very fast in terms of
not only ﬂops but also execution time
Fast and Accurate Computation of Sum and Dot Product – 5

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Error-free transformation
F : a set of t-bit (working precision) ﬂoating point numbers
u := 2−t: relative precision of ﬂoating point arithmetic
For example, t = 53 in IEEE 754 double precision format (u = 2−53)
For a, b ∈ F, we can obtain the following x, y ∈ F:¶

³Knuth (1969) [x, y] = TwoSum(a, b)⇒Dekker (1971) [x, y] = TwoProduct(a, b)µx+y =a+b⇒x+y = a×b

Error-free! (no loss of information)
Fast and Accurate Computation of Sum and Dot Product – 6´

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

TwoSum (Knuth, 1969)[x, y] = TwoSum(a, b)⇒a+b=x+y(|y| ≤ u|x|)⊕, ⊖ : ﬂoating point addition and subtraction¶³

function [x, y] = TwoSum(a, b)
x = a ⊕ b;
c = x ⊖ a;y = (a ⊖ (x ⊖ c)) ⊕ (b ⊖ c);µ´Computational cost: 6 ﬂops
Fast and Accurate Computation of Sum and Dot Product – 7Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

TwoSum’ (Dekker, 1971)[x, y] = TwoSum′(a, b)⇒a+b=x+y(|y| ≤ u|x|)¶³function [x, y] = TwoSum′(a, b)
x = a ⊕ b;
if |a| ≥ |b|
y = b ⊖ (x ⊖ a);else
y = a ⊖ (x ⊖ b);µ´

3 ﬂops + taking 2 absolute values + 1 comparison=⇒ Computational speed slows down
Fast and Accurate Computation of Sum and Dot Product – 8

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Table 1: Comparison of measured computing time for Knuth’s
TwoSum with Dekker’s TwoSum′n ratio (Dekker/Knuth)
100
1.342
1000
2.210
10000
3.046

CPU: Pentium 4 (3.46GHz), GNU Compiler Collection 3.4.2 (g77)=⇒ We prefer to use Knuth’s TwoSum.
Fast and Accurate Computation of Sum and Dot Product – 9

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Split (Dekker, 1971)[aH , aL] = Split(a)⇒a = aH + aL

⌊ 2t ⌋(|aL| ≤ 2|aH |)⊗ : ﬂoating point multiplication¶³function [aH , aL] = Split(a)⌈ 2t ⌉
c = factor ⊗ a; % factor = 2 + 1
aH = c ⊖ (c ⊖ a);aL = a ⊖ aH ;µ´

4 ﬂopsFast and Accurate Computation of Sum and Dot Product – 10

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

TwoProduct (Veltkamp [Dekker, 1971])[x, y] = TwoProduct(a, b)⇒a×b=x+y(|y| ≤ u|x|)¶³function [x, y] = TwoProduct(a, b)
x = a ⊗ b;
[aH , aL] = Split(a); % aH + aL ← a[bH , bL] = Split(b);% bH + bL ← by = aL ⊗ bL ⊖ (((x ⊖ aH ⊗ bH ) ⊖ aL ⊗ bH ) ⊖ aH ⊗ bL);µ´

17 ﬂopsFast and Accurate Computation of Sum and Dot Product – 11

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

TwoProductFMA
FMA: Fused Multiply-Add
By FMA(a, b, c), d := a ∗ b + c can be calculated in 1 ﬂop. Resultd is rounded to the nearest ﬂoating point number.
¶³function [x, y] = TwoProductFMA(a, b)
x = a ⊗ b;
y = FMA(a, b, −x);µ´

17 ﬂops =⇒ 2 ﬂops
Fast and Accurate Computation of Sum and Dot Product – 12

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Error-free vector summation (1)For p ∈ Fn, we can cascade TwoSum as follows:
p2
p1

TwoSum

q2

pn 1

p3
2

TwoSumq3

3

n 2

TwoSum

qn 1

pn
n 1

TwoSumnqn

Fast and Accurate Computation of Sum and Dot Product – 13

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

VecSump′ = VecSum(p)⇒nXi=1pi =n
X

p′i

i=1

(so-called distillation algorithm)¶³function p = VecSum(p)for i = 2 : n[pi, pi−1] = TwoSum(pi, pi−1);µ´

6(n − 1) ﬂops
Fast and Accurate Computation of Sum and Dot Product – 14

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Error-free vector summation (2)VecSum can iteratively be applied:nXnnX
XVecSum′ VecSum′′ VecSum
pi −−−−−−−→
pi −−−−−−−→
pi −−−−−−−→ · · ·i=1

i=1

=⇒nX
i=1

i=1

pi =n
X

p′i =nXp′′i = · · ·i=1

i=1

Error-free!
Fast and Accurate Computation of Sum and Dot Product – 15

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

p2

p3

p4

p5p′2

p′3

p′4

p′5p′′2p′′3p′′4

p1p′1p′′5p′′1p′′′5p′′′
1p′′′
2p′′′
3p′′′
4

Example: n = 5, K = 4.
Fast and Accurate Computation of Sum and Dot Product – 16

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

SumKres = SumK(p, K)⇒res ∼nX

pii=1¶³function res = SumK(p, K)
for k = 1 : K − 1
p = VecSum(p);endµn−1 ¶P
res = pn ⊕ ﬂpi ;µi=1´(6K − 5)n ﬂops
Fast and Accurate Computation of Sum and Dot Product – 17

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Error analysis for SumKn·uγn :=(u := 2−53 in IEEE 754 double precision)
1−n·uPX|pi|cond(pi) := P(condition number of summation)| pi |
Denote by res ∈ F the result obtained by SumK, thenP
X|res − pi|2KP≤ u + 3γn−1 + γ2n−2cond(pi).
| pi |P=⇒ Basically, relative error bound is u + cnu cond( pi)K=⇒ K-fold working precision arithmetic
Fast and Accurate Computation of Sum and Dot Product – 18

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Accurate dot productWe can extend the discussion to dot product using TwoProduct.
For x, y ∈ Fn, using [hi, ri] = TwoProduct(xi, yi)nXi=1xi yi =n
Xi=1

Ã(hi + ri)

=2n
X!
ti

i=1

Error-free! (provided no underﬂow occurs)

Fast and Accurate Computation of Sum and Dot Product – 19

Dagstuhl Seminar 2005

¶T. OGITA, S. M. RUMP and S. OISHI

Dot2

³

function res = Dot2(x, y)[p, s] = TwoProduct(x1, y1);for i = 2 : n[h, r] = TwoProduct(xi, yi);[p, q] = TwoSum(p, h);s = s ⊕ (q ⊕ r);
endres = p ⊕ s;µ´

25n ﬂops
Fast and Accurate Computation of Sum and Dot Product – 20

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Error analysis for Dot2
x, y ∈ Rn⇒

Condition number of xT y is deﬁned by
2|xT | |y|cond(x y) :=.
T|x y|
T

Denote by res ∈ F the result obtained by Dot2, then¯¯T¯ res − x y ¯¯ ≤ u + 1 γn2 cond (xT y).
¯¯ xT y ¯2=⇒ u + cnu2cond(xT y) (doubled working precision arithmetic)
Fast and Accurate Computation of Sum and Dot Product – 21

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

DotK¶

³

function res = DotK(x, y, K)[p, r1] = TwoProduct(x1, y1);for i = 2 : n[h, ri] = TwoProduct(xi, yi);[p, rn+i−1] = TwoSum(p, h);end
r2n = p;res = SumK(r, K − 1);µ´(12K + 1)n ﬂops
Fast and Accurate Computation of Sum and Dot Product – 22

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHIError analysis for DotK
Denote by res ∈ F the result obtained by DotK, then
1 K|res − xT y|2T≤
u+2γ+γcond(xy).
4n−24n−2T|x y|2=⇒ Basically, relative error bound is u + cnuK cond(xT y)=⇒ K-fold working precision arithmetic

Fast and Accurate Computation of Sum and Dot Product – 23

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Numerical example (1)Let us evaluate accuracy of Dot2 and DotK.
• n = 2000 (1000 samples)• Elements of x, y ∈ Fn: random numbers in [−1, 1] s.t. cond(xT y)
is from 1 to 10120|res − xT y|• Plot the relative error
where res is the result of Dot2
T|x y|or DotK (K = 3 : 7)
Fast and Accurate Computation of Sum and Dot Product – 24

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI
Algorithms Dot2 and DotK, K=3:7, n = 2000, 1000 samples
010−5

relative error

10−1010−1510−20

1020

10401060

10condition number80

10100

10

Fast and Accurate Computation of Sum and Dot Product – 25Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Numerical example (2)Let us evaluate speed of matrix-vector product using Dot2.
• n: 100, 200, . . . , 3000
• For A ∈ Fn×n, x, b ∈ Fn, y = Ax − b is calculated• Display the ratio of the measured computing time of ourroutine compared to that of the BLAS routine DGEMV overall dimensions

Fast and Accurate Computation of Sum and Dot Product – 26

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Table 2: Measured ratio of computing time for matrix-vector residual
env. BLAS Dot2 XBLAS (F) XBLAS (C)
∗)1
4.8
9.3
11.3
theor.
1
12.5
18.5
18.5
∗) Pentium IV 2.53 GHz, Fortran (CVF 6.6C)
C (MS VC++ 6.0), BLAS: Intel MKL 7.0

Fast and Accurate Computation of Sum and Dot Product – 27

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Application to linear systems
Let us consider Ax = b with A ∈ Rn×n, b ∈ Rn. Denote by x0 an
approximate solution obtained by LU decomposition ( 23 n3 ﬂops).
Then, Dot2 is suitable for iterative reﬁnement:
1. Calculate residual r̃ ∼ b − Axℓ using Dot2 (O(n2) ﬂops)
2. Solve Ay = r̃ using LU factors (O(n2) ﬂops)
3. Update xℓ+1 = xℓ + ỹ
In addition, this can be combined with self-validating method forlinear systems. We use INLTAB, a fast interval toolbox for MATLAB.
Fast and Accurate Computation of Sum and Dot Product – 28

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHI

Table 3: Results with/without residual iteration (n = 1000)cond(A)
105
1091013
LU
max. rel. err. x0
4.7e-12 2.8e-08 1.9e-04
max. rel. err. xℓ
1.8e-16 1.8e-16 1.8e-16
number of iterations ℓ335
ratio computing time
1.38
1.38
1.65
1Pentium IV 2.53 GHz, Matlab 6.5, INTLAB 4
Remark: For larger dimensions, the additional cost (including
Matlab’s interpretation overhead) relatively decreases.
Fast and Accurate Computation of Sum and Dot Product – 29

Dagstuhl Seminar 2005

T. OGITA, S. M. RUMP and S. OISHIConclusion• Develop and analyze accurate algorithms for sum and dot product• Fast in terms of not only ﬂops but also measured computing time
• Only using one working precision (portable and scalable)

Fast and Accurate Computation of Sum and Dot Product – 30