pyplot (generic function with 1 method)

xxxxxxxxxx
 
begin 
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add("PyPlot")
    Pkg.add("PlutoUI")
    Pkg.add("IterativeSolvers")
    Pkg.add("IncompleteLU")
    Pkg.add("DataFrames")
    
    using IterativeSolvers
    using IncompleteLU
    using PlutoUI
    using PyPlot
    using DataFrames
    
    using LinearAlgebra
    using SparseArrays
    
    function pyplot(f;width=3,height=3)
      clf()
      f()
      fig=gcf()
      fig.set_size_inches(width,height)
      fig
    end
end

9.5 s

Practical iterative methods

3.8 μs

Incomplete LU (ILU) preconditioning

4.8 μs

Idea (Varga, Buleev, $\approx 1960$ : derive a preconditioner not from an additive decomposition but from the LU factorization.

LU factorization has large fill-in. For a preconditioner, just limit the fill-in to a fixed pattern.
Apply the standard LU factorization method, but calculate only those
Result: incomplete LU factors $L$ , $U$ , remainder $R$ :

$A = L U - R$

What about zero pivots which prevent such an algoritm from being computable ?

Theorem (Saad, Th. 10.2): If $A$ is an M-Matrix, then the algorithm to compute the incomplete LU factorization with a given pattern is stable. Moreover, $A = L U - R = M - N$ where $M = L U$ and $N = R$ is a regular splitting.

2.1 ms

Discussion

Generally better convergence properties than Jacobi, Gauss-Seidel
Block variants are possible
ILU Variants:
- ILUM: ("modified"): add ignored off-diagonal entries to main diagonal
- ILUT: ("threshold"): zero pattern calculated dynamically based on drop tolerance
- ILU0: Drop all fill-in
- Incomplete Cholesky: symmetric variant of ILU
Dependence on ordering
Can be parallelized using graph coloring
Not much theory: experiment for particular systems and see if it works well
I recommend it as the default initial guess for a sensible preconditioner

13.8 μs

Further approaches to preconditioning

These are based on ideas which are best explained and developed with multidimensional PDEs in mind.

Multigrid: gives indeed $O (N)$ optimal solver complexity in some situations. This is the holy grail method... I will try to discuss this later in the course.
Domain decomposition - based on the idea the subdivision of the computational domain into a number of subdomains and subsequent repeated solution of the smaller subdomain problems

8.8 μs

Iterative methods in Julia

Julia has some well maintained packages for iterative methods and preconditioning.

IterativeSolvers.jl: various Krylov subspace methods including conjugate gradients
IncompleteLU.jl: Incomplete LU factorizations
AlgebraicMultigrid.jl: Algebraic multigrid methods

7.0 μs

Random sparse M-Matrices

3.8 μs

We will test the methods with random sparse M matrices, so we define a function which gives us a random, strictly diagonally dominant M-Matrix which is not necessarily irreducible. For skew=0 it is also symmetric:

2.8 μs

sprandm (generic function with 1 method)

xxxxxxxxxx
 
function sprandm(n;p=0.5,skew=0)
    A=sprand(n,n,p) # random sparse matrix with positive entries
    for i=1:n       # set diagonal to zero
      A[i,i]=0
    end 
    A=A+(1.0-skew)*transpose(A) # symmetrize if necessary
    d=0.001*rand(n) # define a positive random diagonal vector
    for i=1:n # update to dominance
       d[i]+=sum(A[:,i]) 
    end
    Diagonal(d)-A  # create final matrix
end

55.3 μs

Test the method a bit...

4.5 μs

xxxxxxxxxx
 
N=5

80.0 ns

xxxxxxxxxx
 
A=sprandm(N,p=0.6,skew=1);

859 ms

	x1	x2	x3	x4	x5
	Float64	Float64	Float64	Float64	Float64
1	2.18746	-0.814626	-0.452134	-0.704067	0.0
2	-0.494782	2.59544	-0.183471	0.0	-0.288585
3	0.0	-0.795573	1.17374	-0.228419	-0.0102942
4	-0.993348	-0.985134	0.0	0.933306	0.0
5	-0.698695	0.0	-0.537345	0.0	0.299452

xxxxxxxxxx
 
DataFrame(A)

268 ms

Up to rounding errors, the inverse is nonnegative, as predicted by the theory. There are zero entries because it is not necessarily irreducible. Invertibility is guaranteed by strict diagonal dominance.

2.3 μs

Ainv

5×5 Array{Float64,2}:
 199.183  198.919  198.709  198.892  198.531
 134.768  135.099  134.757  134.647  134.829
 166.991  167.226  167.737  167.027  166.924
 354.249  354.317  353.733  354.883  353.62
 764.396  764.203  764.629  763.781  766.096

xxxxxxxxxx
 
Ainv=inv(Matrix(A))

323 ms

134.64717151847654

xxxxxxxxxx
 
minimum(Ainv)

418 ns

ρ_jacobi (generic function with 1 method)

xxxxxxxxxx
 
function  ρ_jacobi(A)
    B=I(size(A,1))-inv(Diagonal(A))*A;
    maximum(abs.(eigvals(Matrix(B))))
end

35.6 μs

0.9993509631121599

xxxxxxxxxx
 
ρ_jacobi(A)

7.4 ms

Preconditioners

Here, we define two preconditioners which are able to work together with IterativeSolvers.jl.

4.4 μs

Jacobi

2.5 μs

xxxxxxxxxx
 
begin   
    # Data struture: we store the inverse of the main diagonal
    struct JacobiPreconditioner
       invdiag::Vector
    end
    
    # Constructor:
    function JacobiPreconditioner(A::AbstractMatrix)
          n=size(A,1)
        invdiag=zeros(n)
        for i=1:n
            invdiag[i]=1.0/A[i,i]
        end
        JacobiPreconditioner(invdiag)
    end
    
    # Solution of preconditioning system   Mu=v
    # Method name and signature are compatible to IterativeSolvers.jl
    function  LinearAlgebra.ldiv!(u,precon::JacobiPreconditioner,v)
        invdiag=precon.invdiag
        n=length(invdiag)
        for i=1:n
            u[i]=invdiag[i]*v[i]
        end
    u
    end
    
    # In-place solution of preconditioning system
    function  LinearAlgebra.ldiv!(precon::JacobiPreconditioner,v)
    ldiv!(v,precon,v)
    end
    
end

1.6 ms

We can construct a the preconditioner then as follows:

4.9 μs

preconJacobi

JacobiPreconditionerinvdiagFloat641

0.457152

0.385291

0.851977

1.07146

3.33943

xxxxxxxxxx
 
preconJacobi=JacobiPreconditioner(A)

19.9 ms

Float641

0.457152

0.385291

0.851977

1.07146

3.33943

xxxxxxxxxx
 
ldiv!(preconJacobi,ones(N))

1.7 ms

ILU0

For this preconditioner, we need to store the matrix, the inverse of a modified diagonal and the indices of the main diagonal entries in the sparse matrix columns.

3.4 μs

xxxxxxxxxx
 
begin
    
    struct ILU0Preconditioner
    A::AbstractMatrix
    xdiag::Vector
    idiag::Vector
    end
    
    function ILU0Preconditioner(A::AbstractMatrix)
    n=size(A,1)
    colptr=A.colptr
        rowval=A.rowval
        nzval=A.nzval
        idiag=zeros(Int64,n)
        xdiag=zeros(n)
    
    
    # calculate main diagonal indices   
        for j=1:n
            for k=colptr[j]:colptr[j+1]-1
                i=rowval[k]
                if i==j
                    idiag[j]=k
                    break
                end
            end
        end
        
        # calculate modified diagonal 
        for j=1:n
            xdiag[j]=1/nzval[idiag[j]]
            for k=idiag[j]+1:colptr[j+1]-1
                i=rowval[k]
                for l=colptr[i]:colptr[i+1]-1
                    if rowval[l]==j
                        xdiag[i]-=nzval[l]*xdiag[j]*nzval[k]
                        break
                    end
                end
            end
        end
    
        ILU0Preconditioner(A,xdiag,idiag)
    end
    
    function  LinearAlgebra.ldiv!(u,precon::ILU0Preconditioner, v)
    A=precon.A
        colptr=A.colptr
        rowval=A.rowval
    n=size(A,1)
    nzval=A.nzval
        xdiag=precon.xdiag
        idiag=precon.idiag
        T=eltype(v)
    
    # forward substitution  
        for j=1:n
            x=zero(T)
            for k=colptr[j]:idiag[j]-1
                x+=nzval[k]*u[rowval[k]]
            end
            u[j]=xdiag[j]*(v[j]-x)
        end
        
    # backward substitution 
        for j=n:-1:1
        x=zero(T)
        for k=idiag[j]+1:colptr[j+1]-1
            x+=u[rowval[k]]*nzval[k]
        end
            u[j]-=x*xdiag[j]
        end
    u
    end
    
    function  LinearAlgebra.ldiv!(precon::ILU0Preconditioner,v)
    ldiv!(v,precon,v)
    end
​
    SparseArrays.nnz(precon::ILU0Preconditioner)=nnz(precon.A)
​
end

1.4 ms

xxxxxxxxxx
 
preconILU0=ILU0Preconditioner(A);

37.5 ms

Float641

2.93032

2.05551

2.89032

1.68817

3.88721

xxxxxxxxxx
 
ldiv!(preconILU0,ones(N))

2.9 ms

Simple iteration method with interface similar to IterativeSolvers.jl

2.2 μs

simple (generic function with 1 method)

xxxxxxxxxx
 
begin
    function simple!(u,A,b;tol=1.0e-10,log=true,maxiter=100,Pl=nothing)
    res=A*u-b # initial residual
    r0=norm(res) # residual norm
    history=[r0]  # intialize history recording
    for i=1:maxiter 
        u=u-ldiv!(Pl,res) # solve preconditioning system and update solution
        res=A*u-b         # calculate residual 
        r=norm(res)       # residual norm
        push!(history,r)  # record in history 
        if (r/r0)<tol     # check for relative tolerance
        return u,Dict( :resnorm => history)
        end
    end
    return u,Dict( :resnorm =>history )
    end
​
    simple(A,b;tol=1.0e-10, log=true,maxiter=100,Pl=nothing)=simple!(zeros(length(b)),A,b,tol=tol,maxiter=maxiter,log=log,Pl=Pl)
end

85.4 μs

Iterative Method comparison: symmetric problems

2.6 μs

xxxxxxxxxx
 
N1=100

62.0 ns

tol

1.0e-10

xxxxxxxxxx
 
tol=1.0e-10

67.0 ns

xxxxxxxxxx
 
A1=sprandm(N1,p=0.1,skew=0);

68.8 μs

xxxxxxxxxx
 
pyplot() do
    spy(A1)
end

194 ms

xxxxxxxxxx
 
A1Jacobi=JacobiPreconditioner(A1);

9.4 μs

xxxxxxxxxx
 
A1ILU0=ILU0Preconditioner(A1);

32.2 μs

Create also ILU preconditioners from IncompleteLU.jl: These have drop tolerance τ as parameter. The larger τ, the more entries of the LU factors are ignored.

6.4 μs

xxxxxxxxxx
 
A1ILUT_1=IncompleteLU.ilu(A1,τ=0.15);

133 ms

xxxxxxxxxx
 
A1ILUT_2=IncompleteLU.ilu(A1,τ=0.05);

281 μs

xxxxxxxxxx
 
nnz(A1ILU0)

2.6 ms

xxxxxxxxxx
 
nnz(A1ILUT_1)

335 ns

xxxxxxxxxx
 
nnz(A1ILUT_2)

451 ns

xxxxxxxxxx
 
nnz(lu(A1))

604 μs

Create a right hand side for testing

2.7 μs

Float641

0.000965422

0.000284515

9.91317e-5

0.000509845

0.000271775

0.000599232

0.000782524

0.00041395

0.000307366

0.000822928

0.000179343

0.000802294

0.000674138

0.00068135

0.000286788

0.000259899

0.000616573

0.000568139

0.000795105

0.000487304

more91

0.000903778

0.000796121

0.000291694

0.000771596

0.000862553

0.000274352

0.000145037

0.000968574

0.000663811

100

2.6471e-5

xxxxxxxxxx
 
b1=A1*ones(N1)

3.4 μs

So let us run this with Jacobi preconditioner. Theory tells it should converge...

2.1 μs

1Float641

0.00508829

0.0049741

0.00498199

0.00499919

0.00498366

0.00502427

0.00501246

0.00500443

0.00500781

0.00502883

0.00499319

0.00502764

0.00503714

0.00508468

0.00497958

0.00498624

0.00501386

0.00504387

0.00505017

0.0050265

more91

0.00503953

0.00502889

0.00498784

0.00502174

0.00503388

0.00496949

0.00497807

0.00503148

0.0050145

100

0.00497064

2Dict

:resnorm

Float641

0.00591899

0.0054469

0.00540714

0.0054061

0.00540533

0.00540513

0.00540486

0.0054046

0.00540433

0.00540406

0.00540379

0.00540352

0.00540325

0.00540298

0.00540271

0.00540243

0.00540216

0.00540189

0.00540162

0.00540135

more92

0.00538185

0.00538158

0.00538131

0.00538104

0.00538077

0.0053805

0.00538023

0.00537996

100

0.00537968

101

0.00537941

xxxxxxxxxx
 
sol_simple_jacobi,hist_simple_jacobi=simple(A1,b1,tol=tol,maxiter=100,log=true,Pl=A1Jacobi)

205 ms

After 100 steps we are far from the solution, and we need lots of steps to converge, so let us have a look at the spectral radius of the iteration matrix and compare it with the residual reduction in the last iteration step:

2.3 μs

0.99995

0.99995

xxxxxxxxxx
 
 ρ_jacobi(A1),(hist_simple_jacobi[:resnorm][end]/hist_simple_jacobi[:resnorm][end-1])

28.0 ms

It seams we have found a simple spetral radius estimator here ...

2.2 μs

Now for the ILU0 preconditioner:

2.2 μs

1Float641

0.0146098

0.0144986

0.0144998

0.0145256

0.0145013

0.0145409

0.0145176

0.0145135

0.014539

0.0145419

0.0145172

0.0145434

0.0145495

0.0145872

0.0144726

0.0144868

0.0145273

0.0145542

0.0145626

0.0145403

more91

0.0145041

0.0144763

0.0144468

0.0144723

0.0144816

0.0144176

0.0144324

0.0144882

0.014459

100

0.0144285

2Dict

:resnorm

Float641

0.00591899

0.00617601

0.0062019

0.00620316

0.00620246

0.00620158

0.00620068

0.00619978

0.00619887

0.00619797

0.00619706

0.00619616

0.00619526

0.00619435

0.00619345

0.00619255

0.00619164

0.00619074

0.00618984

0.00618894

more92

0.00612429

0.0061234

0.0061225

0.00612161

0.00612072

0.00611982

0.00611893

0.00611804

100

0.00611715

101

0.00611626

xxxxxxxxxx
 
sol_simple_ilu0,hist_simple_ilu0=simple(A1,b1,tol=tol,maxiter=100,log=true,Pl=A1ILU0)

52.2 ms

... the spectral radius estimate is a little bit better...

2.4 μs

0.9998541700844135

xxxxxxxxxx
 
hist_simple_ilu0[:resnorm][end]/hist_simple_ilu0[:resnorm][end-1]

8.3 μs

We have symmetric matrices, so let us try CG:

2.1 μs

Without preconditioning:

2.1 μs

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more91

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

100

1.0

Converged after 29 iterations.

xxxxxxxxxx
 
sol_cg,hist_cg=cg(A1,b1, reltol=tol,log=true,maxiter=100)

96.3 ms

With Jacobi preconditioning:

2.2 μs

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more91

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

100

1.0

Converged after 23 iterations.

xxxxxxxxxx
 
sol_cg_jacobi,hist_cg_jacobi=cg(A1,b1, reltol=tol,log=true,maxiter=100,Pl=A1Jacobi)

78.7 ms

With various variants of ILU preconditioners:

2.2 μs

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more91

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

100

1.0

Converged after 11 iterations.

xxxxxxxxxx
 
sol_cg_ilu0,hist_cg_ilu0=cg(A1,b1, reltol=tol,log=true,maxiter=100,Pl=A1ILU0)

92.5 ms

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more91

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

100

1.0

Converged after 11 iterations.

xxxxxxxxxx
 
sol_cg_ilut_1,hist_cg_ilut_1=cg(A1,b1, reltol=tol,log=true,maxiter=100,Pl=A1ILUT_1)

82.0 ms

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more91

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

100

1.0

Converged after 8 iterations.

xxxxxxxxxx
 
sol_cg_ilut_2,hist_cg_ilut_2=cg(A1,b1, reltol=tol,log=true,maxiter=100,Pl=A1ILUT_2)

111 μs

As we see, all CG variants converge within the given number of iterations steps.
Precoditioning helps
The better the preconditioner, the faster the iteration (though this also depends on the initial value)
The behaviour of the CG residual is not monotone

6.2 μs

xxxxxxxxxx
 
pyplot(width=10) do
    semilogy(hist_simple_ilu0[:resnorm],label="simple jacobi")
    semilogy(hist_simple_jacobi[:resnorm],label="simple ilu0")
    semilogy(hist_cg[:resnorm],label="cg")
    semilogy(hist_cg_jacobi[:resnorm],label="cg jacobi")
    semilogy(hist_cg_ilu0[:resnorm],label="cg ilu0")
    semilogy(hist_cg_ilut_1[:resnorm],label="cg ilut_1")
    semilogy(hist_cg_ilut_2[:resnorm],label="cg ilut_2")
    xlim(0,50)
    xlabel("iteration number k")
    ylabel("residual norm")
    legend(loc="lower right")
    grid()
end

131 ms

Nonsymmetric problems

Here, we skip the simple iteration and look at the performance of some Krylov subspace methods.

4.0 μs

xxxxxxxxxx
 
N2=1000

67.0 ns

xxxxxxxxxx
 
A2=sprandm(N2,p=0.1,skew=1);

4.8 ms

xxxxxxxxxx
 
b2=A2*ones(N2);

105 μs

xxxxxxxxxx
 
A2Jacobi=JacobiPreconditioner(A2);

66.3 μs

xxxxxxxxxx
 
A2ILU0=ILU0Preconditioner(A2);

3.7 ms

xxxxxxxxxx
 
A2ILUT=IncompleteLU.ilu(A2,τ=0.1);

12.1 ms

Try CG:

2.3 μs

1Float641

1.68069

1.483

1.50429

1.65334

1.70164

1.62845

1.50726

1.67786

1.65759

1.72771

1.75626

1.61094

1.70122

1.77005

1.60317

1.50787

1.58845

1.64992

1.49178

1.83337

more991

1.69881

992

1.55268

993

1.61124

994

1.56993

995

1.58721

996

1.5261

997

1.73719

998

1.4985

999

1.63193

1000

1.58343

Not converged after 100 iterations.

xxxxxxxxxx
 
sol2_cg,hist2_cg=cg(A2,b2, reltol=tol,log=true,maxiter=100)

6.5 ms

1Float641

632.834

446.586

462.752

595.862

642.002

574.132

460.705

619.735

599.501

666.189

698.139

554.976

637.467

701.963

546.974

463.757

531.491

595.723

456.943

770.731

more991

633.145

992

516.236

993

555.144

994

519.384

995

538.26

996

487.58

997

678.325

998

461.635

999

582.091

1000

540.247

Not converged after 100 iterations.

xxxxxxxxxx
 
sol2_cg_jacobi,hist2_cg_jacobi=cg(A2,b2, reltol=tol,log=true,maxiter=100,Pl=A2Jacobi)

27.9 ms

1Float641

-0.206584

0.140421

0.110268

-0.156276

-0.247687

-0.114002

0.103738

-0.204957

-0.16342

-0.297926

-0.342681

-0.0874407

-0.246283

-0.36456

-0.0691687

0.0951193

-0.0462991

-0.153633

0.131778

-0.479451

more991

-0.243986

992

0.0201848

993

-0.0976511

994

-0.0167221

995

-0.0406539

996

0.0671789

997

-0.317953

998

0.114745

999

-0.127749

1000

-0.0368361

Not converged after 100 iterations.

xxxxxxxxxx
 
sol2_cg_ILU0,hist2_cg_ILU0=cg(A2,b2, reltol=tol,log=true,maxiter=100,Pl=A2ILU0)

1.4 s

Use the bicgstabl method from IterativeSolvers.jl:

xxxxxxxxxx
 
md"""
Use the `bicgstabl` method from IterativeSolvers.jl:
"""

2.4 μs

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more991

1.0

992

1.0

993

1.0

994

1.0

995

1.0

996

1.0

997

1.0

998

1.0

999

1.0

1000

1.0

Converged after 6 iterations.

xxxxxxxxxx
 
sol2_bicgstab,hist2_bicgstab=bicgstabl(A2,b2,reltol=tol,log=true,max_mv_products=100)

82.0 ms

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more991

1.0

992

1.0

993

1.0

994

1.0

995

1.0

996

1.0

997

1.0

998

1.0

999

1.0

1000

1.0

Converged after 5 iterations.

xxxxxxxxxx
 
sol2_bicgstab_jacobi,hist2_bicgstab_jacobi=bicgstabl(A2,b2,reltol=tol,log=true,max_mv_products=100,Pl=A2Jacobi)

5.0 ms

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more991

1.0

992

1.0

993

1.0

994

1.0

995

1.0

996

1.0

997

1.0

998

1.0

999

1.0

1000

1.0

Converged after 6 iterations.

xxxxxxxxxx
 
sol2_bicgstab_ilu0,hist2_bicgstab_ilu0=bicgstabl(A2,b2,reltol=tol,log=true,max_mv_products=100,Pl=A2ILU0)

417 ms

1Float641

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

more991

1.0

992

1.0

993

1.0

994

1.0

995

1.0

996

1.0

997

1.0

998

1.0

999

1.0

1000

1.0

Converged after 3 iterations.

xxxxxxxxxx
 
sol2_bicgstab_ilut,hist2_bicgstab_ilut=bicgstabl(A2,b2,reltol=tol,log=true,max_mv_products=100,Pl=A2ILUT)

83.6 ms

CG does not converge - the case is also not covered by the theory
Various preconditioners improve the convergence
Is there a bug in the implementation of my ILU0 ?

5.3 μs

xxxxxxxxxx
 
pyplot(width=10) do
    semilogy(hist2_cg[:resnorm],label="cg")
    semilogy(hist2_cg_jacobi[:resnorm],label="cg jacobi")
    semilogy(hist_cg_ilu0[:resnorm],label="cg ilu0")
    semilogy(hist2_bicgstab[:resnorm],label="bicgstab")
    semilogy(hist2_bicgstab_jacobi[:resnorm],label="bicgstab jacobi")
    semilogy(hist2_bicgstab_ilu0[:resnorm],label="bicgstab ilu0")
    semilogy(hist2_bicgstab_ilut[:resnorm],label="bicgstab ilut")
    xlim(0,25)
    xlabel("iteration number k")
    ylabel("residual norm")
    legend(loc="lower right")
    grid()
end

71.9 ms