# Testing Numerical Algorithms in Julia with Matrix Depot

Weijian Zhang

Manchester Numerical Linear Algebra Group Talk

January 13, 2015

## 1. Features in Julia I Found Useful

A language that doesn't affect the way you think about programming is not worth knowing. -- Alan Jay Perlis

• Multiple dispatch

More details can be found at: http://nbviewer.ipython.org/gist/StefanKarpinski/b8fe9dbb36c1427b9f22

In [17]:
function f1(x, y)
x + y
end

Out[17]:
f1 (generic function with 1 method)

In [18]:
function f1(x::Int, y::Int)
x - y
end

Out[18]:
f1 (generic function with 2 methods)

In [19]:
function f1(x::Float64, y::Int)
x^y
end

Out[19]:
f1 (generic function with 3 methods)

In [20]:
f1(2, 3)

Out[20]:
-1

In [23]:
f1(2.0, 6)

Out[23]:
64.0

In [24]:
f1(2.0, 3.0)

Out[24]:
5.0

• Type system
In [8]:
# Type hierarchy
typeof(3)

Out[8]:
Int64

In [11]:
super(Int64)

Out[11]:
Signed

In [13]:
super(Signed)

Out[13]:
Integer

In [14]:
super(Integer)

Out[14]:
Real

In [15]:
super(Real)

Out[15]:
Number

In [16]:
super(Number)

Out[16]:
Any

In [1]:
# abstract type
abstract PolarAlg

# composite types
type NewtonAlg <: PolarAlg
maxiter::Int                # maximum number of iterations
scale::Bool                 # whether to scale
verbose::Bool               # whether to show procedural information
tol::Float64                # convergence tolerance
end

type HalleyAlg <: PolarAlg
maxiter::Int
verbose::Bool
tol::Float64
end

function solve!(alg::NewtonAlg, X::Matrix, U::Matrix, H::Matrix)
# some implementation
println("use Newton's method")
end

function solve!(alg::HalleyAlg, X::Matrix, U::Matrix, H::Matrix)
# some implementation
println("use Halley's method")
end

Out[1]:
solve! (generic function with 2 methods)

In [7]:
newton = NewtonAlg(100, false, false, 1.0e-6);

In [8]:
halley = HalleyAlg(100, false, 1.0e-6);

In [5]:
solve!(newton, eye(3), eye(3), eye(3))

use Newton's method


In [6]:
solve!(halley, eye(3), eye(3), eye(3))

use Halley's method


• Loops can be as fast as C
In [15]:
x = Float64[1:10000000];
y = Float64[20:10000019];

In [16]:
function f2(x, y)
r = exp(-abs(x-y))
return r
end

Out[16]:
f2 (generic function with 1 method)

In [22]:
@time f2(x,y);

elapsed time: 0.46564706 seconds (320000392 bytes allocated, 13.67% gc time)


In [19]:
r = similar(x)
function f2!(r, x, y)
for i = 1:length(x)
r[i] = exp(-abs(x[i] - y[i]))
end
end

Out[19]:
f2! (generic function with 1 method)

In [21]:
@time f2!(r, x, y)

elapsed time: 0.31581529 seconds (80 bytes allocated)


• Call low level library (BLAS, LAPACK, etc)
In [41]:
A = randn(4,4);
B = randn(4000,4000);

In [48]:
# on the first call the functions gets compiled
lu(A);

In [43]:
@time lu(B); # indicate time and the amount of memory allocation

elapsed time: 1.483885573 seconds (384032824 bytes allocated)


In [49]:
lufact(A);

In [45]:
@time lufact(B);

elapsed time: 1.423192229 seconds (128016432 bytes allocated)


In [50]:
LAPACK.getrf!(A);

In [47]:
@time LAPACK.getrf!(B);

elapsed time: 1.40171309 seconds (16288 bytes allocated)



## 2. Testing Numerical Algorithms

• Aim at speed and accuracy

• Error analysis (backward and forward)

• Numerical tests

We need tools that

• support reproducible research
• easy to access and use

I introduce you Matrix Depot, a test matrix collection for Julia.

To install Matrix Depot, simply type

In []:
Pkg.add("MatrixDepot")


Here are the basic usages.

In [52]:
using MatrixDepot

In [53]:
matrixdepot() # see all the matrices in the collection


| symmetric |  inverse  | ill-cond  |  pos-def  |  eigen    |
vand|           |     *     |     *     |           |           |
poisson|     *     |     *     |           |     *     |     *     |
frank|           |           |     *     |           |     *     |
minij|     *     |     *     |           |     *     |     *     |
clement|     *     |     *     |           |           |     *     |
tridiag|     *     |     *     |     *     |     *     |     *     |
circul|     *     |           |           |     *     |     *     |
dingdong|     *     |           |           |           |     *     |
hadamard|           |     *     |           |           |     *     |
moler|     *     |     *     |     *     |     *     |           |
invol|           |     *     |     *     |           |     *     |
fiedler|     *     |     *     |           |           |     *     |
binomial|           |           |           |           |           |
lehmer|     *     |     *     |           |     *     |           |
invhilb|     *     |     *     |     *     |     *     |           |
lotkin|           |     *     |     *     |           |     *     |
wilkinson|     *     |           |           |           |     *     |
triw|           |     *     |     *     |           |           |
rosser|           |           |     *     |           |     *     |
magic|           |     *     |           |           |           |
kahan|           |     *     |     *     |           |           |
pascal|     *     |     *     |     *     |     *     |     *     |
chebspec|           |           |           |           |     *     |
chow|           |           |           |           |     *     |
sampling|           |           |           |           |     *     |
hilb|     *     |     *     |     *     |     *     |           |
cauchy|     *     |     *     |     *     |     *     |           |
pei|     *     |     *     |     *     |     *     |           |
parter|           |           |           |           |     *     |
forsythe|           |     *     |     *     |           |     *     |
randcorr|     *     |           |           |           |           |
neumann|           |           |           |           |     *     |
grcar|           |           |           |           |     *     |


In [54]:
matrixdepot("hilb", 6) # generate a Hilbert matrix

Out[54]:
6x6 Array{Float64,2}:
1.0       0.5       0.333333  0.25      0.2       0.166667
0.5       0.333333  0.25      0.2       0.166667  0.142857
0.333333  0.25      0.2       0.166667  0.142857  0.125
0.25      0.2       0.166667  0.142857  0.125     0.111111
0.2       0.166667  0.142857  0.125     0.111111  0.1
0.166667  0.142857  0.125     0.111111  0.1       0.0909091

In [55]:
matrixdepot("circul", 7) # generate a circul matrix

Out[55]:
7x7 Array{Float64,2}:
1.0  2.0  3.0  4.0  5.0  6.0  7.0
7.0  1.0  2.0  3.0  4.0  5.0  6.0
6.0  7.0  1.0  2.0  3.0  4.0  5.0
5.0  6.0  7.0  1.0  2.0  3.0  4.0
4.0  5.0  6.0  7.0  1.0  2.0  3.0
3.0  4.0  5.0  6.0  7.0  1.0  2.0
2.0  3.0  4.0  5.0  6.0  7.0  1.0

In [56]:
matrixdepot("circul") # see input options and properties

Circul matrix:

Input options:

(type), vec, col_dim: a vector and the column dimension

(type), vec: a vector

(type), dim: the dimension of the matrix

['symmetric', 'pos-def', 'eigen']


In [57]:
matrixdepot("symmetric") # see all symmetric test matrices

Out[57]:
16-element Array{ASCIIString,1}:
"hilb"
"cauchy"
"circul"
"dingdong"
"invhilb"
"moler"
"pascal"
"pei"
"clement"
"fiedler"
"minij"
"tridiag"
"lehmer"
"randcorr"
"poisson"
"wilkinson"

In [58]:
matrixdepot("symmetric", "ill-cond") # symmetric and ill-conditioned

Out[58]:
7-element Array{ASCIIString,1}:
"hilb"
"cauchy"
"invhilb"
"moler"
"pascal"
"pei"
"tridiag"


Given a property, we can loop through all the matrices having this propery

In [59]:
A = eye(4)

print("Identity matrix")

for mat in matrixdepot("symmetric", "ill-cond", "inverse")
print(" x $mat matrix") A = A * full(matrixdepot(mat, 4)) end println(" =") A  Identity matrix x hilb matrix x cauchy matrix x invhilb matrix x moler matrix x pascal matrix x pei matrix x tridiag matrix =  Out[59]: 4x4 Array{Float64,2}: 153.12 -11.919 -15.4345 296.937 109.896 -8.91857 -11.5976 214.433 86.7524 -7.15714 -9.32857 169.702 71.9139 -5.98707 -7.81497 140.876  In [60]: using PolarFact # compute the polar decomposition  In [61]: for i in 1:5 polarfact(rand(8,8), alg = :halley, verbose =:true) # by Halley's method end  Iter. Rel. err. Obj. 1 6.056236e-01 2.516655e+00 2 2.467945e-01 1.409655e+00 3 1.690446e-01 1.035278e+00 4 2.006586e-01 1.441690e-01 5 2.486222e-02 1.009726e-04 6 1.698781e-05 3.009658e-14 7 4.886502e-15 1.163999e-15 Iter. Rel. err. Obj. 1 6.091343e-01 1.756239e+00 2 3.579344e-01 4.065556e-01 3 7.376543e-02 2.269951e-03 4 3.707004e-04 2.764532e-10 5 4.510789e-11 6.006480e-16 Iter. Rel. err. Obj. 1 6.462653e-01 2.466008e+00 2 2.941734e-01 9.094721e-01 3 2.535364e-01 1.317529e-01 4 3.244073e-02 9.924493e-05 5 2.376972e-05 3.583071e-14 6 8.702917e-15 9.506285e-16 Iter. Rel. err. Obj. 1 6.206390e-01 3.045682e+00 2 3.561950e-01 1.456327e+00 3 2.099683e-01 8.098815e-01 4 1.352440e-01 3.180147e-02 5 4.589004e-03 7.784455e-07 6 1.113809e-07 1.013079e-15 Iter. Rel. err. Obj. 1 5.897255e-01 2.969957e+00 2 3.451832e-01 1.257243e+00 3 2.075257e-01 8.133562e-01 4 1.893654e-01 7.336899e-02 5 1.489048e-02 1.627009e-05 6 3.252529e-06 1.165734e-15 7 2.931882e-16 9.298118e-16  Looks like the algorithm converges at about 5 to 7 iterations. But if we test on some ill-condtioned matrices, we will see something differet. In [63]: matrixdepot("ill-cond")  Out[63]: 15-element Array{ASCIIString,1}: "hilb" "cauchy" "frank" "invhilb" "forsythe" "triw" "moler" "pascal" "kahan" "pei" "vand" "invol" "lotkin" "tridiag" "rosser"  In [65]: A = Array(Float64, 8, 8) for mat in matrixdepot("ill-cond") copy!(A, matrixdepot(mat, 8)) println("For$mat matrix: ")
polarfact(A, alg = : halley, verbose = true)
end

For hilb matrix:
Iter.    Rel. err.        Obj.
1     5.601430e-01     1.794116e+00
2     3.704940e-01     1.777826e+00
3     4.347670e-01     1.769705e+00
4     4.723992e-01     1.763472e+00
5     2.005865e-01     1.706739e+00
6     4.278905e-01     1.701279e+00
7     2.070442e-01     1.714559e+00
8     2.098918e-01     1.655083e+00
9     4.508888e-01     1.581553e+00
10     2.091895e-01     1.633941e+00
11     1.333539e-01     1.603491e+00
12     3.495398e-01     1.527235e+00
13     3.655269e-01     1.517672e+00
14     4.715770e-02     1.523321e+00
15     1.338632e-01     1.453862e+00
16     3.374583e-01     1.372003e+00
17     2.881877e-01     1.384344e+00
18     2.635389e-02     1.382941e+00
19     7.967482e-02     1.362638e+00
20     2.312930e-01     1.196586e+00
21     4.750646e-01     4.495227e-01
22     2.239057e-01     5.165468e-03
23     2.582698e-03     4.512566e-09
24     2.256283e-09     2.220446e-16
For cauchy matrix:
Iter.    Rel. err.        Obj.
1     3.194190e-01     1.768889e+00
2     5.714797e-01     1.766394e+00
3     2.107487e-01     1.777069e+00
4     3.548663e-01     1.771669e+00
5     3.752544e-01     1.757173e+00
6     1.973082e-01     1.725540e+00
7     4.647824e-01     1.702592e+00
8     2.852555e-01     1.744404e+00
9     1.655482e-01     1.706240e+00
10     3.980794e-01     1.598188e+00
11     3.096433e-01     1.569067e+00
12     8.507988e-02     1.561310e+00
13     2.396193e-01     1.513854e+00
14     4.172527e-01     1.465708e+00
15     8.921506e-02     1.469800e+00
16     8.903575e-02     1.440458e+00
17     2.504474e-01     1.343648e+00
18     4.069905e-01     1.372590e+00
19     8.995150e-02     1.381922e+00
20     4.576589e-02     1.375042e+00
21     1.367586e-01     1.314849e+00
22     3.612055e-01     9.028514e-01
23     4.290455e-01     9.236526e-02
24     4.614719e-02     2.854394e-05
25     1.427197e-05     8.356881e-16
26     6.398887e-16     4.440892e-16
For frank matrix:
Iter.    Rel. err.        Obj.
1     6.649384e-01     7.603741e+01
2     6.517834e-01     8.141827e+00
3     5.549318e-01     1.330965e+00
4     1.493358e-01     1.192696e+00
5     9.452366e-03     1.191904e+00
6     2.594412e-02     1.187750e+00
7     7.666389e-02     1.151214e+00
8     2.021877e-01     8.808839e-01
9     2.630098e-01     1.510782e-01
10     3.605221e-02     1.850754e-04
11     4.146880e-05     2.791328e-13
12     6.247624e-14     8.844786e-16
For invhilb matrix:
Iter.    Rel. err.        Obj.
1     1.779757e+02     6.715041e+24
2     6.658331e-01     7.480762e+23
3     6.662219e-01     8.326144e+22
4     6.666975e-01     9.249286e+21
5     6.665926e-01     1.027878e+21
6     6.667387e-01     1.141867e+20
7     6.666314e-01     1.268913e+19
8     6.668613e-01     1.409595e+18
9     6.667674e-01     1.565816e+17
10     6.667779e-01     1.739346e+16
11     6.666957e-01     1.932521e+15
12     6.666002e-01     2.147648e+14
13     6.665922e-01     2.386891e+13
14     6.666757e-01     2.652006e+12
15     6.666736e-01     2.946609e+11
16     6.666680e-01     3.274002e+10
17     6.666667e-01     3.637780e+09
18     6.666667e-01     4.041977e+08
19     6.666667e-01     4.491086e+07
20     6.666663e-01     4.990095e+06
21     6.666614e-01     5.544547e+05
22     6.666090e-01     6.160574e+04
23     6.662152e-01     6.844715e+03
24     6.641829e-01     7.601159e+02
25     6.442385e-01     8.404581e+01
26     5.789385e-01     8.938703e+00
27     3.956075e-01     6.806880e-01
28     8.797931e-02     5.434606e-03
29     7.976780e-04     5.132634e-09
30     7.541907e-10     1.138304e-15
For forsythe matrix:
Iter.    Rel. err.        Obj.
1     2.980232e-08     1.000000e+00
For triw matrix:
Iter.    Rel. err.        Obj.
1     5.639175e-01     2.271852e+00
2     2.568783e-01     1.500268e+00
3     1.133687e-01     1.358870e+00
4     2.546539e-01     6.797687e-01
5     1.492881e-01     2.004530e-02
6     3.942868e-03     2.284541e-07
7     4.478627e-08     8.923526e-16
For moler matrix:
Iter.    Rel. err.        Obj.
1     6.663801e-01     5.403350e+01
2     6.616867e-01     5.682165e+00
3     5.663585e-01     1.575035e+00
4     1.266461e-01     1.499537e+00
5     2.440978e-02     1.497795e+00
6     7.307789e-02     1.484546e+00
7     2.150031e-01     1.372038e+00
8     5.036699e-01     7.294252e-01
9     3.719765e-01     2.674173e-02
10     1.338580e-02     5.461016e-07
11     2.730508e-07     4.440902e-16
For pascal matrix:
Iter.    Rel. err.        Obj.
1     6.666680e-01     3.293208e+06
2     6.666692e-01     3.659117e+05
3     6.665939e-01     4.065659e+04
4     6.665099e-01     4.517148e+03
5     6.679726e-01     5.015910e+02
6     6.744717e-01     5.532068e+01
7     6.770768e-01     5.742305e+00
8     5.785485e-01     3.588023e-01
9     1.444438e-01     9.926283e-04
10     4.959823e-04     2.960008e-11
11     1.480009e-11     4.440894e-16
For kahan matrix:
Iter.    Rel. err.        Obj.
1     3.006378e-01     1.444153e+00
2     3.229069e-01     5.473380e-01
3     1.357121e-01     6.869008e-03
4     1.576762e-03     7.654172e-09
5     1.755293e-09     6.336077e-16
For pei matrix:
Iter.    Rel. err.        Obj.
1     6.557377e-01     8.599839e+00
2     5.771798e-01     7.162300e-01
3     2.329699e-01     9.718345e-03
4     4.824011e-03     5.653911e-08
5     2.826955e-08     6.661338e-16
For vand matrix:
Iter.    Rel. err.        Obj.
1     6.666683e-01     6.577410e+11
2     6.666666e-01     7.308234e+10
3     6.666667e-01     8.120260e+09
4     6.666669e-01     9.022511e+08
5     6.666669e-01     1.002501e+08
6     6.666632e-01     1.113890e+07
7     6.666469e-01     1.237655e+06
8     6.666280e-01     1.375170e+05
9     6.667481e-01     1.527937e+04
10     6.678120e-01     1.697383e+03
11     6.705761e-01     1.882644e+02
12     6.745712e-01     2.058795e+01
13     6.542700e-01     1.986638e+00
14     2.958021e-01     7.148929e-02
15     1.498902e-02     1.634970e-05
16     3.484240e-06     8.133583e-16
17     3.227302e-16     5.240474e-16
For invol matrix:
Iter.    Rel. err.        Obj.
1     6.666669e-01     3.095345e+10
2     6.666669e-01     3.439272e+09
3     6.666674e-01     3.821413e+08
4     6.666682e-01     4.246015e+07
5     6.666652e-01     4.717794e+06
6     6.666445e-01     5.241990e+05
7     6.665402e-01     5.824402e+04
8     6.658467e-01     6.471209e+03
9     6.626429e-01     7.185989e+02
10     6.528649e-01     7.940876e+01
11     5.924422e-01     8.401540e+00
12     3.810211e-01     6.127615e-01
13     7.759668e-02     3.837432e-03
14     5.726324e-04     1.614407e-09
15     2.411788e-10     1.126378e-15
For lotkin matrix:
Iter.    Rel. err.        Obj.
1     5.692227e-01     2.249147e+00
2     2.200842e-01     1.747537e+00
3     2.353150e-01     1.748748e+00
4     1.456131e-01     1.782931e+00
5     1.512761e-01     1.709302e+00
6     2.532741e-01     1.721104e+00
7     8.137599e-02     1.753000e+00
8     1.346685e-01     1.694361e+00
9     2.312609e-01     1.625922e+00
10     7.436181e-02     1.620114e+00
11     1.027403e-01     1.588949e+00
12     2.393812e-01     1.447061e+00
13     1.841833e-01     1.391504e+00
14     3.699398e-02     1.382283e+00
15     1.059511e-01     1.343347e+00
16     2.348243e-01     1.331425e+00
17     1.512426e-01     1.364680e+00
18     1.741847e-02     1.363073e+00
19     5.131329e-02     1.336792e+00
20     1.418702e-01     1.128019e+00
21     2.374547e-01     3.312596e-01
22     6.763881e-02     1.817740e-03
23     3.490114e-04     2.014573e-10
24     3.866758e-11     9.480264e-16
For tridiag matrix:
Iter.    Rel. err.        Obj.
1     6.826864e-01     1.609515e+00
2     3.786431e-01     4.528488e-01
3     2.115685e-01     7.114597e-03
4     3.554137e-03     1.473770e-08
5     7.368851e-09     2.220446e-16
For rosser matrix:
Iter.    Rel. err.        Obj.
1     6.666670e-01     1.259239e+07
2     6.666674e-01     1.399155e+06
3     6.666681e-01     1.554617e+05
4     6.666646e-01     1.727360e+04
5     6.666028e-01     1.919366e+03
6     6.659563e-01     2.133400e+02
7     6.599113e-01     2.378807e+01
8     6.081898e-01     2.779588e+00
9     3.245594e-01     1.223730e+00
10     1.460331e-01     1.111469e+00
11     3.051364e-01     5.481643e-01
12     1.619601e-01     1.537983e-02
13     3.889409e-03     1.533311e-07
14     3.859627e-08     8.634586e-16



The lesson is some algorithms perform badly for ill condtioned matrices.

A scaled Newton's method is much more robust.

In [66]:
A = Array(Float64, 8, 8)
for mat in matrixdepot("ill-cond")
copy!(A, matrixdepot(mat, 8))
println("For \$mat matrix: ")
polarfact(A, alg = :newton, verbose = true)
end

For hilb matrix:
Iter.    Rel. err.        Obj.
1     3.935096e+04     6.937286e+09
2     9.990390e-01     1.405558e+04
3     9.593756e-01     3.960797e+01
4     8.191273e-01     1.155700e+00
5     3.154410e-01     3.468750e-02
6     1.694974e-02     4.958326e-05
7     2.479056e-05     4.499344e-10
8     2.249672e-10     2.304525e-20
For cauchy matrix:
Iter.    Rel. err.        Obj.
1     1.129361e+05     2.528751e+10
2     9.992403e-01     3.713719e+04
3     9.758229e-01     5.828685e+01
4     8.396564e-01     1.585939e+00
5     3.877349e-01     1.756765e-02
6     8.679102e-03     1.177750e-05
7     5.888695e-06     2.358062e-11
8     1.179031e-11     8.206269e-23
For frank matrix:
Iter.    Rel. err.        Obj.
1     1.720690e+01     1.055928e+05
2     9.912604e-01     1.486070e+01
3     6.137595e-01     6.185575e-01
4     1.349603e-01     2.313439e-02
5     5.580190e-03     7.616820e-05
6     1.882930e-05     1.004251e-09
7     2.554067e-10     7.377996e-16
For invhilb matrix:
Iter.    Rel. err.        Obj.
1     9.999928e-01     6.937286e+09
2     9.990390e-01     1.405558e+04
3     9.593756e-01     3.960797e+01
4     8.191273e-01     1.155700e+00
5     3.154410e-01     3.468750e-02
6     1.694974e-02     4.958326e-05
7     2.479056e-05     4.499344e-10
8     2.249672e-10     2.304522e-20
For forsythe matrix:
Iter.    Rel. err.        Obj.
1     4.096000e+03     1.677722e+07
2     9.997559e-01     0.000000e+00
3     0.000000e+00     0.000000e+00
For triw matrix:
Iter.    Rel. err.        Obj.
1     2.125000e+00     1.596719e+02
2     8.886894e-01     5.957393e-01
3     1.400441e-01     1.030830e-02
4     2.399435e-03     1.336599e-05
5     2.920694e-06     3.069977e-11
6     6.275015e-12     6.200281e-16
For moler matrix:
Iter.    Rel. err.        Obj.
1     1.014429e+01     5.575228e+04
2     9.922752e-01     2.681202e+00
3     4.850401e-01     9.385416e-02
4     4.396209e-02     4.301413e-04
5     2.149988e-04     1.289403e-08
6     6.447016e-09     2.300304e-17
For pascal matrix:
Iter.    Rel. err.        Obj.
1     5.111182e-01     7.409759e+06
2     9.901372e-01     7.382039e+02
3     9.376860e-01     8.501632e+00
4     6.893363e-01     2.415435e-01
5     1.026379e-01     5.754747e-04
6     2.876179e-04     2.262175e-08
7     1.131088e-08     7.019621e-17
For kahan matrix:
Iter.    Rel. err.        Obj.
1     9.808429e-01     3.929094e+00
2     4.837591e-01     1.586999e-01
3     4.325619e-02     1.495357e-03
4     3.987809e-04     2.093390e-07
5     4.956653e-08     1.021058e-14
For pei matrix:
Iter.    Rel. err.        Obj.
1     7.704899e-01     3.266667e+00
2     4.917670e-01     1.020833e-01
3     4.674139e-02     1.465311e-03
4     7.316569e-04     3.890257e-07
5     1.945128e-07     3.842100e-14
For vand matrix:
Iter.    Rel. err.        Obj.
1     9.930260e-01     2.937385e+08
2     9.975664e-01     2.852154e+03
3     9.304971e-01     1.474986e+01
4     6.661113e-01     8.668548e-01
5     1.808262e-01     4.743964e-02
6     1.251493e-02     3.959454e-04
7     1.074377e-04     3.194464e-08
8     8.669130e-09     6.323389e-16
For invol matrix:
Iter.    Rel. err.        Obj.
1     9.730891e-01     6.964539e+10
2     9.993645e-01     4.433996e+04
3     9.728310e-01     4.130526e+01
4     7.795700e-01     7.087216e-01
5     1.485346e-01     1.242074e-02
6     2.736360e-03     6.594547e-06
7     1.685634e-06     9.162936e-12
8     2.053508e-12     5.123939e-16
For lotkin matrix:
Iter.    Rel. err.        Obj.
1     2.382970e+04     1.125293e+10
2     9.989007e-01     2.387047e+04
3     9.557838e-01     4.912136e+01
4     7.925350e-01     1.801859e+00
5     2.363931e-01     1.058813e-01
6     1.924075e-02     8.674923e-04
7     1.608014e-04     7.964444e-08
8     1.481326e-08     7.042977e-16
For tridiag matrix:
Iter.    Rel. err.        Obj.
1     7.905694e-01     1.246667e+01
2     7.405554e-01     4.641619e-01
3     1.744364e-01     8.118547e-03
4     4.037445e-03     2.925963e-06
5     1.462978e-06     1.704412e-12
6     8.522058e-13     4.589814e-25
For rosser matrix:
Iter.    Rel. err.        Obj.
1     1.633208e+04     8.018671e+15
2     9.999995e-01     3.360144e+03
3     9.791008e-01     3.776066e-02
4     1.395747e-02     2.148643e-04
5     8.063111e-05     9.235922e-09
6     3.466220e-09     2.968025e-16