Accelerate Iterative Methods• Good Algorithms• Good Preconditioners• Parallel Algorithms• Good Implementations• Accurate Computations

Mixed Precision Iterative Methods
using High Precision Arithmetic
Hidehiko Hasegawa
hasegawa@slis.tsukuba.ac.jp
Faculty of Library, Information and Media Science,
University of Tsukuba

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

1

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

Convergence history䣄䣣䣵䣫䣥䢢䣯䣧䣶䣪䣱䣦䣵䢢䢯䢢䣆䣱䣷䣤䣮䣧䢯䣣䣴䣫䣶䣪䣯䣧䣶䣫䣥䢲䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯

Numerical Comparison of
Accelerating Polynomials in
Product-type Iterative Methods䢴䢯䢴䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢶

䢯䢸

䢯䢺

䢯䢳䢲

䢯䢳䢴

Seventh SIAM Conference on
Applied Linear Algebra 2000

H. Hasegawa, K. Abe, and S.-L. Zhang䢲

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢳䢲䢲

䢴䢲䢲

䢵䢲䢲
䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

䢶䢲䢲

䢷䢲䢲

䢸䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

2

Convergence history of Bi-CG part

Convergence of Bi-CG part: Quadruple( reconstruct Bi-CG using alpha and beta in each methods)( reconstruct Bi-CG using alpha and beta in each methods)䣴䣧䣥䣱䣰䣵䣶䣴䣷䣥䣶䣧䣦䢢䣄䣫䣅䣉䢢䣷䣵䣫䣰䣩䢢䣄䣫䣅䣉䢯䣲䣣䣴䣶䢢䣫䣰䢢䣧䣣䣥䣪䢢䣯䣧䣶䣪䣱䣦䣵䢢䣫䣰䢢䣓䣷䣣䣦䣴䣷䣲䣮䣧
䢴

䢲

䢲

䢯䢴

䢯䢴

䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯䣴䣧䣥䣱䣰䣵䣶䣴䣷䣥䣶䣧䣦䢢䣄䣫䣅䣉䢢䣷䣵䣫䣰䣩䢢䣄䣫䣅䣉䢯䣲䣣䣴䣶䢢䣫䣰䢢䣧䣣䣥䣪䢢䣯䣧䣶䣪䣱䣦䣵
䢴

䢯䢶

䢯䢸
䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢺䢯䢸䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢺

䢯䢳䢲

䢯䢳䢲

䢲

䢳䢲䢲

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢴䢲䢲

䢵䢲䢲
䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

䢶䢲䢲

䢷䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

How Bi-CG part works?
• Bi-CGSTAB converges by an effect of MR part( Bi-CG part is still unstable)• GPBi-CG makes Bi-CG part stable• CGS did not converge in Quadruple arithmetic• In Quadruple arithmetic, simple Bi-CG is the best
( Bi-CG is much affected by Rounding errors)• In Quadruple arithmetic, Bi-CG part in BiCGSTAB is bad convergence even if Bi-CG
converges.

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢯䢳䢴

䢸䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

䢲

䢳䢲䢲

䢴䢲䢲

䢵䢲䢲
䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢶䢲䢲

䢷䢲䢲

䢸䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

Convergence history based on one Bi-CG( alpha and beta in Bi-CG are used in all methods)䣣䣮䣲䣪䣣䢢䣣䣰䣦䢢䣤䣧䣶䣣䢢䣣䣴䣧䢢䣤䣣䣵䣧䣦䢢䣱䣰䢢䣶䣪䣣䣶䢢䣱䣨䢢䣄䣫䣅䣉䢢䢴䢲䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯

䢯䢳䢴

䢯䢶

䢯䢴

䢯䢶

䢯䢸
䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢺

䢯䢳䢲

䢯䢳䢴

䢲

䢳䢲䢲

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢴䢲䢲

䢵䢲䢲
䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

䢶䢲䢲

䢷䢲䢲

䢸䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

Convergence history based on one Bi-CG

Convergence history based on one Bi-CG(Quadruple arithmetic is used for Bi-CG)(Quadruple arithmetic is used for ALL)䣄䣣䣵䣫䣥䢢䣯䣧䣶䣪䣱䣦䣵䢢䢯䢢䣓䣷䣣䣦䣴䣷䣲䣮䣧䢢䣣䣴䣫䣶䣪䣯䣧䣶䣫䣥䣣䣮䣲䣪䣣䢢䣣䣰䣦䢢䣤䣧䣶䣣䢢䣤䣣䣵䣧䣦䢢䣱䣰䢢䣄䣫䣅䣉䢢䣫䣰䢢䣓䣷䣣䣦䣴䣷䣲䣮䣧䢯䣣䣴䣫䣶䣪䣯䣧䣶䣫䣥䢴䢴䢲䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯䣔䣧䣮䣣䣶䣫䣸䣧䢢䣴䣧䣵䣫䣦䣷䣣䣮䢢䣰䣱䣴䣯

䢲

䢯䢴

䢯䢶

䢯䢸

䢯䢺

䢯䢳䢴

䢯䢶

䢯䢳䢲

䢯䢳䢴
䢲

䢷䢲

䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢸䢯䢺䣄䣫䣅䣉䢢䢢䢢䢢
䣅䣉䣕䢢䢢䢢䢢䢢
䣄䣫䣅䣉䣕䣖䣃䣄
䣉䣒䣄䣫䣅䣉䢢䢢

䢯䢳䢲

䢯䢴

䢳䢲䢲

䢳䢷䢲

䢴䢲䢲

䢴䢷䢲

䢲

䢳䢲䢲

䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

Seventh SIAM Conference on
Applied Linear Algebra 2000

H. Hasegawa, K. Abe, and S.-L. Zhang

Seventh SIAM Conference on
Applied Linear Algebra 2000

䢴䢲䢲

䢵䢲䢲
䣋䣶䣧䣴䣣䣶䣫䣱䣰䣵

䢶䢲䢲

䢷䢲䢲

䢸䢲䢲

H. Hasegawa, K. Abe, and S.-L. Zhang

How accelerating polynomial works• Qudaruple arithmetic works very well.
• If enough accuracy was provided, Bi-CG was thebest.
• Bi-CGSTAB and GPBi-CG work well.
• In Quadruple arithmetic, sometimes it works asbraking not as accelerating.
• GPBi-CG is robust in both two conditions.
• CGS does not work in both conditions because of“squared”.

Seventh SIAM Conference on
Applied Linear Algebra 2000

H. Hasegawa, K. Abe, and S.-L. Zhang

Utilizing Quadruple-Precision
Floating Point Arithmetic
Operation for the Krylov
Subspace Methods

SIAM Conference on
Applied Linear Algebra 2003

12

BiCG

Gamma = 1.3

SIAM Conference on
Applied Linear Algebra 2003

CGS

13

SIAM Conference on
Applied Linear Algebra 2003

SIAM Conference on
Applied Linear Algebra 2003

BiCG

BiCGSTAB Gamma = 1.3

15

Gamma = 1.3

SIAM Conference on
Applied Linear Algebra 2003

14

Gamma = 2.5

16

BiCGSTAB Gamma = 2.5

SIAM Conference on
Applied Linear Algebra 2003

GPBiCG Gamma = 2.5

17Observations

18

High Precision Arithmetic• Fast and smooth convergence are gained from
More accurate computations.
• Required Mantissa is based on the problems:
BiCG 53 bit for Gamma = 1.3
100 bit for
1.7
200 bit for2.1
200 bit for2.5• Required Mantissa depends on Algorithms:
BiCG
200 bit and 190 iterations
CGS
300 bit and 160
x BiCGSTAB 1500 bit and 210
x GPBiCG 300 bit and 310(Gamma = 2.5 )SIAM Conference on
Applied Linear Algebra 2003

SIAM Conference on
Applied Linear Algebra 2003

19

• Reducing round-off errors• Accelerating algorithms mathematically• Not easy to use

JSIAM Applied Mathematics Seminar, Dec. 27, 201320

High Precision Arithmetic
without any Special Hardware

Important points!

• Symbolic Computation (Computer Algebra)• Variable length Multiple Precision– GMP– MBLAS– exflib• Fixed length Multiple Precision– FORTRAN REAL *16– IEEE– Double-double

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

• Full or Partial
One Precision or Mixed Precisions• Computing Environment
Compiler/Emulation/Interpreter• Program Interface, API

21

Our Solution:

Utilize Accurate Computations
for Iterative Methods• Use Double-double• Use D-D vectors and Double Matrices(Mixed Precision Arithmetic Operations)• Accelerate by SIMD
• Restart with different Precision• Automatic Tuning• Good tools
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

23

JSIAM Applied Mathematics Seminar, Dec. 27, 201322

Advantages• Tough for round-off errors• Small Additional Memory• Small Additional Communications• Much Computations• Applicable for ALL Iterative Methods(even if serial computation such as ILU)

JSIAM Applied Mathematics Seminar, Dec. 27, 201324㻰㼛㼡㼎㼘㼑㻙㻰㼛㼡㼎㼘㼑㻔㻰㻰㻕㻘㻌㻽㼡㼍㼐㻙㻰㼛㼡㼎㼘㼑㻔㻽㻰㻕㻻㼚㼑㻌㻰㻰㻌㼚㼡㼙㼎㼑㼞㻌㼡㼟㼑㼟㻌㼠㼣㼛㻌㼐㼛㼡㼎㼘㼑㻌㼜㼞㼑㼏㼕㼟㼕㼛㼚㻌㼚㼡㼙㼎㼑㼞㼟㻚㻻㼚㼑㻌㻽㻰㻌㼚㼡㼙㼎㼑㼞㻌㼡㼟㼑㼟㻌㼒㼛㼡㼞㻌㼐㼛㼡㼎㼘㼑㻌㼜㼞㼑㼏㼕㼟㼕㼛㼚㻌㼚㼡㼙㼎㼑㼞㼟㻚

Implementation of Fast
Quadruple Arithmetic Operations

㻽㻰㻌㻌Ͻulp (㼡㼚㼕㼠㼟㻌㼕㼚㻌㼠㼔㼑㻌㼘㼍㼟㼠㻌㼜㼘㼍㼏㼑)㻭㼞㼕㼠㼔㼙㼑㼠㼕㼏㻌㼛㼜㼑㼞㼍㼠㼕㼛㼚㻌㼕㼟㻌㼜㼑㼞㼒㼛㼞㼙㼑㼐㻌㼎㼥㻌㼡㼟㼕㼚㼓㻌㼚
㼚㼛㼞㼙㼍㼘 㼐㼛㼡㼎㼘㼑㻌㼜㼞㼑㼏㼕㼟㼕㼛㼚㻌㼛㼜㼑㼞㼍㼠㼕㼛㼚㼟㻚
JSIAM Applied Mathematics Seminar, Dec. 27, 201325

Implementation of Double-double Arithmetic• Quadruple value is stored in two double
floating point numbers

a + b = fl(a + b) + err(a + b)
– a,b䠖double floating point variables– fl(a + b) 䠖addition of a and b in double– err(a+b)䠖(a+b) – fl(a+b): error

double-double arithmetic
Mantissa 52bit+

Exponent
11bit

Round-off Error Free
Double Arithmetic Addition• Round-off error free addition can be donewith two double precision variables:– Double-double arithmetic: a = a.hi + a.lo,|a.hi|>|a.lo|– 8 bits less than IEEE standards– Effective digits are approx. 31 vs. 33 digits.
exponent
11 bit㻰㻚㻌㻴㻚㻌㻮㼍㼕㼘㼑㼥㻘㻌㻽㻰㻌㻔㻯㻗㻗㻌㻛㻌㻲㼛㼞㼠㼞㼍㼚㻙㻥㻜㻌㼐㼛㼡㼎㼘㼑㻙㼐㼛㼡㼎㼘㼑㻌㼍㼚㼐㻌㼝㼡㼍㼐㻙㼐㼛㼡㼎㼘㼑㻌㼜㼍㼏㼗㼍㼓㼑㻕㻘㻭㼢㼍㼕㼘㼍㼎㼘㼑㻌㼍㼠㻌㼔㼠㼠㼜㻦㻛㻛㼏㼞㼐㻚㼘㼎㼘㻚㼓㼛㼢㻛䁺㼐㼔㼎㼍㼕㼘㼑㼥㻛㼙㼜㼐㼕㼟㼠㻛
26

Mantissa 52 bit

IEEE Standards
exponent
15bit

Mantissa 112bit
JSIAM Applied Mathematics Seminar, Dec. 27, 201327

JSIAM Applied Mathematics Seminar, Dec. 27, 201328

Quadruple Addition of a=b+c

Basic Algorithmb.hi• Dekker showed round-off error free
addition in double
3flops.

FAST_TWO_SUM(a,b,s,e)s = a + b
e = b - (s - a)

a

s

b
a
s

Others 6flops.
TWO_SUM(a,b,s,e)s = a + b
v = s - ae = (a - (s - v)) + (b - v)

c.lo

s-a

+

c.hi

=

sh

c.lo

=

b.lo+c.loeh

B. Addition of lower parts
b.lo+

C. Addition of result and error of upper partb.lo+c.lo

b

+

eh

=

sh

e

sh
29

+

eh

=

eh

eh
sh

D. FAST_TWO_SUM of results A and C

s-a
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

c.hi

A. TWO_SUM for upper partsb.hi|a|>=|b|

+

b.lo

a.hi

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

eh

a.lo
30

Quadruple Addition of a=b+c
ADD(a,b,c)
20 flops
TWO_SUM(b.hi,c.hi,sh,eh)
TWO_SUM(b.lo,c.lo,sl,el)
eh = eh + sl
FAST_TWO_SUM(sh,eh,sh,eh)
eh = eh + el
FAST_TWO_SUM(sh,eh,a.hi,a.lo)

a=(a.hi,a.lo), b=(b.hi,b.lo), c=(c.hi,c.lo)

JSIAM Applied Mathematics Seminar, Dec. 27, 201331

JSIAM Applied Mathematics Seminar, Dec. 27, 201332

Minimal Requirement

MuPAT• +/• *• /• SQRT• Input function• Output function (print)• Others
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

䈄 㻹㼡㻼㻭㼀 㻔㻹㼡㼘㼠㼕㼜㼘㼑㻌㻼㼞㼑㼏㼕㼟㼕㼛㼚㻌㻭㼞㼕㼠㼔㼙㼑㼠㼕㼏㻌㼀㼛㼛㼘㼎㼛㼤㻕㻌㼇㻞㼉䇵 㻰㼛㼡㼎㼘㼑㻘㻌㻰㻰㻘㻌㼍㼚㼐㻌㻽㻰㻌㼍㼟㻌㻿㼏㼕㼘㼍㼎 㼠㼛㼛㼘㼎㼛㼤䇵 㻭㼏㼏㼑㼘㼑㼞㼍㼠㼕㼛㼚㻌㼡㼟㼕㼚㼓㻌㻯㻌㼑㼤㼠㼑㼞㼚㼍㼘㻌㼒㼡㼚㼏㼠㼕㼛㼚㼟㼀㻚㻌㻿㼍㼕㼠㼛㻘㻌㻱㻚㻌㻵㼟㼔㼕㼣㼍㼠㼍 㼍㼚㼐㻌㻴㻚㻌㻴㼍㼟㼑㼓㼍㼣㼍㻘㻌㻰㼑㼢㼑㼘㼛㼜㼙㼑㼚㼠㻌㼛㼒㻌㼝㼡㼍㼐㼞㼡㼜㼘㼑㻌㼜㼞㼑㼏㼕㼟㼕㼛㼚㼇㼍㼞㼕㼠㼔㼙㼑㼠㼕㼏㻌㼠㼛㼛㼘㼎㼛㼤㻌㻽㼡㻼㻭㼀 㼛㼚㻌㼟㼏㼕㼘㼍㼎㻘㻌㻵㻯㻯㻿㻭㻞㻜㻝㻜㻘㻌㻼㼞㼛㼏㼑㼑㼐㼕㼚㼓㼟㻌㻼㼍㼞㼠㻌㻵㻵㻘㻌㻔㻞㻜㻝㻜㻕㻿㻚㻌㻷㼕㼗㼗㼍㼣㼍㻘㻌㼀㻚㻌㻿㼍㼕㼠㼛㻘㻌㻱㻚㻌㻵㼟㼔㼕㼣㼍㼠㼍 㼍㼚㼐㻌㻴㻚㻌㻴㼍㼟㼑㼓㼍㼣㼍㻘㻌㻰㼑㼢㼑㼘㼛㼜㼙㼑㼚㼠㻌㼍㼚㼐㻌㼍㼏㼏㼑㼘㼑㼞㼍㼠㼕㼛㼚㻌㼛㼒㻌㼙㼡㼘㼠㼕㼜㼘㼑㼜㼞㼑㼏㼕㼟㼕㼛㼚㻌㼍㼞㼕㼠㼔㼙㼑㼠㼕㼏㻌㼠㼛㼛㼘㼎㼛㼤㻌㻹㼡㻼㻭㼀 㼒㼛㼞㻌㻿㼏㼕㼘㼍㼎 㻘㻌㻶㻿㻵㻭㻹㻌㻸㼑㼠㼠㼑㼞㼟㻘 㼂㼛㼘㻚㻌㻡㻘㻌㼜㼜㻚㻌㻥㻙㻝㻞㻌㻔㻞㻜㻝㻟㻕

33

34

㻿㼏㼍㼘㼍㼞㻘㻌㼂㼑㼏㼠㼛㼞㻘㻌㼍㼚㼐㻌㻹㼍㼠㼞㼕㼤㻌㼍㼞㼑㻌㼠㼞㼑㼍㼠㼑㼐㻌㼍㼟㻌䇾㼏㼛㼚㼟㼠㼍㼚㼠䇿㻌㼐㼍㼠㼍㻌㼠㼥㼜㼑㻌㼕㼚㻌㻿㼏㼕㼘㼍㼎

Ease of Use• By trial and error• No Programming• Interactive• Combination of D, DD, and QD• Any machine

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

JSIAM Applied Mathematics Seminar, Dec. 27, 2013㼏㼛㼚㼟㼠㼍㼚㼠
㼏㼛㼐㼑

35

㼟㼕㼦㼑

㻿㼏㼍㼘㼍㼞

㼍㻌㻩㻌㻝㻧

㻝㽢㻝1㼂㼑㼏㼠㼛㼞

㼍㻌㻩㻌㼇㻝㻧㻞㼉㻧

㻞㽢㻝

1
2

㻹㼍㼠㼞㼕㼤

㼍㻌㻩㻌㼇㻝㻘㻟㻧㻞㻘㻠㼉㻧

㻞㽢㻞

1 3
2 4

JSIAM Applied Mathematics Seminar, Dec. 27, 201336䇾㼟㼜㼍㼞㼟㼑䇿㻌㼐㼍㼠㼍㻌㼠㼥㼜㼑㻌㼕㼚㻌㻿㼏㼕㼘㼍㼎㻱㼤㼠㼑㼚㼟㼕㼛㼚㻌㼛㼒㻌㻰㼍㼠㼍㻌㼠㼥㼜㼑㼟㻌㻒㻌㻻㼜㼑㼞㼍㼠㼛㼞㼟㼏㼛㼚㼟㼠㼍㼚㼠

㻯㻻㻻㼒㼛㼞㼙㼍㼠㻔㻯㼛㼛㼞㼐㼕㼚㼍㼠㼑㻌㼘㼕㼟㼠㻕㼝㼐㻦㻌㼛㼏㼠㼡㼜㼘㼑㼐㼐㻦㻌㼝㼡㼍㼐㼞㼡㼜㼘㼑

2

5

1
4
3

㻞㽢㼏㼛㼚㼟㼠㼍㼚㼠㻠㽢㼏㼛㼚㼟㼠㼍㼚㼠 䇾㼟㼜㼍㼞㼟㼑䇿㻌㻌㼐㼍㼠㼍㻌㼠㼥㼜㼑38㻭㼐㼐㼕㼠㼕㼛㼚㻌㼍㼚㼐㻌㻿㼡㼎㼠㼞㼍㼏㼠㼕㼛㼚

㻰㻰
㻱㼤㼠㻚

㼐㼐

㼏㼛㼚㼟㼠㼍㼚㼠

㼏㼛㼚㼟㼠㼍㼚㼠

㻱㼤㼠㻚㼟㼜㼍㼞㼟㼑㼢㼍㼘㼡㼑

JSIAM Applied Mathematics Seminar, Dec. 27, 201337㻱㼤㼠㼑㼚㼟㼕㼛㼚㻌㼛㼒㻌㼐㼍㼠㼍㻌㼠㼥㼜㼑㼟㻌㼕㼚㻌㼠㼣㼛㻌㼐㼕㼞㼑㼏㼠㼕㼛㼚㼟㼐㼛㼡㼎㼘㼑㼞㼛㼣 㼏㼛㼘㼡㼙㼚¾ 㻻㼜㼑㼞㼍㼠㼛㼞㼟㻌 㻔㻌㻗㻘㻌㻙㻘㻌㻖㻘㻌㻛㻌㻕 㼍㼞㼑㻌㼜㼞㼛㼢㼕㼐㼑㼐㻚¾ 㻾㼑㼐㼡㼏㼕㼚㼓㻌㼙㼑㼙㼛㼞㼥㻌㼟㼜㼍㼏㼑㻦
9 㻏㻌㼛㼒㻌㻺㼛㼚㻌㼦㼑㼞㼛㻌㼕㼟㻌㻡䠂䊻㻣㻚㻡㻑㻌㼛㼒㻌㼏㼛㼚㼟㼠㼍㼚㼠㻌㼐㼍㼠㼍㻌㼠㼥㼜㼑Same operators {+, -, *, /} and functions should be defined
between these data types.
JSIAM Applied Mathematics Seminar, Dec. 27, 2013㼍㻌㻌㻩㻔㻌㻌㻌㻌㻠㻘㻌㻌㻌㻌㻠㻕㻌㼟㼜㼍㼞㼟㼑㻌㼙㼍㼠㼞㼕㼤
㻔㻌㻌㻌㻌㻝㻘㻌㻌㻌㻌㻞㻕㻌㻌㻌㻌㻌㻌㻌㻌㻞㻚㻌㻔㻌㻌㻌㻌㻝㻘㻌㻌㻌㻌㻠㻕㻌㻌㻌㻌㻌㻌㻌㻌㻡㻚㻌㻔㻌㻌㻌㻌㻞㻘㻌㻌㻌㻌㻝㻕㻌㻌㻌㻌㻌㻌㻌㻌㻝㻚㻌㻔㻌㻌㻌㻌㻟㻘㻌㻌㻌㻌㻟㻕㻌㻌㻌㻌㻌㻌㻌㻌㻠㻚㻌㻔㻌㻌㻌㻌㻠㻘㻌㻌㻌㻌㻞㻕㻌㻌㻌㻌㻌㻌㻌㻌㻟㻚㻌

JSIAM Applied Mathematics Seminar, Dec. 27, 2013㼐㼐㼟㼜㼟㼜㼍㼞㼟㼑 39㻲㼡㼚㼏㼠㼕㼛㼚㻌㼚㼍㼙㼑㻌%sp_a_qd%qd_a_sp%qdsp_a_qdsp%qdsp_a_ddsp%ddsp_a_qdsp%qdsp_a_sp%sp_a_qdsp%qdsp_a_qd%qd_a_qdsp%qdsp_a_dd%dd_a_qdsp%qdsp_a_s%s_a_qdsp%ddsp_a_qd%qd_a_ddsp㼛㼜㼑㼞㼍㼠㼕㼛㼚
sp + qd
qd + sp
qdsp + qdsp
qdsp + ddsp
ddsp + qdsp
qdsp + sp
sp + qdsp
qdsp + qd
qd + qdsp
qdsp + dd
dd + qdsp
qdsp + double
double + qdsp
ddsp + qd
qd + ddsp¾ 㻽㻰㻌㼟㼟㼜㼍㼞㼟㼑㻝㻜㻜㻜㻌㼠㼕㼙㼑㼟䠖㼝㼐㼟㼜㻘㻌㻺㻌㻩㻝㻜㻜㻜㻘㻌㻏㻌㼛㼒㻌㻺㼛㼚㻌㼦㼑㼞㼛㻌㻡㻑㻿㼏㼕㼘㼍㼎 䞉䞉䞉㻠㻣㻚㻢㻡㻌㼟㼑㼏㻯㻌㼒㼡㼚㼏㼠㼕㼛㼚㼟㻌䞉䞉䞉 㻥㻞㻚㻟㻡㻌㼟㼑㼏
Scilab only
¾ 㻿㼍㼙㼑㻌㼍㼘㼓㼛㼞㼕㼠㼔㼙㼟㻌㼣㼕㼠㼔㻌㼐㼑㼚㼟㼑

JSIAM Applied Mathematics Seminar, Dec. 27, 201340㻹㼡㼘㼠㼕㼜㼘㼕㼏㼍㼠㼕㼛㼚㻲㼡㼚㼏㼠㼕㼛㼚㻌㼚㼍㼙㼑%sp_m_qd%qd_m_sp%qdsp_m_qdsp%qdsp_m_ddsp%ddsp_m_qdsp%qdsp_m_sp%sp_m_qdsp%qdsp_m_qd%qd_m_qdsp%qdsp_m_dd%dd_m_qdsp%qdsp_m_s%s_m_qdsp%ddsp_m_qd%qd_m_ddsp㼛㼜㼑㼞㼍㼠㼕㼛㼚
sp * qd
qd * sp
qdsp * qdsp
qdsp * ddsp
ddsp * qdsp
qdsp * sp
sp * qdsp
qdsp * qd
qd * qdsp
qdsp * dd
dd * qdsp
qdsp * double
double * qdsp
ddsp * qd
qd * ddsp

Acceleration by C functions¾ 㻰㻰㻌㼟㼜㼍㼞㼟㼑㻌㻝㻘㻜㻜㻜㻌㼠㼕㼙㼑㼟䠖㼐㼐㼟㼜 㻺㻩㻝㻘㻜㻜㻜㻘㻌㻏㻌㼛㼒㻌㻺㼛㼚㻌㼦㼑㼞㼛㻌㼕㼟㻌㻡㻑䠖㼐㼐# of double
MuPAT# of double
DD

QD¾ 㻿㼜㼍㼞㼟㼑㻌㻖㻌㻿㼜㼍㼞㼟㼑㻌㻩㻌㻫㻌㼇㻠㼉㼇㻠㼉㻌㼀㼕㼙㼛㼠㼔㼥㻌㻭㻚㻌㻰㼍㼢㼕㼟㻘㻌㻰㼕㼞㼑㼏㼠㻌㻹㼑㼠㼔㼛㼐㼟㻌㼒㼛㼞㻌㻿㼜㼍㼞㼟㼑㻌㻸㼕㼚㼑㼍㼞㻌㻿㼥㼟㼠㼑㼙㼟㻘㻌㻿㻵㻭㻹㻘㻌㻼㼔㼕㼘㼍㼐㼑㼘㼜㼔㼕㼍㻌㻔㻞㻜㻜㻢㻕㻚䈜
41㻯㻌㼒㼡㼚㼏㼠㼕㼛㼚㼟㻿㼏㼕㼘㼍㼎

㻹㼡㻼㻭㼀㼏㼛㼚㼟㼠㼍㼚㼠㼐㼐㼝㼐㻗㻘 㻙㻘㻌㻖㻘㻌㻛㻗㻘 㻙㻘㻌㻖㻘㻌㻛㻗㻘 㻙㻘㻌㻖㻘㻌㻛㼒㼡㼚㼏㼠㼕㼛㼚㼟㼒㼡㼚㼏㼠㼕㼛㼚㼟㼒㼡㼚㼏㼠㼕㼛㼚㼟㼟㼜㼍㼞㼟㼑㼐㼐㼟㼜㼝㼐㼟㼜㻗㻘 㻙㻘㻌㻖㻘㻌㻛㻗㻘 㻙㻘㻌㻖㻘㻌㻛㼒㼡㼚㼏㼠㼕㼛㼚㼟㼒㼡㼚㼏㼠㼕㼛㼚㼟㻗㻘 㻙㻘㻌㻖㻘㻌㻛

1

1

0.016

0.014

0.013

11

24

27

MuPAT

0.21 (12.8)

0.39 (28.4)

0.39 (30.6)

MuPAT with C

0.26 (16.4)

0.31 (22.8)

0.32 (24.9)# of double

Use C functions㻻㼢㼑㼞㼢㼕㼑㼣㻌㼛㼒㻌㼚㼑㼣㻌㻹㼡㻼㻭㼀1d㻿㼏㼕㼘㼍㼎䞉䞉䞉㻝㻝㻟㻡㻚㻣㻤 㼟㼑㼏㻯㻌㼒㼡㼚㼏㼠㼕㼛㼚䞉䞉䞉 㻝㻜㻚㻥㻞 㼟㼑㼏

JSIAM Applied Mathematics Seminar, Dec. 27, 2013㻿㼑㼏㻔㻌㼞㼍㼠㼕㼛㻌㼢㼟㻌㼏㼛㼚㼟㼠㼍㼚㼠㻌㻕

Repeated 10^6 times

91

217

649

MuPAT

2.91 (181.7)

4.21 (309.7)

21.29 (1663.5)

MuPAT with C

0.34 (21.1)

0.39 (28.3)

0.39(30.3)

Intel Core i5 2.5GHz, 4GB, Windows 7,Scilab version 5.3.3
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

42

Memory Consumption㻱㼤㼠㻚㼒㼡㼚㼏㼠㼕㼛㼚㼟

JSIAM Applied Mathematics Seminar, Dec. 27, 201343

JSIAM Applied Mathematics Seminar, Dec. 27, 201344

Result of matrix operations(Memory)

Result of matrix operations(Computation Time, 100 times)

6%
11%
15%
40%
63%
99%

JSIAM Applied Mathematics Seminar, Dec. 27, 201345

JSIAM Applied Mathematics Seminar, Dec. 27, 201346

BiCG for ill-conditioned
Problems
Lis & Lis-test
a Library of Iterative Solvers forlinear systems

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

JSIAM Applied Mathematics Seminar, Dec. 27, 201347

48

Lis has more than 10*13*11 combinations㻼㼞㼑㼏㼛㼚㼐㻚

㻿㼛㼘㼢㼑㼞㼟

㻿㼠㼛㼞㼍㼓㼑㻌㻲㼛㼞㼙㼍㼠

Jacobi㻯㻳㻯㻾㻿㻦㻌㻯㼛㼙㼜㼞㼑㼟㼟㼑㼐㻌㻾㼛㼣

SSOR

㻮㼕㻯㻳

㻯㻯㻿㻦㻌㻯㼛㼙㼜㼞㼑㼟㼟㼑㼐㻌㻯㼛㼘㼡㼙㼚ILU(k)

㻯㻳㻿Hybrid㻮㼕㻯㻳㻿㼀㻭㻮

㻹㻿㻾㻦㻌㻹㼛㼐㼕㼒㼕㼑㼐㻌㻯㼛㼙㼜㼞㼑㼟㼟㼑㼐㻌㻿㼜㼍㼞㼟㼑㻌㻾㼛㼣

I+S㻮㼕㻯㻳㻿㼀㻭㻮㻔㼘㻕㻰㻵㻭㻦㻌㻰㼕㼍㼓㼛㼚㼍㼘

SAINV㻳㻼㻮㼕㻯㻳

㻱㻸㻸㻦㻌㻱㼘㼘㼜㼍㼏㼗㻙㻵㼠㼜㼍㼏㼗㻌㼓㼑㼚㻚㻌㼐㼕㼍㼓㻚

SA-AMG㻮㼕㻯㻳㻿㼍㼒㼑

㻶㻰㻿㻦㻌㻶㼍㼓㼓㼑㼐㻌㻰㼕㼍㼓㼛㼚㼍㼘

Crout ILU

㻻㼞㼠㼔㼛㼙㼕㼚㻔㼙㻕㻯㻻㻻㻦㻌㻯㼛㼛㼞㼐㼕㼚㼍㼠㼑

additiveschwarz㻳㻹㻾㻱㻿㻔㼙㻕

㻰㻺㻿㻦㻌㻰㼑㼚㼟㼑㼀㻲㻽㻹㻾㻮㻿㻾㻦㻌㻮㼘㼛㼏㼗㻌㻿㼜㼍㼞㼟㼑㻌㻾㼛㼣User
defined㻶㼍㼏㼛㼎㼕

㻮㻿㻯㻦㻌㻮㼘㼛㼏㼗㻌㻿㼜㼍㼞㼟㼑㻌㻯㼛㼘㼡㼙㼚㻳㼍㼡㼟㼟㻙㻿㼑㼕㼐㼑㼘㼂㻮㻾㻦㻌㼂㼍㼞㼕㼍㼎㼘㼑㻌㻮㼘㼛㼏㼗㻌㻾㼛㼣

Steps
1. Initialize
2. Make matrix
3. Make vector
4. Define Solver
5. Set Values
6. Set conditions
7. Execute
8. Finalize/,6B0$75,;$/,6B9(&725E[/,6B62/9(5VROYHULQWLWHUGRXEOHWLPHVLWLPHVSWLPHVOLVBLQLWLDOL]H DUJFDUJY OLVBPDWUL[BFUHDWH /,6B&200B:25/' $ OLVBYHFWRUBFUHDWH /,6B&200B:25/' E OLVBYHFWRUBFUHDWH /,6B&200B:25/' [ OLVBVROYHUBFUHDWH VROYHU OLVBLQSXW $E[DUJY>@ OLVBYHFWRUBVHWBDOO E OLVBVROYHUBVHWBRSWLRQ& VROYHU OLVBVROYH $E[VROYHU OLVBVROYHUBJHWBLWHUV VROYHU LWHU OLVBVROYHUBJHWBWLPHV VROYHU WLPHVLWLPHV SWLPHV SULQWI LWHU GWLPH H S H
L H gQLWHUWLPHVSWLPHVLWLPHV OLVBILQDOL]H 㻿㻻㻾
JSIAM Applied Mathematics Seminar, Dec. 27, 201349

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

Design of Fast Quad. Operations for Lis

Acceleration• Same API with Double• Double: Input 䠄A䠈b䠈䡔0 䠅• Double: Output• Double: Creation of Preconditioner M• Fast Quad.: Iterative solution x
All working variables• Fast Quad.: Application of Preconditioner
Mu=v• SIMD is used for vectors䠄dot, axpy, matvec䠅

JSIAM Applied Mathematics Seminar, Dec. 27, 201350

– SSE2: 2 Multiply-and-add in same time
– AVX: 4 Multiply-and-add in same time
– AVX2: 4 Fused Multiply-and-add in same time

• 2 or 4 FMA in a loop with loop unrolling– pd instruction of SSE2 can be used for all• Code tuning– Alignment– Some hand optimization

51

JSIAM Applied Mathematics Seminar, Dec. 27, 201352

Architecture of intel core i7 2600k

Vector operations (BLAS1)
Kogakuin University

Kogakuin University

Normal instruction (SISD)
SSE2
AVX
Floating pointadder

x1 x2 x3 x4+ + + +
y1 y2 y3 y4x1 x2 x3 x4
㽢 㽢 㽢 㽢y1 y2 y3 y4

256bit load 㽢2
256bit store 㽢1

Floating point
multiplier

256bit load 㽢2
256bit store 㽢1

SIMD register

Operation

Memory access(Load + Store)

The number ofdouble precision operations(add+sub : mult)axpy

y = ax + y

4+2

35 (26:9)

axpyz

z = ax + y

4+2

35 (26:9)

xpay

y = x + ay

4+2

35 (26:9)

dot

val = x 䞉 y

4+0

35 (26:9)

nrm2

val = ||x2||

2+0

31 (24:7)

scale

x = ax2+2

24 (15:9)

y Addition and Multiplication in parallel
x, y, z : DD vector; a, val : DD variabley Peak performance
3.4G 㽢 4 (AVX) 㽢 2 (adder + multiplier) = 27.2 GFLOPS / core
= 108.8 GFLOPS (4core)
10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

53

Performances of DD vector operations (in cache)
Kogakuin University“GFLOPS” := (# of double precision op. * N) / elapsed time
10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

Performances of multi-threading (axpy)

Kogakuin University

120

x3.4

x3.7
x3.4

x3.4

x3.4

40
56%

53%

52%

59%

54%

51%

27.2
20
62%

63%

62%

65%

65%

Peak
(1 core)

axpyz

xpay

dot

nrm2

40
30

Bounded bymemory bandwidth (21.2GB/s)

20

109%
82%

76%

76%

0

scale

2

Name of vector operations (N = 105)
10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

50

10

62%

0
axpy

Performance [GFLOPS]

Performance [GFLOPS]

x3.7

1 thread2 threads
3 threads
4 threads

L3 cache size

60

1 thread4 threads

80
60

70

Peak
(4 core)

108.8
100

54

2.6

4

6

8

Size of vector N (x105)55

10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 201356

Reducing memory access

Performance of SpMV

Kogakuin University

50• CRS (compressed row storage) format is used
Performance [GFLOPS]

45

Bytes / flops of y = Ax
Bytes / flops
yD = ADxD

14 (28 bytes / 2 flops)

yDD = ADDxDD

1.5 (52 bytes / 35 flops)

䖂 yDD = ADxDD

1.3 (44 bytes / 33 flops)• yDD= ADxDD

M

S

25
M

M

L

S

15
10

45.0

Performance [GFLOPS]

L3 cache size

30

x3.2

20
76%

58

Kogakuin University

238%

Not bounded bymemory bandwidth (21.2GB/s)

Name of matrix (nnz/row)

Performance of Transposed SpMV
50.0

40
Performance[GFLOPS]

30

57

252%

15

S S

10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

Kogakuin University

25

S

35

Small < 8 MB (in cache)Medium < 15 MB
Large

Performance of SpMV (bandmatrix (CRS), bandwidth = 32)

x3.2

M L

M

0

10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

35

S L
S S L

S

S

L S L

5

à DD arithmetic is bounded by memory accessà Input matrix A will be given by double precisionà Data size is a half

45

S

1 thread4 threads

40

20

Kogakuin University

75%

40.0

transposed SpMV 1 threadtransposed SpMV 4 threads
SpMV 1 thread
SpMV 4 threads

35.0
30.0
25.0
20.0
15.0
10.0
5.0

10

1 thread4 threads

5
0
0.125

0.19 0.25

0.5
1Size of matrix N (x105)2

10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 2013

0.0

4

Name of matrix(nnz/row)
59

10TH INTL. CONF. on PARALLEL PROCESSING and APPLIED MATHEMATICS 201360

Comparison of Real *16 vs. Fast
Quadruple with BiCG

x 7.1

olm1000 䠄n=1000,cond=3.0E+6䠅

㻝㻚㻱㻗㻜㻟

㻝㻚㻱㻗㻜㻟

x 3.2

㻰 㻻 㼁 㻮 㻸㻱

Matrix A(CRS)㻲㻻 㻾 㼀 㻾 㻭 㻺 㻌㻽 㼁 㻭 㻰

Lis QUAD

FORTRAN QUAD

4(n+nnz)+8nnz 4(n+nnz)+8nnz

㻝㻚㻱㻙㻜㻣

㻝㻚㻱㻙㻜㻟
㻝㻚㻱㻙㻜㻡
㻝㻚㻱㻙㻜㻣
㻝㻚㻱㻙㻜㻥
㻝㻚㻱㻙㻝㻝

㻝㻚㻱㻙㻝㻟

16n

Vector x

8n

16n

16n

Workings

6*8n

6*16n

6*16n

121.9MB175.8MB221.6MBJSIAM Applied Mathematics Seminar, Dec. 27, 2013㻝㻚㻱㻙㻝㻟
㻜

㻡㻜

㻝㻜㻜
཯᚟ᅇᩘ

㻝㻡㻜

㻞㻜㻜

㻜

㻝㻜㻜

㻞㻜㻜

㻟㻜㻜
㻠㻜㻜
཯᚟ᅇᩘ

㻡㻜㻜

㻢㻜㻜

㻣㻜㻜

• Almost same accuracy (At most 10%)!
61

Convergence History of A4 with
Preconditioned BiCG

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

62

Lis-test for evaluation㻶㼍㼏㼛㼎㼕㻙
㻮㼕㻯㻳㻔㻰㻻㼁㻮㻸㻱㻕㻵㻸㼁㻔㻜㻕㻙
㻮㼕㻯㻳㻔㻰㻻㼁㻮㻸㻱㻕㻿㻿㻻㻾㻙㻮㼕㻯㻳㻔㻰㻻㼁㻮㻸㻱㻕㻵㻸㼁㻯㻙㻮㼕㻯㻳㻔㻰㻻㼁㻮㻸㻱㻕㻶㼍㼏㼛㼎㼕㻙㻮㼕㻯㻳㻔㻽㼁㻭㻰㻕㻵㻸㼁㻔㻜㻕㻙㻮㼕㻯㻳㻔㻽㼁㻭㻰㻕㻿㻿㻻㻾㻙㻮㼕㻯㻳㻔㻽㼁㻭㻰㻕
㻡㻜㻜
㻝㻜㻜㻜
㻝㻡㻜㻜
㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㼕㼠㼑㼞㼍㼠㼕㼛㼚㼟㻝㻚㻱㻙㻜㻡

㻝㻚㻱㻙㻜㻝

㻝㻚㻱㻙㻝㻝

8n㻜㻝㻚㻱㻙㻜㻟

4(n+nnz)+16nnz

8n㻝㻚㻜㻜㻱㻗㻜㻡
㻝㻚㻜㻜㻱㻗㻜㻟
㻝㻚㻜㻜㻱㻗㻜㻝
㻝㻚㻜㻜㻱㻙㻜㻝
㻝㻚㻜㻜㻱㻙㻜㻟
㻝㻚㻜㻜㻱㻙㻜㻡
㻝㻚㻜㻜㻱㻙㻜㻣
㻝㻚㻜㻜㻱㻙㻜㻥
㻝㻚㻜㻜㻱㻙㻝㻝
㻝㻚㻜㻜㻱㻙㻝㻟

㻝㻚㻱㻙㻜㻝

㻸㼕㼟㻌㻽㼁㻭㻰
㻲㻻㻾㼀㻾㻭㻺㻌㻽㼁㻭㻰㻝㻚㻱㻗㻜㻝

㻝㻚㻱㻙㻜㻥

Vector b

sum㻸㼕㼟㻌㻽㼁㻭㻰
㻲㻻㻾㼀㻾㻭㻺㻌㻽㼁㻭㻰㻝㻚㻱㻗㻜㻝㻸 㼕㼟 㻌 㻽 㼁 㻭 㻰

DOUBLE㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙rdb2048l (n=2048,cond=1.8E+3䠅┦ᑐṧᕪ䝜䝹䝮

㻝㻤㻜
㻝㻢㻜
㻝㻠㻜
㻝㻞㻜
㻝㻜㻜
㻤㻜
㻢㻜
㻠㻜
㻞㻜
㻜

┦ᑐṧᕪ䝜䝹䝮

ᐇ⾜᫬㛫㻔⛊㻕Time for 50 BiCG Iterations
Poisson 䠄n=10^6, CRS䠅, Xeon 2.8GHz㻞㻜㻜㻜• Over 2K combinations:
10 Preconditioners x 13 Solvers x 11 Storage formats x 2 precisions• Not necessary to install. Run from USB• Prepare Matrix data as text file with Matrix
Market’ exchange format• Run in parallel if the PC has multi-core• To click, solutions, history, etc are computed㻵㻸㼁㻯㻙㻮㼕㻯㻳㻔㻽㼁㻭㻰㻕

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

63

JSIAM Applied Mathematics Seminar, Dec. 27, 201364

Comparison is done easily!

Lis-test: GUI for Library Lis
JSIAM Applied Mathematics Seminar, Dec. 27, 201365

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

66

Basic Idea of Restart

Algorithms• Until Now:

(1)Solve Ax*=b with some initial value x0(2)Solve Ax=b with an initial value x*– In general, (1) and (2) have same spaces, same
methods, and same precisions– (1) and (2) have same spaces, same methods but
different precisions(combination of Double and Fast Quadruple).

JSIAM Applied Mathematics Seminar, Dec. 27, 201367

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

68

SWITCH Algorithm

PERIODIC algorithm• Restart with different precision arithmetic• A Fast Quadruple is used each k iterations– Current solution xk is passed at the restart– Upper and Lower part of Double-Double var.
are stored in different arrays– Only Upper part is used for
for(k=0;k<matitr;k++){Double Precision
The first stepif( nrm2<restart_tol ) break;– Two Stages are}
performed by
Clear all values except x
for(k=k+1;k<maxtr;k++) {Different Precisions– All values are passed at the change– No cost at the change of Q Æ D– Lower part is cleared at the change of D Æ Q㼒㼛㼞㻔㼗㻩㻜㻧㼗㻨㼙㼍㼤㼕㼠㼞㻧㼗㻗㻗㻕㻌㼧㼕㼒㻔㻌㼗㻑㼕㼚㼠㼑㼞㼢㼍㼘㻨㼚㼡㼙 㻕㻌㼧㻲㼍㼟㼠㻌㻽㼡㼍㼐㼞㼡㼜㼘㼑㻌㼕㼟㻌㼡㼟㼑㼐㼩㻌㼑㼘㼟㼑㻌㼧
㻸㼛㼣㼑㼞㻌㼜㼍㼞㼠㻌㼕㼟㻌㼏㼘㼑㼍㼞㼑㼐㻰㼛㼡㼎㼘㼑㻌㼕㼟㻌㼡㼟㼑㼐㼩The second stepif( nrm2<tol ) break;}
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

JSIAM Applied Mathematics Seminar, Dec. 27, 201369

TeoplitzȖ=1.3, n=10^5

Convergence History
Ȗ=1.3

㻟㻚㻥㻤
㻠㻚㻤㻥

㼨㼨㼎㻙㻭㼤㼨㼨
㻞㻚㻠㻣㻱㻙㻝㻜
㻞㻚㻡㻟㻱㻙㻝㻜
㻞㻚㻡㻝㻱㻙㻝㻜
㻥㻚㻟㻠㻱㻙㻝㻝

㻝㻚㻜㻱㻗㻜㻜

㻝㻚㻜㻱㻗㻜㻜

㻟㻚㻝㻢㻱㻙㻝㻜
㻞㻚㻜㻞㻱㻙㻝㻜

㻰㻻㼁㻮㻸㻱
㻽㼁㻭㻰㻰㻽㻙㻿㼃㻵㼀㻯㻴㻔䃔㻩㻝㻱㻙㻝㻝㻕㻝㻚㻜㻱㻙㻜㻞
㻝㻚㻜㻱㻙㻜㻠
㻝㻚㻜㻱㻙㻜㻢
㻝㻚㻜㻱㻙㻜㻤
㻝㻚㻜㻱㻙㻝㻜

㻰㻻㼁㻮㻸㻱
㻽㼁㻭㻰㻰㻽㻙㻿㼃㻵㼀㻯㻴㻔䃔㻩㻝㻱㻙㻤㻕㻝㻚㻜㻱㻙㻜㻞
㻝㻚㻜㻱㻙㻜㻠
㻝㻚㻜㻱㻙㻜㻢
㻝㻚㻜㻱㻙㻜㻤
㻝㻚㻜㻱㻙㻝㻜
㻝㻚㻜㻱㻙㻝㻞

㻝㻚㻜㻱㻙㻝㻞
㻜– Epsilon is restart criterion of DQ-SWITCH– Num: Quad. Ops. Used num times per 10 iterations
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

Ȗ=1.4
㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙

㼟㼑㼏㻚
㻢㻚㻢㻜
㻞㻚㻟㻟
㻝㻚㻤㻞
㻝㻚㻢㻣

㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙㻲㻹㻭㻞㼋㻿㻿㻱㻞㻿㼃㻵㼀㻯㻴 䃔㻩㻝㻚㻜㻱㻙㻜㻥䃔㻩㻝㻚㻜㻱㻙㻝㻜䃔㻩㻝㻚㻜㻱㻙㻝㻝㻼㻱㻾㻵㻻㻰㻵㻯 㼚㼡㼙㻩㻝
㼚㼡㼙㻩㻞
㼚㼡㼙㻩㻟
㼚㼡㼙㻩㻠
㼚㼡㼙㻩㻡
㼚㼡㼙㻩㻢

㼕㼠㼑㼞㻚
㼠㼛㼠㼍㼘
㼐㼛㼡㼎㼘㼑
㻝㻝㻟
㻜
㻥㻡
㻣㻠
㻥㻡
㻤㻢
㻝㻜㻟
㻝㻜㻜
㻙
㻙
㻙
㻙
㻙
㻙
㻙
㻙
㻝㻜㻣
㻡㻞
㻝㻝㻤
㻠㻢

70

㻡㻜
㻝㻜㻜
㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㼕㼠㼑㼞㼍㼠㼕㼛㼚㼟

㻝㻡㻜

㻜

㻡㻜

㻝㻜㻜

㻝㻡㻜

㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㼕㼠㼑㼞㼍㼠㼕㼛㼚㼟• DQ-SWITCH is good convergence
71

JSIAM Applied Mathematics Seminar, Dec. 27, 201372㼍㼕㼞㼒㼛㼕㼘㼋㻞㼐㻰㻻 㼁㻮㻸㻱㻽 㼁㻭㻰㻿㼃㻵㼀㻯㻴 䃔 㻩㻝㻚㻜㻱㻙㻝㻜䃔 㻩㻝㻚㻜㻱㻙㻝㻝䃔 㻩㻝㻚㻜㻱㻙㻝㻞

Epsilon dependency
䃒㻩㻝㻚㻟
㼕㼠㼑㼞㻚
䃔
㼠㼛㼠㼍㼘
㼐㼛㼡㼎㼘㼑
㼝㼡㼍㼐
㼟㼑㼏㻚
㼠㼛㼠㼍㼘
㻽㼁㻭㻰
㻝㻝㻟
㻞㻚㻤㻤
㻰㻽㻙㻿㼃㻵㼀㻯㻴 㻝㻚㻜㻜㻱㻙㻜㻟
㻝㻝㻠
㻞
㻝㻝㻞
㻞㻚㻤㻣
㻝㻚㻜㻜㻱㻙㻜㻠
㻝㻜㻥
㻝㻝
㻥㻤
㻞㻚㻡㻥
㻝㻚㻜㻜㻱㻙㻜㻡
㻝㻜㻡
㻞㻟
㻤㻞
㻞㻚㻞㻢
㻝㻚㻜㻜㻱㻙㻜㻢
㻝㻜㻠
㻟㻡
㻢㻥
㻞㻚㻜㻝
㻝㻚㻜㻜㻱㻙㻜㻣
㻥㻡
㻠㻣
㻠㻤
㻝㻚㻡㻤
㻝㻚㻜㻜㻱㻙㻜㻤
㻥㻠
㻢㻝
㻟㻟
㻝㻚㻞㻥
㻝㻚㻜㻜㻱㻙㻜㻥
㻥㻡
㻣㻠
㻞㻝
㻝㻚㻜㻤
㻝㻚㻜㻜㻱㻙㻝㻜
㻥㻡
㻤㻢
㻥
㻜㻚㻤㻢
㻝㻚㻜㻜㻱㻙㻝㻝
㻝㻜㻟
㻥㻤
㻡
㻜㻚㻤㻠

䃒㻩㻝㻚㻠
㼕㼠㼑㼞㻚
㼐㼛㼡㼎㼘㼑
㼝㼡㼍㼐
㼟㼑㼏㻚
㻝㻡㻡
㻟㻚㻥㻠
㻝㻡㻢
㻞
㻝㻡㻠
㻟㻚㻥㻠
㻝㻡㻞
㻝㻡
㻝㻟㻣
㻟㻚㻢㻞
㻝㻠㻢
㻟㻝
㻝㻝㻡
㻟㻚㻝㻢
㻝㻟㻤
㻠㻣
㻥㻝
㻞㻚㻢㻣
㻝㻞㻟
㻢㻡
㻡㻤
㻝㻚㻥㻢
㻝㻝㻥
㻤㻟
㻟㻢
㻝㻚㻡㻟
㻙㻙
㻙㻙
㻙㻙

• Choice of appropriate epsilon is important• Small epsilon reduces much computation time
• Smaller epsilon makes divergence
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

73

㻝㻜㻜㻜

㻝㻡㻜㻜

㼨㼨㼎㻙㻭㼤㼨㼨
㻟㻚㻡㻞㻱㻙㻝㻜
㻝㻚㻠㻥㻱㻙㻝㻜
㻝㻚㻡㻥㻱㻙㻝㻜
㻥㻚㻞㻞㻱㻙㻝㻝
㻝㻚㻠㻢㻱㻙㻝㻜
㼨㼨㼎㻙㻭㼤㼨㼨
㻞㻚㻥㻢㻱㻙㻜㻥
㻠㻚㻞㻡㻱㻙㻝㻝
㻝㻚㻣㻝㻱㻙㻝㻜
㻠㻚㻞㻜㻱㻙㻝㻜
㻠㻚㻠㻣㻱㻙㻝㻜

74

Poisson 䠄n=10^6, CRS䠅, Xeon 2.8GHz
x 7.1

㻝㻤㻜
㻝㻢㻜
㻝㻠㻜
㻝㻞㻜
㻝㻜㻜
㻤㻜
㻢㻜
㻠㻜
㻞㻜
㻜

x 3.2
x 1.2
㻰 㻻 㼁 㻮 㻸㻱

㻡㻜㻜

㼟㼑㼏㻚
㻞㻚㻜㻟
㻣㻚㻟㻝
㻟㻚㻢㻣
㻞㻚㻠㻞
㻞㻚㻟㻞
㼟㼑㼏㻚
㻟㻚㻠㻞
㻝㻜㻚㻡㻟
㻠㻚㻡㻣
㻠㻚㻜㻣
㻠㻚㻝㻤

Computation Time㻰㻻㼁㻮㻸㻱
㻽㼁㻭㻰㻰㻽㻙㻿㼃㻵㼀㻯㻴㻔䃔㻩㻝㻱㻙㻠㻕㻜㼨㼨㼎㻙㻭㼤㼨㼨
㻟㻚㻞㻡㻱㻙㻜㻤
㻡㻚㻟㻢㻱㻙㻝㻜
㻟㻚㻝㻡㻱㻙㻝㻜
㻟㻚㻝㻟㻱㻙㻝㻜
㻟㻚㻡㻢㻱㻙㻝㻜JSIAM Applied Mathematics Seminar, Dec. 27, 2013

Computationtimeᐇ⾜᫬㛫㻔⛊㻕㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙

㼣㼍㼚㼓㻟
㼠㼛㼠㼍㼘
㼐㼛㼡㼎㼘㼑
㻰㻻㼁㻮㻸㻱
㻠㻣㻢
㻠㻣㻢
㻽㼁㻭㻰
㻟㻣㻞
㻿㼃㻵㼀㻯㻴 䃔㻩㻝㻚㻜㻱㻙㻝㻜㻠㻢㻜㻟㻢㻝䃔㻩㻝㻚㻜㻱㻙㻝㻝㻠㻡㻥㻠㻠㻠䃔㻩㻝㻚㻜㻱㻙㻝㻞㻠㻣㻥㻠㻣㻢㼘㼍㼚㼓㼡㼍㼓㼑
㼠㼛㼠㼍㼘㼐㼛㼡㼎㼘㼑
㻰㻻㼁㻮㻸㻱
㻟㻥
㻟㻥
㻽㼁㻭㻰
㻟㻢
㻿㼃㻵㼀㻯㻴 䃔㻩㻝㻚㻜㻱㻙㻝㻜㻟㻤㻟㻠䃔㻩㻝㻚㻜㻱㻙㻝㻝㻟㻣㻟㻡䃔㻩㻝㻚㻜㻱㻙㻝㻞㻠㻜㻟㻥– İ is restaring criterion of SWITCH

㼟㼑㼏㻚
㻝㻤㻚㻢㻠
㻢㻥㻚㻟㻥
㻞㻠㻚㻞㻡
㻞㻝㻚㻢㻢
㻞㻞㻚㻤㻣

• QUAD and SWITCH improve 2 digits for solution’ quality• SWITCH is 20% overhead on the double, however robust

A4: Electronics Effect㻝㻚㻜㻱㻗㻜㻤
㻝㻚㻜㻱㻗㻜㻢
㻝㻚㻜㻱㻗㻜㻠
㻝㻚㻜㻱㻗㻜㻞
㻝㻚㻜㻱㻗㻜㻜
㻝㻚㻜㻱㻙㻜㻞
㻝㻚㻜㻱㻙㻜㻠
㻝㻚㻜㻱㻙㻜㻢
㻝㻚㻜㻱㻙㻜㻤
㻝㻚㻜㻱㻙㻝㻜
㻝㻚㻜㻱㻙㻝㻞

㼕㼠㼑㼞㻚
㼠㼛㼠㼍㼘
㼐㼛㼡㼎㼘㼑
㻠㻡㻢㻣
㻠㻡㻢㻣
㻟㻤㻟㻤
㻠㻠㻜㻞
㻠㻜㻥㻝
㻠㻟㻟㻝
㻠㻝㻣㻢
㻠㻣㻜㻥
㻠㻡㻢㻣

Mixed.

㻸 㼕㼟 㻌 㻽 㼁 㻭 㻰㻲㻻 㻾 㼀 㻾 㻭 㻺 㻌㻽 㼁 㻭 㻰㻞㻜㻜㻜㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㼕㼠㼑㼞㼍㼠㼕㼛㼚㼟
JSIAM Applied Mathematics Seminar, Dec. 27, 201375

JSIAM Applied Mathematics Seminar, Dec. 27, 201376

Auto Restart with Different Precisions

Auto Restart of DQ-SWITCH• Convergent history shows three patterns:< nrm2_min )• Compute deviation if( nrm2
nrm2_min = nrm2; x_bak = x;
= nrm2;of residual norm nrm_bak[k%10]
if( k>=10 ) {– 䠄C)Converge– 䠄D)Diverge– 䠄S)Stagnate㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙㻝㻚㻜㻱㻗㻜㻜

㻰㻻㼁㻮㻸㻱
㻽㼁㻭㻰㻰㻽㻙㻿㼃㻵㼀㻯㻴㻔䃔㻩㻝㻱㻙㻝㻝㻕㻝㻚㻜㻱㻙㻜㻞
㻝㻚㻜㻱㻙㻜㻠
㻝㻚㻜㻱㻙㻜㻢

• 䠄D)
• 䠄S)

㻝㻚㻜㻱㻙㻜㻤
㻝㻚㻜㻱㻙㻝㻜
㻝㻚㻜㻱㻙㻝㻞

• To Detect 䠄D) and (S),restart at the point

㻜

㻡㻜
㻝㻜㻜
㻺㼡㼙㼎㼑㼞㻌㼛㼒㻌㼕㼠㼑㼞㼍㼠㼕㼛㼚㼟v = 0.0; c = 0;
for(i=0;i<10;i++) {t = nrm_bak[i] - nrm_bak[(k-9)%10];t = t / nrm_bak[(k-9)%10];v = v + t*t;if( nrm_bak[(k-9)%10] <= nrm_bak[i] )c = c+1;}v = v / 10;if( v<=0.1 || (c==10 && v>=100) ) break;if( nrm2<tol ) break;v t 10 2
v d 10 1㻝㻡㻜77Electoronics BiCG with ILU(0)㻾㼑㼘㼍㼠㼕㼢㼑㻌㼞㼑㼟㼕㼐㼡㼍㼘㻌㻞㻙㼚㼛㼞㼙

2}

JSIAM Applied Mathematics Seminar, Dec. 27, 2013㻝㻚㻜㻱㻗㻜㻢
㻝㻚㻜㻱㻗㻜㻠
㻝㻚㻜㻱㻗㻜㻞
㻝㻚㻜㻱㻗㻜㻜
㻝㻚㻜㻱㻙㻜㻞
㻝㻚㻜㻱㻙㻜㻠
㻝㻚㻜㻱㻙㻜㻢
㻝㻚㻜㻱㻙㻜㻤
㻝㻚㻜㻱㻙㻝㻜
㻝㻚㻜㻱㻙㻝㻞

v

1 p § nrm(i )  nrm(1) ·¦ ¨ nrm(1) ¸¸
p i 1 ¨©
¹

JSIAM Applied Mathematics Seminar, Dec. 27, 201378

Mixed Precision Iterative Methods• Complicated problems are solved with Mixed or QUAD.
• Overhead of the mixed precision iterative methods is20%㻰㻻㼁㻮㻸㻱
㻽㼁㻭㻰㻰㻽㻙㻿㼃㻵㼀㻯㻴㻔䃔㻩㻝㻱㻙㻠㻕㻭㻰㻽㻙㻿㼃㻵㼀㻯㻴
㻜㻚㻜

㻞㻚㻜

㻠㻚㻜

㻢㻚㻜

• SWITCH is Good at least 2 digits with 20% more– D Æ Q: easy, robust, however depends on timing of restart
㻤㻚㻜

㻝㻜㻚㻜

• Auto restart of DQ-SWITH㻝㻞㻚㻜– Deviation is used to detect “Diverge” or “Stagnate”

㻱㼤㼑㼏㼡㼠㼕㼛㼚㻌㼠㼕㼙㼑㼟㻌㻔㼕㼚㻌㼟㼑㼏㼛㼚㼐㼟㻕• Divergence or Stagnation is detected.
• Computation time is reduced.
JSIAM Applied Mathematics Seminar, Dec. 27, 201379

JSIAM Applied Mathematics Seminar, Dec. 27, 201380

Test Problems

Parallel Issues for Fast Quad.

• Poisson• Depends on the implementation of Ax,
ATx, M-1x, M-Tx, and Matrix Storage
Format• Data transfered is almost same
• Heavy Computation
Æ Suitable for Distributed Parallel• Less round-off errors
Æ lighter preconditioner (easy to parallelize)– 2 dimension, FDM– N=10^6, nnz=5x10^6

rdb2048l• rdb2048l (Chemical engineering)– MatrixMarket, n=2048, nnz=12032, cond = 1.8x10^3• olm1000 (Hydrodynamics)– MatrixMarket, n=1000, nnz=3996, cond = 3x10^6olm1000• A4 (Electronic potential)•– n=23,994, nnz=214,060
Cryg10000 (CRYSTAL GROWTH EIGENMODES )– UF Sparse Matrix Collection, n=10000, nnz=49699

A4• circuit_3 (Circuit Simulation)– n=12,127, nnz=48,137
cryg10000

JSIAM Applied Mathematics Seminar, Dec. 27, 201381

• University of Florida Sparse Matrix Collection
dimension

airfoil_2d
wang3
language

14,21426,064399,130Size of memory
Double Lis Quad
3.9MB
4.7MB
259,688
3.7MB
5.1MB
177,168
61.1MB
1,216,334 39.8MB

nnz

JSIAM Applied Mathematics Seminar, Dec. 27, 201382

Application Program Interface

Comparison of Double and DQ-SWITCH

MatrixJSIAM Applied Mathematics Seminar, Dec. 27, 201383

• Data Types (Precision)• Matrix Storage Format• Algorithms• Function names• Computing Environments

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

84

The traditional way of using libraries
1. Preparation of matrices and vectors using
library-specific data structures
2. Function calls with a function's name andits arguments in a prescribed order
As a result...
• User programs will depend on a specificlibrary

SILC– Not easy to replace the library by another
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

85

SILC: Simple Interface for Library
Collections

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

86

Solving a system of linear equations
Ax=b• In the traditional way (using LAPACK in C)• Basic ideas

double *A, *b;
int kl, ku, lda, ldb, nrhs, info, *ipiv;– Data transfer and a request for computation– Mathematical expressions for the request– A separate memory space for the computation

dgbtrf (N, N, kl, ku, A, lda, ipiv, &info); /* LU factorization */if (info == 0)
dgbtrs ('N', N, kl, ku, nrhs, A, lda, ipiv, b, ldb, &info); /* solve */Input data
User program"x = Aᨘb"

Results• In SILC

Separate
memory space

silc_envelope_t A, b, x;
SILC_PUT ("A", &A);
SILC_PUT ("b", &b);
SILC_EXEC ("x = Aᨘb"); /* call a solver (e.g., dgbtrf & dgbtrs) */SILC_GET (&x, "x");Library collections

JSIAM Applied Mathematics Seminar, Dec. 27, 201387

JSIAM Applied Mathematics Seminar, Dec. 27, 201388

SILC servers in different computing
environments

Main benefits of using SILC
• Source-level independence between user
programs and matrix computation libraries• A user program (client) that solves Ax = b

PC

C

S

– Easy access to alternative solvers and matrix
storage formats, possibly in other libraries
PC– Instant porting to other computing environments C
without any modification in user programs

PC• You need to prepare only the smallest
amount of data

– Where A is a tridiagonal matrix in the CRS format– Run in the notebook PC of Environment (a)– In a 100-Base TX local-area network
Environment

S

PC

SMP

C

S

Specification(a) A notebook PC(b) SGI Altix3700– Temporary buffers are automatically allocated• Language-independent mathematical expressions– Applicable in many programming languages (C, Fortran,
Python, MATLAB)

JSIAM Applied Mathematics Seminar, Dec. 27, 2013

• SILC servers in (c) and (d) achieved better performance
because of parallel computation
Execution time (in seconds)

10,000C(d) SGI Altix3700

16 threads

Same as (b)

S

C

S

C

S

JSIAM Applied Mathematics Seminar, Dec. 27, 201390

• Data structures– Data types: scalar, vector, matrix, cubic array– Precisions: integer, real, complex (single or double)– Matrix storage formats: dense, banded, CRS– Binary arithmetic operators (+, ෥, *, /, %)– Solutions of systems of linear equations (Aᨘb)– Conjugate transposes (A'), complex conjugates (A~)– Built-in functions

C

S

C

S

100

10

C

• Mathematical expressions

S(c) OpenPower 710 (4 threads)(d) Altix3700 (16 threads)

1 thread

S

Functionalities– Data size: 0.46MB (N=10,000) to 4.27MB (N=80,000)(b) Altix3700 (1 thread)Intel Itanium2 1.3GHz u 32,
32GB memory, Red Hat Linux
Advanced Server 2.1

C

4 threads

89• About 0.1 second of data communications over the LAN

1,000

N/A

(c) IBM eServer
IBM Power5 1.65GHz u 2
OpenPower 710 (4 logical CPUs), 1GB memory,
SuSE Linux Enterprise Server 9

Experimental results(a) Notebook PC

OpenMPIntel Pentium M 733 1.1GHz,
768MB memory,
Fedora Core 3• Ex. "sqrt(b' * b)" is the 2-norm of vector bC1
10,00020,00040,000– Subscript

S• Ex. "A[1:5,1:5]" is a 5u5 submatrix of A

80,000Dimension N

JSIAM Applied Mathematics Seminar, Dec. 27, 201391

JSIAM Applied Mathematics Seminar, Dec. 27, 201392Collaborators and Acknowledgement

Conclusion• Accurate Computation– Powerful tool for “Iterative Methods”– Another choice for designing Algorithm– Tool for analysis• MuPAT: Ease of Use of D-D and Q-D• Lis: Iterative Solvers with Fast D-D• Lis-test: the simplest tool
• SILC: General Purpose API
JSIAM Applied Mathematics Seminar, Dec. 27, 2013

• Lis & Lis-test
H. Kotakemori
A. Fujii (Kogakuin U)
K. Nakajima (U Tokyo)
A. Nishida(U Kyushu)• SILCT. Kajiyama (Universidade
Nova de Lisboa)
A. Nukada (TITECH)
R. Suda (U Tokyo)
A. Nishida (U Kyushu)• Acceleration on AVX
T. Hishinuma (Kogakuin U.)
K. Asakawa (Kogakuin U)
A. Fujii (Kogakuin U)
T. Tanaka (Kogakuin U)• MuPAT
T. Saitoh (Tokyo U of
Science)
S. Kikkawa (TUS)
E. Ishiwata (TUS)

Lis and SILC are parts of SSI project which is funded by JST/CREST
93

JSIAM Applied Mathematics Seminar, Dec. 27, 201394

Appendix A :Algorithms for DD and QD Arithmetics

Algorithm A.4, Two-Prod, produces a nonoverlapping expansion p + e such that a × b = p + e,
where p is an approximation to a × b and e represents the round-oﬀ error in the calculation of p,in [14, p. 326].

We describe the details of the algorithms for
DD and QD arithmetics. The procedures of algorithms are based on Knuth [5], Dekker [3], Priest[9], Shewchuk [14], Bailey [1] and Hida et al. [4].

A.1

Algorithm A.4 Two-Prod(a, b)

Preliminaries for DD and QD arithmetics

In this section, we introduce some algorithms of
ﬂoating-point arithmetic.
Assuming that | a | ≥ | b |, Algorithm A.1, FastTwo-Sum, produces a nonoverlapping expansions + e such that a + b = s + e, where s is an approximation to a + b and e represents the round-oﬀ
error in the calculation of s, in [14, p. 312].

1: p ← a ⊗ b
2: [ah , al ] ←Split(a)3: [bh , bl ] ←Split(b)4: e ← ((ah ⊗ bh ⊖ p) ⊕ ah ⊗ bl ⊕ al ⊗ bh ) ⊕ al ⊗ bl
5: return (p, e)

Algorithm A.5, Two-Sqr, produces a nonoverlapping expansion p+e such that a2 = p+e, where
p is an approximation to a2 and e represents theround-oﬀ error in the calculation of p.
Algorithm A.5 Two-Sqr(a)1: p ← a ⊗ a
2: [ah , al ] ←Split(a)3: e ← ((ah ⊗ ah ⊖ p) ⊕ (ah ⊗ al ) ⊗ 2.0) ⊕ al ⊗ al
4: return (p, e)

Algorithm A.6, Three-Sum, produces a nonoverlapping expansion d + e + f such that a + b + c =d + e + f , in [4].

Algorithm A.2, Two-Sum, is similar to Algorithm A.1, but Algorithm A.2 does not require
the condition of | a | ≥ | b |.

Algorithm A.6 Three-Sum(a, b, c)

Algorithm A.1 Fast-Two-Sum(a, b) : Assume that

1: [t0 , t1 ] ← Two-Sum(a, b)
2: [d, t2 ] ← Two-Sum(t0 , c)
3: [e, f ] ← Two-Sum(t1 , t2 )4: return (d, e, f )|a| ≥ |b|
1: s ← a ⊕ b
2: v ← s ⊖ a
3: e ← b ⊖ v
4: return (s, e)

Algorithm A.7, Three-Sum2, produces two doubleprecision numbers d = (a ⊕ b) ⊕ c and e = (a + b +c) − s, in [4].

Algorithm A.2 Two-Sum(a, b)
1: s ← a ⊕ b
2: v ← s ⊖ a
3: e ← (a ⊖ (s ⊖ v)) ⊕ (b ⊖ v)
4: return (s, e)

Algorithm A.7 Three-Sum2(a, b, c)
1: [t0 , t1 ] ← Two-Sum(a, b)
2: [d, t2 ] ← Two-Sum(t0 , c)
3: e = t1 ⊕ t2
4: return (d, e)

Algorithm A.3, Split, produces a 26 bit value
ah and a nonoverlapping 26 bit value al such that| ah | > | al | and a = ah + al , in [14, p. 325].

Supposing that a0 , a1 , a2 , a3 and a4 construct a
ﬁve-term expansion with limited overlapping bits,with a0 being the most signiﬁcant component. Then
Algorithm A.8, Renormalize, produces a four-term
nonoverlapping expansion b(qd) = b0 + b1 + b2 + b3 .

Algorithm A.3 Split(a)1: t ← 134217729 ⊗ a2: v ← t ⊖ a
3: ah ← t ⊖ v
4: al ← a ⊖ ah
5: return (ah , al )A–1

double number a(dd) and returns the double-double
number c(dd) = c0 + c1 . If you want to add adouble-double number b(dd) to a double-precision
number a, dd d add (b0 , b1 , a) returns the result.
d dd add is same as dd d add.

Algorithm A.8 Renormalize(a0 , a1 , a2 , a3 , a4 )1: [s, t3 ] ← Fast-Two-Sum(a3 , a4 )2: [s, t2 ] ← Fast-Two-Sum(a2 , s)
3: [s, t1 ] ← Fast-Two-Sum(a1 , s)
4: [b0 , t0 ] ← Fast-Two-Sum(a0 , s)
5: [s, t2 ] ← Fast-Two-Sum(t2 , t3 )6: [s, t1 ] ← Fast-Two-Sum(t1 , s)
7: [b1 , t0 ] ← Fast-Two-Sum(t0 , s)
8: [s, t1 ] ← Fast-Two-Sum(t1 , t2 )9: [b2 , t0 ] ← Fast-Two-Sum(t0 , s)
10: b3 = t0 ⊕ t1
11: return (b0 , b1 , b2 , b3 )Algorithm A.10 dd d add (a0 , a1 , b)
1: [s, e] ← Two-Sum(a0 , b)
2: e ← e ⊕ a1
3: [c0 , c1 ] ← Fast-Two-Sum(s, e)
4: return (c0 , c1 )Algorithm A.9, Renormalize2, produces a fourterm nonoverlapping expansion b(qd) = b0 + b1 +b2 + b3 . This algorithm is similar to Algorithm
A.8 except for the number of arguments.

Algorithm A.11, dd dd add, shows the procedure for adding a double-double number b(dd) to adouble-double number a(dd) and returns the doubledouble number c(dd) = c0 + c1 .

Algorithm A.9 Renormalize2(a0 , a1 , a2 , a3 )Algorithm A.11 dd dd add (a0 , a1 , b0 , b1 )1: [s, t2 ] ← Fast-Two-Sum(a2 , a3 )2: [s, t1 ] ← Fast-Two-Sum(a1 , s)
3: [b0 , t0 ] ← Fast-Two-Sum(a0 , s)
4: [s, t1 ] ← Fast-Two-Sum(t1 , t2 )5: [b1 , t0 ] ← Fast-Two-Sum(t0 , s)
6: [b2 , b3 ] ← Fast-Two-Sum(t0 , t1 )7: return (b0 , b1 , b2 , b3 )1: [s, e] ← Two-Sum(a0 , b0 )2: e ← e ⊕ (a1 ⊕ b1 )3: [c0 , c1 ] ← Fast-Two-Sum(s, e)
4: return (c0 , c1 )A.2.2
Table 11 shows the number of double-precision
arithmetic operations for Algorithm A.1 ∼ A.9.

Algorithm A.12 (A.13), dd d sub (d dd sub),shows the procedure for subtracting a double- precision number b (a double-double number b(dd) )from a double-double number a(dd) (a double-precision number a) and returns the double-double
number c(dd) = c0 + c1 .

Table 11: Number of double-precision arithmetic
operations for Algorithm A.1 ∼ A.9
Algorithm
Fast-Two-Sum (A.1)
Two-Sum (A.2)
Split (A.3)
Two-Prod (A.4)
Two-Sqr (A.5)
Three-Sum (A.6)
Three-Sum2 (A.7)
Renormalize (A.8)
Renormalize2 (A.9)

A.2

⊕, ⊖
3
6
3
10
7
18
13
28
18

⊗
0
0
1
7
5
0
0
0
0

Total
3
6
4
17
12
18
13
28
18

Algorithm A.12 dd d sub (a0 , a1 , b)
1: s ← −b
2: [c0 , c1 ] ← dd d add (a0 , a1 , s)
3: return (c0 , c1 )Algorithm A.13 d dd sub (a, b0 , b1 )1: s0 ← −b0
2: s1 ← −b1
3: [c0 , c1 ] ← d dd add (a, s0 , s1 )4: return (c0 , c1 )Algorithms for DD arithmetic

In this section, we show the algorithms of fourarithmetic operations for double-double numbers.

A.2.1

subtraction

Algorithm A.14, dd dd sub, shows the procedure for subtracting a double-double number b(dd)from a double-double number a(dd) and returns
the double-double number c(dd) = c0 + c1 .

addition

Algorithm A.10, dd d add, shows the procedure
for adding a double-precision number b to a doubleA–2

Algorithm A.14 dd dd sub (a0 , a1 , b0 , b1 )Algorithm A.18 d dd div (a, b0 , b1 )1: s0 ← −b0
2: s1 ← −b1
3: [c0 , c1 ] ← dd dd add (a0 , a1 , s0 , s1 )4: return (c0 , c1 )1: c0 ← a ⊘ b0
2: [p, e] ← Two-Prod(c0 , b0 )3: c1 ← ((a ⊖ p) ⊖ e ⊖ c ⊗ b1 ) ⊘ b0
4: [c0 , c1 ] ← Fast-Two-Sum(c0 , c1 )5: return (c0 , c1 )A.2.3

multiplication

Supposing that b0 ̸= 0. Algorithm A.19, dd dd
div, shows the procedure for dividing a doubledouble number a(dd) by a double-double numberb(dd) and returns the double-double number c(dd) =c0 + c1 .

Algorithm A.15, dd d mul, shows the procedure
for multiplying a double-double number a(dd) by adouble-precision number b and returns the doubledouble number c(dd) = c0 + c1 . d dd mul is sameas dd d mul.

Algorithm A.19 dd dd div (a0 , a1 , b0 , b1 )Algorithm A.15 dd d mul (a0 , a1 , b)1: c0 ← a0 ⊘ b0
2: [p, e] ← Two-Prod(c0 , b0 )3: c1 ← ((a0 ⊖ p) ⊖ e ⊕ a1 ⊖ c ⊗ b1 ) ⊘ b0
4: [c0 , c1 ] ← Fast-Two-Sum(c0 , c1 )5: return (c0 , c1 )1: [p, e] ← Two-Prod(a0 , b)
2: e ← e ⊕ (a1 ⊗ b)
3: [c0 , c1 ] ← Fast-Two-Sum(p, e)
4: return (c0 , c1 )Algorithm A.16, dd dd mul, shows the procedure for multiplying a double-double number a(dd)by a double-double number b(dd) and returns the
double-double number c(dd) = c0 + c1 .

Table 12 shows the number of double precision
arithmetic operations for double-double arithmetic.

A.3

Algorithms for QD arithmetic

Algorithm A.16 dd dd mul (a0 , a1 , b0 , b1 )A.3.1

addition

1: [p, e] ← Two-Prod(a0 , b0 )2: e ← e ⊕ (a0 ⊗ b1 )3: e ← e ⊕ (a1 ⊗ b0 )4: [c0 , c1 ] ← Fast-Two-Sum(p, e)
5: return (c0 , c1 )A.2.4

Algorithm A.20, qd d add, shows the procedure
for adding a double-precision number b to a quaddouble number a(qd) and returns the quad-double
number c(qd) = c0 + c1 + c2 + c3 .

division

Algorithm A.20 qd d add (a0 , a1 , a2 , a3 , b)

Supposing that b ̸= 0 and b0 ̸= 0. Algorithm
A.17 (A.18), dd d div (d dd div), shows the procedure for dividing a double-double number a(dd)(a double-precision number a) by a double-precision
number b (a double-double number b(dd) ) and returns the double-double number c(dd) = c0 + c1 .
Algorithm A.17 dd d div (a0 , a1 , b)1: c0 ← a0 ⊘ b
2: [p, e] ← Two-Prod(c0 , b)
3: c1 ← ((a0 ⊖ p) ⊖ e ⊕ a1 ) ⊘ b
4: [c0 , c1 ] ← Fast-Two-Sum(c0 , c1 )5: return (c0 , c1 )A–3

1: [c0 , e] ← Two-Sum(a0 , b)
2: [c1 , e] ← Two-Sum(a1 , e)
3: [c2 , e] ← Two-Sum(a2 , e)
4: [c3 , e] ← Two-Sum(a3 , e)
5: [c0 , c1 , c2 , c3 ] ← Renormalize(c0 , c1 , c2 , c3 , e)
6: return (c0 , c1 , c2 , c3 )Algorithm A.21, qd dd add, shows the procedure for adding a double-double number b(dd) to aquad-double number a(qd) and returns the quaddouble number c(qd) = c0 + c1 + c2 + c3 .

Table 12: Number of double precision arithmetic operations for double-double arithmetic

Addition
Subtraction
Multiplication
Division

Algorithm
dd d add
dd dd add
dd d sub, d dd sub
dd dd subdd d mul
dd dd mul
dd d div
d dd div
dd dd div

⊕&⊖
10
11
10
11
14
15
16
16
17

⊗
0
0
0
0
8
9
7
8
8

⊘
0
0
0
0
0
0
2
2
2

Total
10
11
10
11
22
24
25
26
27

Algorithm A.21 qd dd add (a0 , a1 , a2 , a3 , b0 , b1 )Algorithm A.23 qd d sub (a0 , a1 , a2 , a3 , b)

1: [c0 , e0 ] ← Two-Sum(a0 , b0 )2: [c1 , e1 ] ← Two-Sum(a1 , b1 )3: [c1 , e0 ] ← Two-Sum(c1 , e0 )4: [c2 , e1 ] ← Two-Sum(a2 , e1 )5: [c2 , e0 ] ← Two-Sum(c2 , e0 )6: [e0 , e1 ] ← Two-Sum(e0 , e1 )7: [c3 , e0 ] ← Two-Sum(a3 , e0 )8: e0 = e0 ⊕ e1
9: [c0 , c1 , c2 , c3 ] ← Renormalize(c0 , c1 , c2 , c3 , e0 )10: return (c0 , c1 , c2 , c3 )1: s ← −b
2: [c0 , c1 , c2 , c3 ] ← qd d add (a0 , a1 , a2 , a3 , s)
3: return (c0 , c1 , c2 , c3 )Algorithm A.24 d qd sub (a, b0 , b1 , b2 , b3 )1: s0 ← −b0
2: s1 ← −b1
3: s2 ← −b2
4: s3 ← −b3
5: [c0 , c1 , c2 , c3 ] ← d qd add (a, s0 , s1 , s2 , s3 )6: return (c0 , c1 , c2 , c3 )Algorithm A.22, qd qd add, shows the procedure for adding a quad-double number b(qd) to aquad-double number a(qd) and returns the quaddouble number c(qd) = c0 + c1 + c2 + c3 .

Algorithm A.25 (A.26), qd dd sub (dd qd sub),shows the procedure for subtracting a double-double number b(dd) (a quad-double number b(qd) ) froma quad-double number a(qd) (a double-double number a(dd) ) and returns the quad-double numberc(qd) = c0 + c1 + c2 + c3 .

Algorithm A.22 qd qd add (a0 , a1 , a2 , a3 , b0 , b1 , b2 , b3 )1: [c0 , e0 ] ← Two-Sum(a0 , b0 )2: [c1 , e1 ] ← Two-Sum(a1 , b1 )3: [c2 , e2 ] ← Two-Sum(a2 , b2 )4: [c3 , e3 ] ← Two-Sum(a3 , b3 )5: [c1 , e0 ] ← Two-Sum(c1 , e0 )6: [c2 , e0 , e1 ] ← Three-Sum(c2 , e1 , e0 )7: [c3 , e0 ] ← Three-Sum2(c3 , e2 , e0 )8: e0 = e0 ⊕ e1 ⊕ e3
9: [c0 , c1 , c2 , c3 ] ← Renormalize(c0 , c1 , c2 , c3 , e0 )10: return (c0 , c1 , c2 , c3 )Algorithm A.25 qd dd sub (a0 , a1 , a2 , a3 , b0 , b1 )1: b0 ← −b0
2: b1 ← −b1
3: [c0 , c1 , c2 , c3 ] ← qd dd add (a0 , a1 , a2 , a3 , b0 , b1 )4: return (c0 , c1 , c2 , c3 )Algorithm A.26 dd qd sub (a0 , a1 , b0 , b1 , b2 , b3 )A.3.2

subtraction

Algorithm A.23 (A.24), qd d sub (d qd sub), shows the procedure for subtracting a double-precision number b (a quad-double number b(qd) ) froma quad-double number a(qd) (a double-precision
number a) and returns the quad-double numberc(qd) = c0 + c1 + c2 + c3 .
A–4

1: b0 ← −b0
2: b1 ← −b1
3: b2 ← −b2
4: b3 ← −b3
5: [c0 , c1 , c2 , c3 ] ← dd qd add (a0 , a1 , b0 , b1 , b2 , b3 )6: return (c0 , c1 , c2 , c3 )Algorithm A.27 shows the procedure for subtracting a quad-double number b(qd) from a quad-double number a(qd) and returns the quad-double
number c(qd) = c0 + c1 + c2 + c3 .

by a quad-double number b(qd) and returns thequad-double number c(qd) = c0 + c1 + c2 + c3 .

Algorithm A.27 qd qd sub (a0 , a1 , a2 , a3 , b0 , b1 , b2 , b3 )Algorithm A.30 qd qd mul(a0 , a1 , a2 , a3 , b0 , b1 , b2 , b3 )1: b0 ← −b0
2: b1 ← −b1
3: b2 ← −b2
4: b3 ← −b3
5: [c0 , c1 , c2 , c3 ] ← qd qd add (a0 , a1 , a2 , a3 , b0 , b1 , b2 , b3 )6: return (c0 , c1 , c2 , c3 )1: [p0 , q0 ] ← Two-Prod(a0 , b0 )2: [p1 , q1 ] ← Two-Prod(a0 , b1 )3: [p2 , q2 ] ← Two-Prod(a1 , b0 )4: [p1 , p2 , q0 ] ← Three-Sum(p1 , p2 , q0 )5: [p3 , q3 ] ← Two-Prod(a0 , b2 )6: [p4 , q4 ] ← Two-Prod(a1 , b1 )7: [p5 , q5 ] ← Two-Prod(a2 , b0 )8: [p2 , q1 , q2 ] ← Three-Sum(p2 , q1 , q2 )9: [p3 , p4 , p5 ] ← Three-Sum(p3 , p4 , p5 )10: [p2 , p3 ] ← Two-Sum(p2 , p3 )11: [p4 , q1 ] ← Two-Sum(p4 , q1 )12: [p3 , p4 ] ← Two-Sum(p3 , p4 )13: p4 ← q2 ⊕ p5 ⊕ q1 ⊕ p4
14: p3 ← p3 ⊕ (a0 ⊗ b3 ) ⊕ (a1 ⊗ b2 ) ⊕ (a2 ⊗ b1 ) ⊕ (a3 ⊗ b0 ) ⊕q0 ⊕ q3 ⊕ q4 ⊕ q5
15: [c0 , c1 , c2 , c3 ] ← Renormalize(p0 , p1 , p2 , p3 , p4 )16: return (c0 , c1 , c2 , c3 )A.3.3

multiplication

Algorithm A.28, qd d mul, shows the procedure
for multiplying a quad-double number a(qd) by adouble-precision number b and returns the quaddouble number c(qd) = c0 + c1 + c2 + c3 .
Algorithm A.28 qd d mul (a0 , a1 , a2 , a3 , b)
1: [p0 , q0 ] ← Two-Prod(a0 , b)
2: [p1 , q1 ] ← Two-Prod(a1 , b)
3: [p2 , q2 ] ← Two-Prod(a2 , b)
4: [p1 , q0 ] ← Two-Sum(p1 , q0 )5: [p2 , q0 , q1 ] ← Three-Sum(p2 , q0 , q1 )6: p3 ← a3 ⊗ b
7: [p3 , q0 ] ← Three-Sum2(p3 , q0 , q2 )8: q0 ← q0 ⊕ q1
9: [c0 , c1 , c2 , c3 ] ← Renormalize(p0 , p1 , p2 , p3 , q0 )10: return (c0 , c1 , c2 , c3 )A.3.4

division

Supposing that b ̸= 0 and b0 ̸= 0. Algorithm
A.31 (A.32), qd d div (d qd div), shows the procedure for dividing a quad-double number a(qd) (adouble-precision number a) by a double-precision
number b (a quad-double number b(qd) ) and returns the quad-double number c(qd) = c0 + c1 +c2 + c3 .

Algorithm A.29, qd dd mul, shows the procedure for multiplying a quad-double number a(qd)by a double-double number b(dd) and returns thequad-double number c(qd) = c0 + c1 + c2 + c3 .

Algorithm A.31 qd d div (a0 , a1 , a2 , a3 , b)1: c0 ← a0 ⊘ b
2: [t0 , t1 ] ← Two-Prod(c0 , b)
3: [r0 , r1 , r2 , r3 ] ← qd dd sub (a0 , a1 , a2 , a3 , t0 , t1 )4: c1 ← r0 ⊘ b
5: [t0 , t1 ] ← Two-Prod(c1 , b)
6: [r0 , r1 , r2 , r3 ] ← qd dd sub (r0 , r1 , r2 , r3 , t0 , t1 )7: c2 ← r0 ⊘ b
8: [t0 , t1 ] ← Two-Prod(c2 , b)
9: [r0 , r1 , r2 , r3 ] ← qd dd sub (r0 , r1 , r2 , r3 , t0 , t1 )10: c3 ← r0 ⊘ b
11: [c0 , c1 , c2 , c3 ] ← Renormalize2(c0 , c1 , c2 , c3 )12: return (c0 , c1 , c2 , c3 )Algorithm A.29 qd dd mul (a0 , a1 , a2 , a3 , b0 , b1 )1: [p0 , q0 ] ← Two-Prod(a0 , b0 )2: [p1 , q1 ] ← Two-Prod(a0 , b1 )3: [p2 , q2 ] ← Two-Prod(a1 , b0 )4: [p3 , q3 ] ← Two-Prod(a2 , b0 )5: [p4 , q4 ] ← Two-Prod(a1 , b1 )6: [p1 , p2 , q0 ] ← Three-Sum(p1 , p2 , q0 )7: [p2 , p3 , p4 ] ← Three-Sum(p2 , p3 , p4 )8: [q1 , q2 ] ← Two-Sum(q1 , q2 )9: [p2 , q1 ] ← Two-Sum(p2 , q1 )10: [p3 , q2 ] ← Two-Sum(p3 , q2 )11: [p3 , q1 ] ← Two-Sum(p3 , q1 )12: p4 ← p4 ⊗ q2 ⊗ q1
13: q3 ← q3 ⊕ q4 ⊕ (a3 ⊗ b0 ) ⊕ (a2 ⊗ b1 )14: [p3 , q0 ] ← Three-Sum2(p3 , q0 , q3 )15: p4 ← p4 ⊕ q0
16: [c0 , c1 , c2 , c3 ] ← Renormalize(p0 , p1 , p2 , p3 , p4 )17: return (c0 , c1 , c2 , c3 )Algorithm A.30, qd qd mul, shows the procedure for multiplying a quad-double number a(qd)A–5

Algorithm A.32 d qd div (a, b0 , b1 , b2 , b3 )Algorithm A.35 qd qd div (a0 , a1 , a2 , a3 , b0 , b1 , b2 , b3 )1: c0 ← a ⊘ b0
2: [t0 , t1 , t2 , t3 ] ← d qd mul (c0 , b0 , b1 , b2 , b3 )3: [r0 , r1 , r2 , r3 ] ← d qd sub (a, t0 , t1 , t2 , t3 )4: c1 ← r0 ⊘ b0
5: [t0 , t1 , t2 , t3 ] ← d qd mul (c1 , b0 , b1 , b2 , b3 )6: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )7: c2 ← r0 ⊘ b0
8: [t0 , t1 , t2 , t3 ] ← d qd mul (c2 , b0 , b1 , b2 , b3 )9: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )10: c3 ← r0 ⊘ b0
11: [c0 , c1 , c2 , c3 ] ← Renormalize2(c0 , c1 , c2 , c3 )12: return (c0 , c1 , c2 , c3 )1: c0 ← a0 ⊘ b0
2: [t0 , t1 , t2 , t3 ] ← d qd mul (c0 , b0 , b1 , b2 , b3 )3: [r0 , r1 , r2 , r3 ] ← qd qd sub (a0 , a1 , a2 , a3 , t0 , t1 , t2 , t3 )4: c1 ← r0 ⊘ b0
5: [t0 , t1 , t2 , t3 ] ← d qd mul (c1 , b0 , b1 , b2 , b3 )6: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )7: c2 ← r0 ⊘ b0
8: [t0 , t1 , t2 , t3 ] ← d qd mul (c2 , b0 , b1 , b2 , b3 )9: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )10: c3 ← r0 ⊘ b0
11: [c0 , c1 , c2 , c3 ] ← Renormalize2(c0 , c1 , c2 , c3 )12: return (c0 , c1 , c2 , c3 )Supposing that b0 ̸= 0. Algorithm A.33 (A.34),
qd dd div (dd qd div), shows the procedure fordividing a quad-double number a(qd) (a doubledouble number a(dd) ) by a double-double numberb(dd) (a quad-double number b(qd) ) and returns thequad-double number c(qd) = c0 + c1 + c2 + c3 .

arithmetic operations for quad-double arithmetic.
One quad-double arithmetic operation needs tens
or hundreds of double-precision operations, thenthe computation time may require hundreds times
greater than double-precision arithmetic.

A.4

Algorithm A.33 qd dd div (a0 , a1 , a2 , a3 , b0 , b1 )1: c0 ← a0 ⊘ b0
2: [t0 , t1 ] ← d dd mul (c0 , b0 , b1 )3: [r0 , r1 , r2 , r3 ] ← qd dd sub (a0 , a1 , a2 , a3 , t0 , t1 )4: c1 ← r0 ⊘ b
5: [t0 , t1 ] ← d dd mul (c1 , b0 , b1 )6: [r0 , r1 , r2 , r3 ] ← qd dd sub (r0 , r1 , r2 , r3 , t0 , t1 )7: c2 ← r0 ⊘ b
8: [t0 , t1 ] ← d dd mul (c2 , b0 , b1 )9: [r0 , r1 , r2 , r3 ] ← qd dd sub (r0 , r1 , r2 , r3 , t0 , t1 )10: c3 ← r0 ⊘ b
11: [c0 , c1 , c2 , c3 ] ← Renormalize2(c0 , c1 , c2 , c3 )12: return (c0 , c1 , c2 , c3 )The other algorithms for QD

Algorithm A.36, qd sqr, shows the procedure
for squaring a quad-double number a(qd) and returns the quad-double number c(qd) = c0 + c1 +c2 + c3 .
Algorithm A.36 qd sqr (a0 , a1 , a2 , a3 )Algorithm A.34 dd qd div (a0 , a1 , b0 , b1 , b2 , b3 )1: c0 ← a0 ⊘ b0
2: [t0 , t1 , t2 , t3 ] ← d qd mul (c0 , b0 , b1 , b2 , b3 )3: [r0 , r1 , r2 , r3 ] ← dd qd sub (a0 , a1 , t0 , t1 , t2 , t3 )4: c1 ← r0 ⊘ b0
5: [t0 , t1 , t2 , t3 ] ← d qd mul (c1 , b0 , b1 , b2 , b3 )6: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )7: c2 ← r0 ⊘ b0
8: [t0 , t1 , t2 , t3 ] ← d qd mul (c2 , b0 , b1 , b2 , b3 )9: [r0 , r1 , r2 , r3 ] ← qd qd sub (r0 , r1 , r2 , r3 , t0 , t1 , t2 , t3 )10: c3 ← r0 ⊘ b0
11: [c0 , c1 , c2 , c3 ] ← Renormalize2(c0 , c1 , c2 , c3 )12: return (c0 , c1 , c2 , c3 )Supposing that b0 ̸= 0. Algorithm A.35, qd qd
div, shows the procedure for dividing a quad-double
number a(qd) by a quad-double number b(qd) andreturns the quad-double number c(qd) = c0 + c1 +c2 + c3 .
Table 13 shows the number of double precision
A–6

1: [p0 , q0 ] ← Two-Sqr(a0 )2: [p1 , q1 ] ← Two-Prod(a0 , a1 )3: p1 ← p1 ⊗ 2.0
4: q1 ← q1 ⊗ 2.0
5: [p2 , q2 ] ← Two-Prod(a0 , a2 )6: p2 ← p2 ⊗ 2.0
7: q2 ← q2 ⊗ 2.0
8: [p3 , q3 ] ← Two-Sqr(a1 )9: [p1 , q0 ] ← Two-sum(p1 , q0 )10: [q0 , q1 ] ← Two-sum(q0 , q1 )11: [p2 , p3 ] ← Two-sum(p2 , p3 )12: [s0 , t0 ] ← Two-Sum(q0 , p2 )13: [s1 , t1 ] ← Two-Sum(q1 , p3 )14: [s1 , t0 ] ← Two-Sum(s1 , t0 )15: t0 ← t0 ⊕ t1
16: [s1 , t0 ] ← Fast-Two-Sum(s1 , t0 )17: [p2 , t0 ] ← Fast-Two-Sum(s0 , s1 )18: [p3 , q0 ] ← Fast-Two-Sum(t1 , t0 )19: p4 ← 2.0 ⊗ a0 ⊗ a3
20: p5 ← 2.0 ⊗ a1 ⊗ a2
21: [p4 , p5 ] ← Two-Sum(p4 , p5 )22: [q2 , q3 ] ← Two-Sum(q2 , q3 )23: [t0 , t1 ] ← Two-Sum(p4 , q2 )24: t1 ← t1 ⊕ p5 ⊕ q3
25: [p3 , p4 ] ← Two-Sum(p3 , t0 )26: p4 ← p4 ⊕ q0 ⊕ t1
27: [c0 , c1 , c2 , c3 ] ← Renormalize(p0 , p1 , p2 , p3 , p4 )28: return (c0 , c1 , c2 , c3 )

Table 13: Number of double precision arithmetic operations for quad-double arithmetic
Addition

Subtraction

Multiplication

Division

Algorithm
qd d add
qd dd add
qd qd add
qd d sub, d qd sub
qd dd sub, dd qd sub
qd qd subqd d mul
qd dd mulqd qd mul
qd d divd qd div
qd dd div
dd qd div
qd qd div

⊕, ⊖
52
71
91
52
71
91
96
154
171
261
540
273
569
579

⊗
0
0
0
0
0
0
22
35
46
21
66
24
66
66

⊘
0
0
0
0
0
0
0
0
0
4
4
4
4
4

Total
52
71
91
52
71
91
118
189
217
286
610
301
639
649

Supposing that n is a power of 2. Algorithm
A.37, mul pwr dd, shows the procedure for multiplying each component of double-double number
by n.

Supposing that n is an integer. Algorithm A.39
shows the procedure for computing an n-th powerof a quad-double number a(qd) and returns thequad-double number.

Algorithm A.37 mul pwr dd (a0 , a1 , n)Algorithm A.39 qd pow (a0 , a1 , a2 , a3 , n) : Assume

1: b0 ← a0 ⊗ n
2: b1 ← a1 ⊗ n
3: return (b0 , b1 )that n ≥ 0
1: if n = 0 then2:
return 1.0
3: end if
4: if n = 1 then5:
return (a0 , a1 , a2 , a3 )6: end if
7: r0 ← a0
8: r1 ← a1
9: r2 ← a2
10: r3 ← a3
11: s0 ← 1.0
12: N ← n
13: while N > 0 do
14:
if N is odd then15:[s0 , s1 , s2 , s3 ]← qd qd mul(s0 , s1 , s2 , s3 , r0 , r1 , r2 , r3 )16:
end if17:
N ←N ⊘2
if N > 0 then18:
19:
[r0 , r1 , r2 , r3 ] ← qd sqr (r0 , r1 , r2 , r3 )end if20:
21: end while
22: return (r0 , r1 , r2 , r3 )Supposing that n is a power of 2. Algorithm
A.38, mul pwr qd, shows the procedure for multiplying each component of quad-double number
by n.
Algorithm A.38 mul pwr qd (a0 , a1 , a2 , a3 , n)1: b0 ← a0 ⊗ n
2: b1 ← a1 ⊗ n
3: b2 ← a2 ⊗ n
4: b3 ← a3 ⊗ n
5: return (b0 , b1 , b2 , b3 )A–7

References (Not full):[1] T. J. Dekker, A Floating-point technique for extending the available precision, Numerische Mathematik,Vol. 18, pp. 224-242 (1971) .
[2] D. H. Bailey, QD, MPFUN, http://crd.lbl.gov/˜dhbailey/mpdist/[3] Y. Hida, X. S. Li and D. H. Bailey, Quad-double arithmetic: algorithms, implementation, and application,
Technical Report LBNL-46996, Lawrence Berkeley National Laboratory, Berkeley, CA 94720 (2000).
[4] Y. Hida, X. S. Li and D. H. Bailey, Algorithms for quad-double precision floating point arithmetic, Proc. of
the 15th IEEE Symposium on Computer Arithmetic, pp. 155-162 (2001).
[5] R. Brent, A fortran multiple-precision arithmetic package, ACM TOMS, Vol. 4, No. 1, pp. 71-81 (1978).
[6] J. D. Hogg and J. A. Scott, A fast robust mixed precision solver for the solution of sparse symmetriclinear systems, Technical Report, RAL-TR-2008-023 (2008).
[7] T. Saito, E. Ishiwata, and H. Hasegawa, Development of quadruple precision arithmetic toolbox QuPAT
on Scilab, Proceedings of the 2010 International Conference on Computational Science and its Applications(ICCSA 2010), Part II, Lecture Notes in Computer Science 6017, pp. 60-70, Springer-Verlag (2010).
[8] T. Saito, E. Ishiwata and H. Hasegawa, Analysis of the GCR method with mixed precision arithmetic
using QuPAT, Journal of Computational Science, Volume 3, Issue 3, pp. 87-91, Elsevier (2012).
[9] S. Kikkawa, T. Saito, E. Ishiwata and H. Hasegawa, Development and acceleration of multiple precision
arithmetic toolbox MuPAT for Scilab, JSIAM Letters Vol. 5, pp.9-12 (2013).
[10] MuPAT: http://www.mi.kagu.tus.ac.jp/qupat.html[11] Hisashi Kotakemori, Hidehiko Hasegawa, and Akira Nishida, Perfomance Evaluation of Parallel
Iterative Method Library using OpenMP, Proc. of the 8th International Conference on High Performance
Computing in Asia Pacific Region (HPC Asia 2005), p. 432-436, Nov., 2005, Beijing, China.
[12] Hisashi Kotakemori, Akihiro Fujii, Hidehiko Hasegawa, and Akira Nishida, Implementation of Fast Quad
Precision Operation and Acceleration with SSE2 for Iterative Solver Library, IPSJ Transactions on
advanced Computing Systems, Vol. 1, No. 1, pp. 73-84 (2008) in Japanese.
[13] Toshiaki Hishinuma, Akihiro Fujii, Teruo Tanaka, and Hidehiko Hasegawa, AVX Acceleration of Sparse
Matrix-Vector Multiplication in Double-Double, Proc. on High Performance and Computing Symposium
2013, pp. 23-31 (2013) in Japanese.
[14] R. Barrett, M. Berry, T. F. Chan, J.W. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romineand H. Van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods,
2nd Edition, SIAM (1994).
[15] Lis: http://www.ssisc.org/lis/index.en.html[16] Tamito Kajiyama, Akira Nukada, Hidehiko Hasegawa, Reiji Suda, and Akira Nishida, SILC: A Flexible
and Environment Independent Interface to Matrix Computation Libraries, Lecture Notes in Computer
Science 3911, pp. 928-935, Springer-Verlag (2006) (PPAM2005).
[17] SILC: http://www.ssisc.org/silc/index.html