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