versioninfo()

using BenchmarkTools, MatrixDepot, IterativeSolvers, LinearAlgebra, SparseArrays

# Wathen matrix of dimension 30401 x 30401
A = matrixdepot("wathen", 100)

using UnicodePlots
spy(A)

# sparsity level
count(!iszero, A) / length(A)

# rhs
b = ones(size(A, 1))
# solve Ax=b by CG
xcg = cg(A, b);
@benchmark cg($A, $b)

using Preconditioners
@time p = CholeskyPreconditioner(A, 2)

# solver Ax=b by PCG
xpcg = cg(A, b, Pl=p)
# same answer?
norm(xcg - xpcg)

# PCG is >10 fold faster than CG
@benchmark cg($A, $b, Pl=$p)

using BenchmarkTools, IterativeSolvers, LinearAlgebra, Random, SparseArrays

Random.seed!(280) # seed
n, p = 10000, 5000
X = sprandn(n, p, 0.001) # iid standard normals with sparsity 0.01
β = ones(p)
y = X * β + randn(n)

β̂_qr = X \ y
# least squares by QR
@benchmark $X \ $y

β̂_lsqr = lsqr(X, y)
@show norm(β̂_qr - β̂_lsqr)
# least squares by lsqr
@benchmark lsqr($X, $y)

β̂_lsmr = lsmr(X, y)
@show norm(β̂_qr - β̂_lsmr)
# least squares by lsmr
@benchmark lsmr($X, $y)

using LinearMaps, IterativeSolvers

# Overwrite y with A * x
# left difference assuming periodic boundary conditions
function leftdiff!(y::AbstractVector, x::AbstractVector) 
    N = length(x)
    length(y) == N || throw(DimensionMismatch())
    @inbounds for i in 1:N
        y[i] = x[i] - x[mod1(i - 1, N)]
    end
    return y
end

# Overwrite y with A' * x
# minus right difference
function mrightdiff!(y::AbstractVector, x::AbstractVector) 
    N = length(x)
    length(y) == N || throw(DimensionMismatch())
    @inbounds for i in 1:N
        y[i] = x[i] - x[mod1(i + 1, N)]
    end
    return y
end

# define linear map
D = LinearMap{Float64}(leftdiff!, mrightdiff!, 100; ismutating=true) 

@show size(D)
v = ones(size(D, 2)) # vector of all 1s
@show D * v
@show D' * v;

Matrix(D)

using SparseArrays
sparse(D)

using UnicodePlots
spy(sparse(D))

using Arpack
Arpack.svds(D, nsv=3)

using LinearAlgebra
# check solution against the direct method for SVD
svdvals(Matrix(D))

Arpack.eigs(D'D, nev=3, which=:LM)