Introduction to verified computation

September 2015

Waseda University
Graduate School of Fundamental Science and Engineering
Major in Pure and Applied Mathematics
Research on Numerical Analysis

Kouta Sekine

c Copyright by Kouta Sekine, 2015⃝
All Rights Reserved.

ii

Contents

Chapter 1. Preliminary . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1. Usage of gfortran . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2. Mathematical Notations . . . . . . . . . . . . . . . . . . . . . . . . .

2

Chapter 2. What’s verified computation? . . . . . . . . . . . . . . .

5

2.1. Numerical Computations and Its Error . . . . . . . . . . . . . . . . .

6

2.2. Verified Computations and Interval Arithmetic on Real Numbers . .

8Chapter 3. Floating-Point Number and Its Interval Arithmetic

11

3.1. IEEE 754 : Standard for Floating-Point Arithmetic . . . . . . . . . .

123.2. Interval Arithmetic on Floating-Point Numbers . . . . . . . . . . . .

14

3.3. How to Change Rounding Rules . . . . . . . . . . . . . . . . . . . .

14

Chapter 4. Verification Theory for System of Linear Equations

19

Chapter 5. How to Install Some Packages . . . . . . . . . . . . . . .

25

5.1. Install BLAS and Lapack . . . . . . . . . . . . . . . . . . . . . . . .

26

5.2. Install Slib . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Chapter 6. Programming of Verified Theory . . . . . . . . . . . . .

31

6.1. Interval Arithmetic for Matrix on floating-point number . . . . . . .

32

6.2. Programming of Verified Theory for System of Linear Equations . . .

34

iii

Chapter 1

Preliminary11.1. Usage of gfortranIn this lecture, we will use the gfortran which is the fortran compiler. If you havenever used the Fortran language, you try it now!
Exercise 1.1. Write the Algorithm 1 and run it. Let’s filename be “hello.f90”.
Algorithm 1 hello.f90
program hello
write(*,*) ’Hello world!’stop
end program hello

You can compile the following command for “hello.f90” by using gfortran on the
Gnome terminal:$ gfortran hello.f90
Finally, you execute the program using the following command:$ ./a.out
If display “Hello world!”, then succeed.

1.2. Mathematical Notations
In this section, we prepare some mathematical notations.
• Let R be the set of real numbers.
• The closed interval denoted by [a, b] which is the set of real numbers givenby [a, b] = {x ∈ R|a ≤ x ≤ b}.
• Let I denote the n × n identity matrix.
• Let O denote the n × n matrix of all zeros and 0 denote the n-vector of all
zeros.
• Let e the n-vector of all ones i.e. e := (1, . . . , 1)T ∈ Rn .
2

• Inequalities for matrices are understood componentwise, e.g. for real n × nmatrices A = (aij ) and B = (bij ) the notation A ≤ B means aij ≤ bij for all(i, j).
• The notation |A| means |A| = (|aij |) ∈ Rn×n , the nonnegative matrix consisting of componentwise absolute values of A.
• Similar notation is applied to real vectors.
• Let ∥ · ∥ denote the norm of vector and matrix.
• For a n dimensional vector x = (x1 , · · · , xn )T , the maximum norm is definedby ∥x∥∞ := max1≤i≤n |xi |.
• For a m × n matrix A = (aij ) , the maximum norm is defined by ∥A∥∞ :=∑max1≤i≤m nj=1 |aij |.
• Let IR be set of real intervals.
• Let F be set of floating-point numbers defined by IEEE 754 standard.
• Let IF be set of intervals for floating-point numbers.

3

Chapter 2

What’s verified computation?

5

2.1. Numerical Computations and Its Error
Numerical computation solved to some problems. For example, system of linear
equationsAx = b, A ∈ Rn×n , b ∈ Rn ,(1)

we can eﬃciently obtain an approximate solution x̃ of (1) using some software packages1, even if n is 10000. Of course, I do not want to solve to this problem by handcalculation!! However, when floating-point arithmetic is used for numerical computations, then the solution includes rounding errors. In general, rounding errors are
small. However, in some cases solutions are aﬀected by rounding error. Let’s feelrounding errors.

Exercise 2.1. Write the Algorithm 2 and run it.

Algorithm 2 sqrtpow.f90
program sqrtpow
integer :: i, nreal(8) :: xwrite(*,*) ’Please input a number n=’read(*,*) nx = 2d0 ! Substitute 2 × 100 for x, and ! sign means comment out in fortran.
write(*,*) ’Exact : ’, xdo i=1,nx = sqrt(x) ! sqrt(x) means the square root of x.
end dodo i=1,nx = x**2d0 ! x**2d0 means the square of x.
end dowrite(*,*) ’Approximate : ’, x
end program sqrtpow

1In this lecture, we will use the Lapack(Linear Algebra PACKage) with the BLAS(Basic Linear

Algebra Subprogram). See Section 5.1.

6

Example of execution:$ gfortran sqrtpow.f90$ ./a.out
Please input a number n =0
Exact : 2.0000000000000000
Approximate : 2.0000000000000000$ ./a.out
Please input a number n =25
Exact : 2.0000000000000000
Approximate : 2.0000000066771721$ ./a.out
Please input a number n =50
Exact : 2.0000000000000000
Approximate : 1.6487212645509468$ ./a.out
Please input a number n =55
Exact : 2.0000000000000000
Approximate : 1.0000000000000000

7

Algorithm 2 first repeatedly run the square root of x by Do statement. Next, theprogram repeatedly run the square of x by Do statement. If real numbers, exactly
the same as before. However, because we use floating-point numbers, results include
rounding errors.

2.2. Verified Computations and Interval Arithmetic on Real
Numbers
Numerical solutions include rounding errors. Therefore, we do not know how
accurate computed solutions are. Verified computations solve this problem. The
essence of verified computation is interval. For example, floating-point numbers can’t
describe 1/3 = 0.3333 · · · . In verified computation, we describe interval [0.33, 0.34]
which enclose 1/3. Since intervals contain the exact solutions, we can understandhow the accuracy of approximate solutions.
How do you calculate four arithmetic operators in intervals?

Definition 2.2 (Interval Arithmetic on real numbers). Let IR be set of realintervals. For x = [x, x̄] ∈ IR and y = [y, ȳ] ∈ IR, we define the following intervalarithmetic:x + y = [x + y, x̄ + ȳ]
x − y = [x − ȳ, x̄ − y]x · y = [min(x̄ · ȳ, x̄ · y, x · ȳ, x · y), max(x̄ · ȳ, x̄ · y, x · ȳ, x · y)]x/y = [min(x̄/ȳ, x̄/y, x/ȳ, x/y), max(x̄/ȳ, x̄/y, x/ȳ, x/y)] for 0 ̸∈ y
For ◦ = {+, −, ·, /}, Definition 2.2 satisfyx̃ ◦ ỹ ∈ x ◦ y, ∀x̃ ∈ x, ∀ỹ ∈ y.
8

Example 2.3. Let x = [−2, 4] and y = [1, 2]. Calculate four arithmetic operators
using the interval arithmetic.
x + y = [−2, 4] + [1, 2] = [−2 + 1, 4 + 2] = [−1, 6]
x − y = [−2, 4] − [1, 2] = [−2 − 2, 4 − 1] = [−4, −3]
x · y = [−2, 4] · [1, 2] = [min(−2, −4, 4, 8), max(−2, −4, 4, 8)] = [−4, 8]
x/y = [−2, 4]/[1, 2] = [min(−2, −1, 4, 2), max(−2, −1, 4, 2)] = [−2, 4]
Exercise 2.4. Let x = [2, 4] and y = [−2, −1]. Calculate four arithmetic operators using the interval arithmetic.

9Chapter 3

Floating-Point Number and Its
Interval Arithmetic

11

3.1. IEEE 754 : Standard for Floating-Point Arithmetic
The IEEE 754 : Standard for Floating-Point Arithmetic is a technical standard for
floating-point number. Floating-point unit of many CPU use the IEEE 754 standard.
For more detail, see
http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4610935
In this chapter, I will introduce to formats of floating-point number in IEEE 754.
Next, I will describe rounding rules.
IEEE 754 standard have single and double formats which are binary floating-point
basic formats. Representations of floating-point numbers in the binary formats are(−1)s × 2e × m, wherea) s is 1-bit sign (0 or 1)b) e is any integer satisfying emin ≤ e ≤ emax.
c) m is a number represented by the form d0 .d1 d2 · · · dp−1 where di is 0 or 1.
The value of emax and p are given in Table 3.1 (emin shall be 1 − emax).
Table 3.1. Verification results.
Parameter single double
p
24
53
emax
127
1023

In Figure 3.1, we displays the representation of Single format. We note that startof the significand field is d1 because we can put d0 = 1 i.e. normalization.

Sign sExponent e (8bit)

Significand field p-1 (23bit)

s e0 e1 e2 e3 e4 e5 e6 e7 d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13d14d15 d16d17 d18d19d20d21 d22d23

Figure 3.1. Representation of Single format.

12

Example 3.1. Convert from the following 10 decimal numbers to the floating
point number. Here, p = 4 and e have 3 bit.
7.25 = (111.01)2= (−1)0 × 22 × (1.1101)2
where (·)2 means 2 decimal numbers. Answer is
0 010 1101
Let F be set of floating-point numbers defined by IEEE 754 standard. You can
notice that if p = 3 in Example 3.1 , then we can not represent the floating point
number. Here, we have F ⊂ R. IEEE 754 standard defines five rounding rules. In
this lecture, we introduce following three rounding rules:a) Rounding to Nearest :・Rounding to Nearest rounds to the nearest value.
・This operator is denoted by 2 : R → Fb) Rounding toward +∞ :・Rounding toward +∞ directed rounding towards positive infinity.
・This operator is denoted by △ : R → Fc) Rounding toward −∞ :䕱䕧・Rounding toward −∞ directed rounding towards negative infinity.
・This operator is denoted by ▽ : R → F

From Figure 3.2 , the real number r is enclosed by interval [▽r,△r].
䕱r

䕧r

䕱r

0

䕧r

0

Real number r

Real number r

Figure 3.2. Rounding (· is a floating-point
number).
䕧䕱

13

3.2. Interval Arithmetic on Floating-Point Numbers
In Section 2.2, we define the interval arithmetic on real numbers. However, insome cases, real numbers are not represented floating-point numbers. We now extend
Definition 2.2 to use a computer.

Definition 3.2 (Interval Arithmetic on floating-point numbers). Let IF be setof intervals for floating-point numbers. For x = [x, x̄] ∈ IF and y = [y, ȳ] ∈ IF, wedefine the following interval arithmetic:x + y = [▽(x + y), △(x̄ + ȳ)]x − y = [▽(x − ȳ), △(x̄ − y)]x·y= [min(▽(x̄ · ȳ), ▽(x̄ · y), ▽(x · ȳ), ▽(x · y)), max(△(x̄ · ȳ), △(x̄ · y), △(x · ȳ), △(x · y))]x/y= [min(▽(x̄/ȳ), ▽(x̄/y), ▽(x/ȳ), ▽(x/y)), max(△(x̄/ȳ), △(x̄/y), △(x/ȳ), △(x/y))] for 0 ̸∈ y

3.3. How to Change Rounding Rules
If your fortran compiler support the oﬃcial Fortran 2003 standard, then you canchange rounding rules using fortran. For example, after Intel Fortran Compiler version
13 and gfortran version 5.1 are. If not support Fortran 2003 standard, then we change
rounding rules using the C language. In this lecture, we will change rounding rules
using the C language.
We now show Algorithm 3 which include setround function and getround function.
The setround function is setting rounding rules and this argument is −1, 0 and 1.
setround(0) : Rounding to Nearest
setround(1) : Rounding toward + ∞setround(−1) : Rounding toward − ∞14The getround function returns current rounding rules (but not setting). Return valueis:
0 : Rounding to Nearest
1 : Rounding toward + ∞−1 : Rounding toward − ∞There functions can use the fortran language.
Algorithm 3 cround.c#include<stdio.h>#include<fenv.h>
void setround(int a){
if (a == -1) fesetround(FE DOWNWARD);
else if (a == 0) fesetround(FE TONEAREST);else if (a == 1) fesetround(FE UPWARD);
else printf(”ERROR : setround : Please input -1, 0, 1¥n”);}
void getround(int *a){volatile double e, x, y, z;
e = 10e-30;
x = 1.0 + e;
y = 1.0 - e;
if (x == y) *a = 0;else {z = (-1.0) + e;
if (x == 1 ||z == -1) *a = -1;
else if (y == 1) *a = 1;
else *a = 2;}}
void setround (int *a){setround(*a);}
void getround (int *a){getround(a);}15We next make static library which is cround.a file using following commands:$ gcc − c cround.c$ ar cr cround.a cround.o
Exercise 3.3. Write the Algorithm 4 and run it. Here, you need to use thefollowing command:$ gfortran roundtest.f90 cround.aAlgorithm 4 roundtest.f90
program roundtest
real(8) :: a, b, c
integer :: rnd
call setround(0)a = -1d0
b = 3d0
call getround(rnd)write(*,*) ’Rounding rule is ’, rndcall setround(1)c = a/b
write(*,*) ’Rounding to +infinity’
write(*,’(b64)’) ccall setround(-1)
c = a/b
write(*,*) ’Rounding to -infinity’
write(*,’(b64)’) ccall setround(0)stop
end program roundtest

16

Example of execution:$ gfortran roundtest.f90 cround.a$ ./a.out
Rounding rule is 0
Rounding to + infinity
1011111111010101010101010101010101010101010101010101010101010101
Rounding to − infinity
1011111111010101010101010101010101010101010101010101010101010110

17

Chapter 4

Verification Theory for Systemof Linear Equations

19

In this chapter, we will introduce verification theories for system of linear equations
Ax = b, A ∈ Rn×n , b ∈ Rn .

(2)

We first present the Banach perturbation lemma.
Lemma 4.1 (Banach’s perturbation lemma). Let X and Y be Banach spaces. A
bounded linear operator D : X → Y has the bounded inverse operator D−1 : Y → X.
If a bounded linear operator B : X → Y holds∥D−1 B∥ < 1,(3)

then the bounded linear operator D + B : X → Y has the bounded inverse operator(D + B)−1 : Y → X and we have∥(D + B)−1 ∥ ≤∥D−1 ∥
1 − ∥D−1 B∥(4)

We next present the following theorem for regularity of a matrix A ∈ Rn×n .
Lemma 4.2. Let I denote the n × n identity matrix. Let A be an n × n real matrix
and R be some approximate inverse of A. If∥RA − I∥ < 1(5)

is satisfied, then (RA)−1 exists and we have∥(RA)−1 ∥ ≤1
1 − ∥RA − I∥

(6)

Proof
We use Lemma 4.1 as D = I and B = RA − I. Since D−1 B = I(RA − I), thecondition (3) in Lemma 4.1 is the same as the condition (5). Therefore, the matrix
RA(= D + B) has the inverse matrix (RA)−1 and we have (6).
2
20

We next present the verification theorem for system of linear equations.
Theorem 4.3. Let I denote the n × n identity matrix. Let A be an n × n realmatrix and b an n real vector. Let x∗ be a exact solution of Ax = b and x̃ an its
approximate solution. Let R be some approximate inverse of A. If∥RA − I∥ < 1
is satisfied, then the exact solution of Ax = b exists and we have∥x∗ − x̃∥ ≤∥R(b − Ax̃)∥.
1 − ∥RA − I∥

(7)

Proof
We rewrite (2) as
RAx = Rb.

(8)

Since RA is nonsingular from Lemma 4.2, (8) has an exact solutionx∗ = (RA)−1 Rb.
We havex∗ − x̃ = (RA)−1 Rb − x̃ = (RA)−1 (Rb − RAx̃)= (RA)−1 R(b − Ax̃).
Finally, from (6), we have∥x∗ − x̃∥ = ∥(RA)−1 R(b − Ax̃)∥≤ ∥(RA)−1 ∥∥R(b − Ax̃)∥≤∥R(b − Ax̃)∥.
1 − ∥RA − I∥221Remark 4.4. From Theorem 4.6,∥x∗ − x̃∥ ≤∥R(b − Ax̃)∥∥R∥∥b − Ax̃∥≤.
1 − ∥RA − I∥
1 − ∥RA − I∥

(9)The notation |A| means |A| = (|aij |) ∈ Rn×n , the nonnegative matrix consistingof componentwise absolute values of A. Similar notation is applied to real vectors.
Finally, we present an componentwise error bounds |x∗ − x̃| for Ax = b. Before that
we show the following lemma.
Lemma 4.5. Let e the n-vector of all ones i.e. e := (1, . . . , 1)T ∈ Rn . For anyx ∈ Rn ,|x| ≤ ∥x∥∞ e

(10)

Proof
Form the definition of infinity norm, we have|xj | ≤ ∥x∥∞ , 1 ≤ j ≤ n.
Therefore,|x| ≤ ∥x∥∞ e.
2
Theorem 4.6. Let I denote the n × n identity matrix. Let A be an n × n realmatrix and b an n real vector. Let x∗ be a exact solution of Ax = b and x̃ an its
approximate solution. Let R be some approximate inverse of A. If∥RA − I∥ < 1
is satisfied, then the exact solution of Ax = b exists and we have|x∗ − x̃| ≤ |R(b − Ax̃)| +∥R(b − Ax̃)∥∞|RA − I|e.
1 − ∥RA − I∥∞
22

(11)

Proof
We havex∗ − x̃ = RAx̃ − RAx̃ + Rb − Rb + x∗ − x̃
= R(b − Ax̃) + RAx̃ − Rb + x∗ − x̃
= R(b − Ax̃) + RAx̃ − RAx∗ + x∗ − x̃
= R(b − Ax̃) − RA(x∗ − x̃) + x∗ − x̃
= R(b − Ax̃) + (I − RA)(x∗ − x̃)
From (10) and (11),|x∗ − x̃| = |R(b − Ax̃) + (I − RA)(x∗ − x̃)|≤ |R(b − Ax̃)| + |(I − RA)(x∗ − x̃)|≤ |R(b − Ax̃)| + |I − RA||x∗ − x̃|≤ |R(b − Ax̃)| + |RA − I|∥x∗ − x̃∥∞ e≤ |R(b − Ax̃)| +∥R(b − Ax̃)∥∞|RA − I|e1 − ∥RA − I∥∞223Chapter 5

How to Install Some Packages

25

5.1. Install BLAS and Lapack
BLAS (Basic Linear Algebra Subprograms) are routines that provide basic vector
and matrix operations (e.g. dot product and matrix multiplication). Because BLAS
are eﬃcient and portable, linear algebra software use this (e.g. Lapack). Now, we
can choose some BLAS:a) Reference BLAS :・It is a freely-available software package but slow.

b) Intel MKL :・It is very fast.
・You need to charge.
・MATLAB use this.

c) ATLAS :・It is free.
・ATLAS is faster than Reference BLAS.
・However, I feel slow.

d) Open BLAS :・It is free.
・Open BLAS is as fast as Intel MKL.

In this lecture, we install the Open BLAS.
1) Download :・You can download zip file or ter.gz file from the following URL.

http://www.openblas.net/
2) Unzip :・You need to unzip to Home directory.

3) Edit Make.rule file :・You need to edit the Make.rule file.

26

(Before)# Fortran compiler. Default is g77# FC = gfortran⇓(After)# Fortran compiler. Default is g77
FC = gfortran(Before)# If you don’t need CBLAS interface, please comment it in.
# NO CBLAS = 1⇓(After)# If you don’t need CBLAS interface, please comment it in.
NO CBLAS = 1(Before)# If you don’t need LAPACK, please comment it in.
# If you set NO LAPACK=1, the library automatically sets NO LAPACKE=1.
# NO LAPACK = 1⇓(After)# If you don’t need LAPACK, please comment it in.
# If you set NO LAPACK=1, the library automatically sets NO LAPACKE=1.
NO LAPACK = 1

27

(Before)# If you need to synchronize FP CSR between threads (for x86/x86 64 only).
# CONSISTENT FPCSR = 1⇓(After)# If you need to synchronize FP CSR between threads (for x86/x86 64 only).
CONSISTENT FPCSR = 1(Before)# gfortran option for LAPACK# enable this flag only on 64bit Linux and if you need a thread safe lapack library# FCOMMON OPT = -frecursive⇓(After)# gfortran option for LAPACK# enable this flag only on 64bit Linux and if you need a thread safe lapack library
FCOMMON OPT = -frecursive4) make command :・You type “make” in Gnome terminal.

5) Move the libopenblas ???.so file:・In home directory, you make a directory that name “lib”.
・libopenblas ???.so file in OpenBLAS directory is. moved to lib directory.
・libopenblas ???.so file change to name “oblas.so”.

28

Next, we install Lapack. Lapack provides rotines for solving system of linear
equations, eigenvalue problems, and singular value problems.
1) Download :・You can download tgz file from the following URL.

http://www.netlib.org/lapack/2) Unzip :・You need to unzip to Home directory.

3) Edit make.inc.example file :・You may edit the make.inc.example file.

(Before)BLASLIB = ../../librefblas.a⇓(After)BLASLIB = /̃lib/oblas.so4) Change the filename of make.inc.example:
・make.inc.example file change to name “make.inc”.

5) make command :・You type “make” in Gnome terminal.

6) Move the liblapack.a file:・liblapack.a file in Lapack directory is. moved to lib directory.
・liblapack.a file change to name “lapack.a”.

29

5.2. Install Slib
Slib is not yet open my interval arithmetic library.
1) Download :・You can download zip file from the following URL (Japanese site...).

http://271.jp/slib/index.php
・Here, you need the following account and password.

Account:guest
Password:niigata
2) Unzip :・You need to unzip to Home directory.

3) make command :・You type “make” in Gnome terminal.

4) Move the libslib.a file:・libslib.a file in Slib directory is. moved to lib directory.
・libslib.a file change to name “slib.a”.

30Chapter 6

Programming of Verified Theory31

6.1. Interval Arithmetic for Matrix on floating-point number
In Section 3.2, we define the interval arithmetic on floating-point number. We
extend Definition 3.2 to use a matrix.

Definition 6.1 (Add and Subtract : Interval Arithmetic for matrix). Let IF be
set of intervals for floating-point numbers. For X = [X, X̄] ∈ IFn×m and Y = [Y, Ȳ ] ∈IFn×m , we define the following interval arithmetic:
X + Y = [▽(X + Y), △(X̄ + Ȳ )]X − Y = [▽(X − Ȳ ), △(X̄ − Y)]

Next, we want to show the matrix multiplication. However, matrix multiplication
is calculated by many add and multiply. Therefor, we define midrad form of intervals.
Let x be an interval of real number. Real numbers xm and xr are defined byxm =x̄ − x+ x,
2

xr = xm − x.
Then, we rewrite the interval x as
⟨xm , xr ⟩ = [x, x̄].
Here, xm and xr means midpoint and radius of interval x, respectively.
Next, let x be an interval of floating-point number. floating-point numbers xm andxr are defined by(
xm = △)x̄ − x+x ,2xr = △ (xm − x) ,32

where △(·) means that all calculation in curly bracket is upward. Here, we have[x, x̄] ⊂ ⟨xm , xr ⟩ ⊂ [▽(xm − xr ), △(xm + xr )].
The midrad form is applied to matrix multiplication for interval.
Definition 6.2 (Add and Subtract : Interval Arithmetic for matrix). Let IF
be set of intervals for floating-point numbers. For X = ⟨Xm , Xr ⟩ ∈ IFn×m andY = ⟨Ym , Yr ⟩ ∈ IFm×n , we define the following interval arithmetic:
X · Y = [▽(Xm · Ym − C), △(Xm · Ym + C)]
where
C = △((|Xm | + Xr ) · Yr + Xr · |Ym |).
Exercise 6.3. Write the Algorithm 5 and run it. Here, you need to use thefollowing command:$ gfortran matint.f90 ∼ /lib/slib.a ∼ /lib/lapack.a ∼ /lib/oblas.a
You need to make the matint.f90 in slib directory.

33

Algorithm 5 matint.f90
program matint
use slibreal(8), dimension(:,:), allocatable :: Aup, Adown, Bup, Bdown, Cup, Cdown
real(8), dimension(:,:), allocatable :: Amid, Arad, Bmid, Brad, C
integer :: n=5
call setround(0)call rand(Aup,n)call rand(Bup,n)Adown = Aup
Bdown = Bup!Addcall setround(1)Cup = Aup + Bup
call setround(-1)
Cdown = Adown + Bdown
call setround(0)!Multiplication
call setround(1)Amid = (Aup - Adown)/2d0 + Adown
Arad = Amid - Adowm
Bmid = (Bup - Bdown)/2d0 + Bdown
Brad = Bmid - Bdowm
C = ((abs(Amid) + Arad) .gemm. Brad) + (Arad .gemm. abs(Bmid))
Cup = (Amid .gemm . Bmid) + C
call setround(-1)
Cdown = (Amid .gemm . Bmid) - C
call setround(0)end program matint

6.2. Programming of Verified Theory for System of Linear
Equations34

Algorithm 6 Vlinear1.f90
program vlinear1
use slib
implicit nonereal(8), allocatable, dimension(:,:) :: A, b, x, R, I, Gd, Gu
real(8), allocatable, dimension(:,:) :: rd, ru, Gnorm, Rnorm, resnorm, D, Ainv, error
integer :: ncall setround(0)n = 1000
call rand(A,n)
call rand(b,n,1)
call inv(A,R)
x = R .gemm. bcall eye(I,n)
call setround(-1)
Gd = abs((R.gemm.A) - I)
rd = abs((A .gemm. x) - b)call setround(1)Gu = abs((R.gemm.A) - I)
ru = abs((A .gemm. x) - b)call output(ru)Gu = max(Gu,Gd)ru = max(ru,rd)call norm(Gu,0,Gnorm)
call norm(R,0,Rnorm)
call norm(ru,0,resnorm)
call setround(-1)
D = 1d0 - Gnorm
if (D(1,1) > 0d0) thencall setround(1)Ainv = Rnorm/D
error = Ainv*resnorm
write(*,*) error
else
write(*,*) ’Verification is failed...’end ifcall setround(0)end program vlinear135