juliasmoothoptimizers / krylov.jl

A Julia Basket of Hand-Picked Krylov Methods

Krylov.jl: A Julia basket of hand-picked Krylov methods

Documentation	CI	Coverage	DOI	Downloads

How to Cite

If you use Krylov.jl in your work, please cite it using the metadata given in CITATION.cff.

BibTeX

@article{montoison-orban-2023,
  author  = {Montoison, Alexis and Orban, Dominique},
  title   = {{Krylov.jl: A Julia basket of hand-picked Krylov methods}},
  journal = {Journal of Open Source Software},
  volume  = {8},
  number  = {89},
  pages   = {5187},
  year    = {2023},
  doi     = {10.21105/joss.05187}
}

Content

This package provides implementations of certain of the most useful Krylov method for a variety of problems:

1. Square or rectangular full-rank systems

Ax = b

should be solved when b lies in the range space of A. This situation occurs when

• A is square and nonsingular,
• A is tall and has full column rank and b lies in the range of A.

2. Linear least-squares problems

minimize ‖b - Ax‖

should be solved when b is not in the range of A (inconsistent systems), regardless of the shape and rank of A. This situation mainly occurs when

• A is square and singular,
• A is tall and thin.

Underdetermined systems are less common but also occur.

If there are infinitely many such x (because A is column rank-deficient), one with minimum norm is identified

minimize ‖x‖   subject to   x ∈ argmin ‖b - Ax‖.

3. Linear least-norm problems

minimize ‖x‖   subject to   Ax = b

should be solved when A is column rank-deficient but b is in the range of A (consistent systems), regardless of the shape of A. This situation mainly occurs when

• A is square and singular,
• A is short and wide.

Overdetermined systems are less common but also occur.

4. Adjoint systems

Ax = b   and   Aᴴy = c

where A can have any shape.

5. Saddle-point and Hermitian quasi-definite systems

[M     A]  [x] = [b]
[Aᴴ   -N]  [y]    [c]

where A can have any shape.

6. Generalized saddle-point and non-Hermitian partitioned systems

[M   A]  [x] = [b]
[B   N]  [y]    [c]

where A can have any shape and B has the shape of Aᴴ. A, B, b and c must be all nonzero.

Krylov solvers are particularly appropriate in situations where such problems must be solved but a factorization is not possible, either because:

• A is not available explicitly,
• A would be dense or would consume an excessive amount of memory if it were materialized,
• factors would consume an excessive amount of memory.

Iterative methods are recommended in either of the following situations:

• the problem is sufficiently large that a factorization is not feasible or would be slow,
• an effective preconditioner is known in cases where the problem has unfavorable spectral structure,
• the operator can be represented efficiently as a sparse matrix,
• the operator is fast, i.e., can be applied with better complexity than if it were materialized as a matrix. Certain fast operators would materialize as dense matrices.

Features

All solvers in Krylov.jl have in-place version, are compatible with GPU and work in any floating-point data type.

How to Install

Krylov can be installed and tested through the Julia package manager:

julia> ]
pkg> add Krylov
pkg> test Krylov

Bug reports and discussions

If you think you found a bug, feel free to open an issue. Focused suggestions and requests can also be opened as issues. Before opening a pull request, start an issue or a discussion on the topic, please.

If you want to ask a question not suited for a bug report, feel free to start a discussion here. This forum is for general discussion about this repository and the JuliaSmoothOptimizers organization, so questions about any of our packages are welcome.

krylov.jl's Issues

Update CR

Update CR to incorporate @DahitoMA's modifications.

Use operators with preallocation when possible

test_cg.jl gives an example of operator with preallocation. All methods should use one whenever possible. This requires updating variants.jl and checking that there are no side effects due to reusing storage.

Preconditioners

We should add preconditioners in the obvious places (e.g., cg()).

Preconditioners roadmap

• CG
• CG_LANCZOS
• CG_LANCZOS_SHIFT_SEQ
• CR
• MINRES
• SYMMLQ
• CGLS
• CGNE
• CRAIG
• CRAIGMR
• CRMR
• LSLQ
• LSMR
• LSQR
• CRLS

zero rhs - return is inconsistent

At the end of the method, it returns x and stats, but in some places, it returns only x (e.g., cg when b is zero). Is this by design, or should I change it?

USYMLQR

It will be great to have this method for solving saddle-point and SQD systems.
It will also complete the family of "LQR" methods with BiLQR and TriLQR.
A basic functional version would be perfect.

Add a boolean option history to store or not the residual and Aresidual norms history

Krylov methods

• CG
• CR
• SYMMLQ
• CG-LANCZOS
• CG-LANCZOS-SHIFT-SEQ
• MINRES
• MINRES-QLP
• DQGMRES
• DIOM
• USYMLQ
• USYMQR
• TRICG
• TRIMR
• TRILQR
• CGS
• BICGSTAB
• BILQ
• QMR
• BILQR
• CGLS
• CRLS
• CGNE
• CRMR
• LSLQ
• LSQR
• LSMR
• LNLQ
• CRAIG
• CRAIGMR

Preallocate workspace for solvers

Krylov methods

• CG
• CR
• SYMMLQ
• CG-LANCZOS
• CG-LANCZOS-SHIFT-SEQ
• MINRES
• MINRES-QLP
• DQGMRES
• DIOM
• USYMLQ
• USYMQR
• TRICG
• TRIMR
• TRILQR
• CGS
• BICGSTAB
• BILQ
• QMR
• BILQR
• CGLS
• CRLS
• CGNE
• CRMR
• LSLQ
• LSQR
• LSMR
• LNLQ
• CRAIG
• CRAIGMR

For MINRES, this could take the form (not tested)

function minres(A, b :: AbstractVector{T}; kwargs...) where T <: AbstractFloat
  x = similar(b)
  r1 = similar(b)
  window = kwargs.get(:window, 5)
  err_vec = Vector{T}(undef, window)
  minres!(A, b, x, r1, err_vec; kwargs...)
end

function minres!(A, b, x, r1, err_vec;
                M=opEye(), λ :: T=zero(T), atol :: T=√eps(T)/100,
                rtol :: T=√eps(T)/100, etol :: T=√eps(T),
                window :: Int=5, itmax :: Int=0, conlim :: T=1/√eps(T),
                verbose :: Int=0) where T <: AbstractFloat
  ...
  x .= 0
  r1 .= b
  err_vec .= 0
  ...
end

The objective is for minres!() to not allocate at all. We could also define

abstract type KrylovSolver end
mutable struct MinresSolver <: KrylovSolver
  x
  r1
  err_vec
  # maybe default parameter values...
  function MinresSolver(A, b)
    new(...)
  end
end

and pass an instance of MinresSolver around.

@amontoison @geoffroyleconte

GPU Documention Broken

Hi, the GPU example for a general square system appears to be broken.

using CUDA, Krylov
using CUDA.CUSPARSE, SparseArrays

# LU ≈ A for CuSparseMatrixCSC matrices
P = ilu02(A_gpu, 'O')

# Solve Py = x
function ldiv!(y, P, x)
  copyto!(y, x)                        
  sv2!('N', 'L', 'N', 1.0, P, y, 'O')     # <---- 3rd entry (`N`) not recognized 
  sv2!('N', 'U', 'U', 1.0, P, y, 'O')    # <---- 3rd entry (`U`) not recognized 
  return y
end

# Operator that model P⁻¹
y = similar(b_gpu); n = length(b_gpu); T = eltype(b_gpu)
opM = LinearOperator(T, n, n, false, false, x -> ldiv!(y, P, x))

# Solve an unsymmetric system with an incomplete LU preconditioner on GPU
(x, stats) = bicgstab(A_gpu, b_gpu, M=opM)

The arrowed lines appear to be the issue. Is this equivalent to sv2!('N', 'L', 1.0, P, y, 'O') and sv2!('N', 'U', 1.0, P, y, 'O') respectively?

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

LU with pivoting -- DIOM

An option was added in DIOM to give (or not) the possibility of row pivoting (useful if H_k is singular) #169 .
But we should derive a residual norm estimate suitable for this case.

LQstats

We should add a LQstats structure for LSLQ, LNLQ, SYMMLQ, BiLQ. #138 .

CRAIG should stop based on the backward error

The stopping tests should be as in LSQR, etc., and should include the backward error rNorm / sqrt(β₁^2 + Anorm² * xNorm²).

ERROR: LoadError: UndefVarError: distribute not defined

Hi,

When I try to use the cg_lanczos_shift_par function I get the error:

ERROR: LoadError: UndefVarError: distribute not defined

Thank you very much for your help!

IterativeSolvers.jl

Did you have a look at https://github.com/JuliaLang/IterativeSolvers.jl?
There seems to be a lot of good things here (especially in Dominique's branch) that could be added there

Implement transfer QMR -> BiCG

See Freund & Nachtigal p.329. What criterion should we use for the transfer? Residual? Error? Other?

Flag for enforcing positive definiteness or linesearch

Preconditioned USYMLQ and USYMQR

With the trick found for storing vₖ in MINRES-QLP (#200), we should be able to add preconditioning in elliptic norms for USYMLQ and USYMQR.

Constructors KrylovSolver

Thanks for the amazing work including the KrylovSolver.

I just find having MinresSolver(A, b) as a constructor not practical since only the size of A and the type of b is needed.
For instance, in an optimization solver, we typically want to initialize the KrylovSolver at the initialization step where the size and the type is known but we don't have access yet to the matrix and the rhs.

Maybe we could have a MinresSolver(m, n, S) by default and an additional one with MinresSolver(A, b). What do you think?

CRMR tests sometimes pass, sometimes fail

Compare
https://travis-ci.org/JuliaOptimizers/Krylov.jl/builds/128475745
and
https://travis-ci.org/jagot/Krylov.jl/builds/128475739

which both test the same commit.

cg returns NaN for singular problem

Q = [1.0 0.0; 0.0 0.0]
g = -Q * ones(2) + [0.0; 1.0]
x, ks = Krylov.cg(Q, -g)

Output:

([NaN, NaN], 
Simple stats
  solved: false
  inconsistent: false
  residuals:  [ 1.4e+00  1.4e+00      NaN      NaN      NaN ]
  Aresiduals: []
  status: maximum number of iterations exceeded
)

Preconditioned LSMR

I noticed that https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/master/src/lsmr.jl has preconditioning, but couldn't figure out how to use the preconditioning to solve a simple LLS problem. Do you have a MWE that shows how to use preconditioners (e.g. something like modifying)

N = 10
σ = 0.1
X = rand(N,3)
c = rand(3)
y = X * c + σ * randn(N)
c_estimate = X \ y
@show cond(X' * X)
Krylov.lsmr(X, y)  # how to add preconditioning

Also, and this is more of a theoretical question.... but does preconditioning typically help if the normal equations are ill-conditioned? Are there typical tricks which work for those setups, etc.

I thought someone suggested that algebraic multigrid helped precondition the sorts of matrices that come about in LLS, but I could be wildly wrong.

Standardize stopping tests

All methods should have stopping tests based on tolerances, on the floating-point arithmetic, and (if possible) on backward-error measures.

High memory allocation when using preconditioners

I have noticed some higher-than-expected memory allocations when using a preconditioner other than opEye().

I only tested minres and cg so far, and the culprit seems to be a matrix-vector product with M, namely, in cg and minres.
If I replace, e.g., v = M * r2 in minres by mul!(v, M, r2), the allocations are more reasonable (see example below).

using Krylov, LinearAlgebra, SparseArrays, Random
using BenchmarkTools

Random.seed!(0)
m = 2^10
A = sprand(m, m, 0.01)
S = A+A'
b = ones(m)

M1 = Krylov.opEye()
M2 = 1.0I
M3 = Diagonal(ones(m))

# Warm-up will be done during calibration
@btime minres($S, $b)         # No preconditioner
@btime minres($S, $b, M=$M1)  # M = opEye
@btime minres($S, $b, M=$M2)  # M = I (uniform scaling)
@btime minres($S, $b, M=$M3)  # M = Diagonal(1.0, ..., 1.0)

Current version

  19.906 ms (34 allocations: 114.13 KiB)
  19.902 ms (34 allocations: 114.13 KiB)
  21.170 ms (1059 allocations: 8.24 MiB)
  20.397 ms (54 allocations: 122.78 KiB)

and if I replace with mul!(v, M, r2)

  20.147 ms (34 allocations: 114.13 KiB)
  20.006 ms (34 allocations: 114.13 KiB)
  20.654 ms (35 allocations: 122.25 KiB)
  20.344 ms (54 allocations: 122.78 KiB)

Same behavior with cg.

EDIT: making the above modification to minres causes some tests to fail

problems running the examples

I am trying to run the example with a general square system on a gpu and having some problems.

I was copying and pasting the lines from above and found out that I needed to make A_gpu a square matrix, so I made it 200 x 200 using the sprand command above. This seemes to work better but then I get an error about zero pivot, see below.

Very sorry for not understanding this better but is there an example you might suggest that would work better?

julia> P = ic02(A_gpu, 'O')
ERROR: CUSPARSEError: zero pivot (code 9, CUSPARSE_STATUS_ZERO_PIVOT)
Stacktrace:
 [1] throw_api_error(::CUDA.CUSPARSE.cusparseStatus_t) at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/cusparse/error.jl:21
 [2] macro expansion at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/cusparse/error.jl:32 [inlined]
 [3] cusparseXcsric02_zeroPivot(::Ptr{Nothing}, ::Ptr{Nothing}, ::Base.RefValue{Int32}) at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/utils/call.jl:93
 [4] macro expansion at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/cusparse/wrappers.jl:776 [inlined]
 [5] macro expansion at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/utils/call.jl:211 [inlined]
 [6] ic02!(::CuSparseMatrixCSC{Float64}, ::Char) at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/cusparse/wrappers.jl:767
 [7] ic02(::CuSparseMatrixCSC{Float64}, ::Char) at /home/fpoulin/.julia/packages/CUDA/dZvbp/lib/cusparse/wrappers.jl:958
 [8] top-level scope at REPL[24]:1

minres_qlp performs scalar operations on GPU

Something is happening inside minres_qlp that seems to not be GPU-friendly

using Random, LinearAlgebra, CUDA, Krylov

Random.seed!(0)
n = 64
A = Matrix(Symmetric(rand(n, n)))
cA = CuArray(A)

b = ones(n)
cb = CuArray(b)

x, stats = minres_qlp(A, b);
cx, cstats = minres_qlp(cA, cb);

The last call leads to this warning:

┌ Warning: Performing scalar operations on GPU arrays: This is very slow, consider disallowing these operations with `allowscalar(false)`
└ @ GPUArrays ~/.julia/packages/GPUArrays/eVYIC/src/host/indexing.jl:43

Subsequent timings look like

julia> @time x, stats = minres_qlp(A, b);
  0.000091 seconds (28 allocations: 8.391 KiB)

julia> CUDA.@time cx, stats = minres_qlp(cA, cb);
  0.088764 seconds (91.31 k CPU allocations: 3.673 MiB) (6 GPU allocations: 3.000 KiB, 0.06% gc time of which 79.26% spent allocating)

---

julia> versioninfo()
Julia Version 1.5.1
Commit 697e782ab8 (2020-08-25 20:08 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-6600 CPU @ 3.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)

Estimates of the error in the LNLQ method

Documentation and examples

• Documentation
• Examples

Complex linear systems

Currently, we are only able to solve linear systems with Float32 and Float64 types.
It will be relevant to generalized them to Float16, BigFloat (#99) and Complex types.

Krylov methods

• CG (#490)
• CR (#508)
• SYMMLQ (#506)
• CG-LANCZOS without shifts (#509)
• CG-LANCZOS with shifts (#509)
• MINRES (#504)
• MINRES-QLP (#505)
• DQGMRES (#491)
• DIOM (#489)
• GMRES (#481)
• FOM (#486)
• USYMLQ (#496)
• USYMQR (#495)
• TRICG (#492)
• TRIMR (#494)
• TRILQR (#498)
• CGS (#507)
• BICGSTAB (#481)
• BILQ (#497)
• QMR (#481)
• BILQR (#517)
• CGLS (#500)
• CRLS (#501)
• CGNE (#503)
• CRMR (#502)
• LSLQ (#515)
• LSQR (#514)
• LSMR (#513)
• LNLQ (#510)
• CRAIG (#511)
• CRAIGMR (#512)
• GPMR (#483)

`@kzeros` on derivative Array types

I noticed some solvers return an error when the right-hand side is a derivative type (using Julia 1.6). For instance:

using Krylov
a = rand(10)
b = view(a, 1:4)
Krylov.kzeros(typeof(b), 4)

return

julia> Krylov.kzeros(typeof(b), 4)
ERROR: MethodError: no method matching SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}(::UndefInitializer, ::Int64)
Closest candidates are:
  SubArray{T, N, P, I, L}(::Any, ::Any, ::Any, ::Any) where {T, N, P, I, L} at subarray.jl:19
Stacktrace:
 [1] kzeros(S::Type{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, n::Int64)
   @ Krylov ~/.julia/packages/Krylov/cVOOc/src/krylov_utils.jl:174
 [2] top-level scope
   @ REPL[22]:1

Problem with Matrix-vector operation in CRMR

I tried to understand why the test with a preconditioner didn't passed with Julia 0.7. And the problem seems to come from the new version of LinearOperators.

With Julia 0.6 :

CRMR: system of 2 equations in 10 variables
Aprod     ‖A'r‖       ‖r‖
    1  5.48e-02  3.01e-03
    3  1.35e-02  8.68e-03
    5  2.30e-16  3.70e-17
CRMR: residual: 2.6e-15
‖x - xmin‖ = 8.8e-17

With Julia 0.7

CRMR: system of 2 equations in 10 variables
Aprod     ‖A'r‖       ‖r‖
    1  5.48e-02  3.01e-03
    3  4.25e-02  6.58e-03
    5  1.14e-01  1.84e-02
    7  1.20e+00  1.94e-01
    9  1.36e+02  2.20e+01
   11  1.75e+06  2.84e+05
   13  2.91e+14  4.71e+13
   15  8.04e+30  1.30e+30
   17  6.12e+63  9.91e+62
   19  3.55e+129  5.75e+128
   21       Inf  1.93e+260
   23       NaN       NaN
   25       NaN       NaN
CRMR: residual:     NaN

The Aᵗp operation returns a wrong value (without preconditioner).
I also tested with a CRMR version that allow matrix instead of LinearOperator with julia 0.7. And this time we have the good result.

crmr2(A, b, M=M, verbose=true)
CRMR: system of 2 equations in 10 variables
Aprod     ‖A'r‖       ‖r‖
    1  5.48e-02  3.01e-03
    3  1.35e-02  8.68e-03
    5  2.30e-16  3.70e-17

Allow Github Actions failures with Julia nightly

Cf JuliaSmoothOptimizers/SolverTest.jl#6

Broadcast

Broadcast is not completely optimized on GPU (#242) but it could give some performance benefit by combining multiple CUBLAS calls into a single kernel in the future. Broadcast allows the code to be generic and we should track the performance of it.

QMR for inconsistent systems

Seems like we need to develop a bound on ||A'r||, and maybe compute r explicitly.

Add function signature in documentation

The current documentation of user-facing functions, e.g., Krylov.cg, Krylov.minres, etc., does not include these functions' signature.

It makes it hard for newbies like me to know how exactly to call a function, without going through the source code.

Here is an example of what I had in mind, for Krylov.cg:

---

cg(A, b; kwargs...)

The conjugate gradient method to solve the symmetric linear system Ax=b.

The method does not abort if A is not definite.

A preconditioner M may be provided in the form of a linear operator and is
assumed to be symmetric and positive definite.
M also indicates the weighted norm in which residuals are measured.

Arguments

• A: Left-hand side
• b::AbstractVector: Right-hand side

Key-word arguments

• M: user-provided preconditioner
• atol: absolute tolerance
• rtol: relative tolerance
• itmax::Int: maximum number of iterations
  ...

---

Even just the signature would be useful.

Sign of the regularization in minres

I suspect there is a bug with the sign of the regularization in minres. Following the doc, it should solve the shifted linear system

(A + λ I) x = b

Taking A = I, b = 1, and λ=0.5 should return a vector with 2/3. However, it returns 2, which would correspond to put -λ.

using Krylov, LinearAlgebra
In = diagm(0 => ones(3))
b = ones(3)
(x, stats) = minres(In, b, λ = 0.5) # return ≈2*ones(3)
(x, stats) = minres(In, b, λ = -0.5) # return ≈2/3*ones(3)

Precision is lost when using matrices

Precision is lost when using matrices.

A = rand(BigFloat, 3, 3)
A = A * A'
b = rand(BigFloat, 3)
x = cg(A, b)[1]
eltype(x)

returns Float64. I couldn't fix, because just adding types in variants.jl breaks other places.

Preconditioned CR

With #84, preconditioned CR doesn't work anymore even if we are not in a linesearch or trust-region context : lineseach = false and radius = 0.
#54

Substitute A.tprod by A'

I did it in a branch https://github.com/amontoison/Krylov.jl/tree/transpose_product
But we need to solve JuliaSmoothOptimizers/LinearOperators.jl#92 before.

If b is an array of Integers, doesn't work.

Not sure which changes yet, but it stopped working.

Allocations tests

We should check the allocations of each solver (old and future ones). It's a good way to find mistakes and verify that implementations perform well.

Allocations tests roadmap

Methods already in Krylov.jl

• CG_LANCZOS
• CG_LANCZOS_SHIFT_SEQ
• CG
• CR
• CGLS
• CRLS
• CGNE
• CRAIG
• CRMR
• CRAIGMR
• LSLQ
• LSMR
• LSQR
• MINRES
• SYMMLQ
• DQGMRES
• DIOM
• CGS

Future methods

• USYMLQ
• USYMQR
• TFQMR
• GMRES
• FOM

Performance loss with Julia 0.6

Example on matrix bcsstk18 from the UFL collection in MatrixMarket format:

julia> @benchmark xs, stats = cg_lanczos_shift_seq(A, b, [1,2,3])
BenchmarkTools.Trial: 
  memory estimate:  22.68 GiB
  allocs estimate:  908205
  --------------
  minimum time:     7.600 s (15.76% GC)
  median time:      7.600 s (15.76% GC)
  mean time:        7.600 s (15.76% GC)
  maximum time:     7.600 s (15.76% GC)
  --------------
  samples:          1
  evals/sample:     1

julia> VERSION
v"0.5.2"

but

julia> @benchmark xs, stats = cg_lanczos_shift_seq(A, b, [1,2,3])
BenchmarkTools.Trial: 
  memory estimate:  23.06 GiB
  allocs estimate:  2949084
  --------------
  minimum time:     10.638 s (12.18% GC)
  median time:      10.638 s (12.18% GC)
  mean time:        10.638 s (12.18% GC)
  maximum time:     10.638 s (12.18% GC)
  --------------
  samples:          1
  evals/sample:     1

julia> VERSION
v"0.6.0"

Krylov 1.0

Epic created through ZenHub to control what is needed for the Krylov 1.0 release.

Use mul! for in-place solvers

Krylov methods

• CG
• CR
• SYMMLQ
• CG-LANCZOS
• CG-LANCZOS-SHIFT-SEQ
• MINRES
• MINRES-QLP
• DQGMRES
• DIOM
• USYMLQ
• USYMQR
• TRICG
• TRIMR
• TRILQR
• CGS
• BICGSTAB
• BILQ
• QMR
• BILQR
• CGLS
• CRLS
• CGNE
• CRMR
• LSLQ
• LSQR
• LSMR
• LNLQ
• CRAIG
• CRAIGMR

Add documentation about preconditioners

It will be useful to show how to build relevant preconditioners for our Krylov methods in the documentation with the help of linear operators.

• Jacobi
• SPAI
• Incomplete L |D| Lᵀ
• ILU
• LLᵀ
• LDLᵀ

Allow more flexibility in verbose option

Currently, Krylov.jl allows to control the verbosity using the verbose option, common to each method. If set to true, the algorithm prints out the trace at each iteration.
I am wondering if we could allow finer-grained control to the verbose option? E.g. to print the trace every 10 or 50 iterations? I see two ways forward:

• add a new argument verbose_it::Int specifying how often we should print out the trace if verbose=true.
• update verbose to allow it to take non-negative integer values (that would incur a breaking change, though)

First Release

Krylov doesn't have a release yet (on METADATA), we should make one supporting 0.6 and 0.7.
What needs to be included?

• #62
• #63
• #67

Statuses

All statuses should be symbols.

Missing methods

julia> craig(A', b)
ERROR: MethodError: no method matching craig(::Adjoint{Float64,Array{Float64,2}}, ::Array{Float64,1})

Add iter on stats

Add iter on stats (cg.jl, at least).

Eliminate random tests