diff --git a/Manifest.toml b/Manifest.toml
index 32b5805..3a72c6b 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,5 +1,20 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[Arpack]]
+deps = ["Arpack_jll", "Libdl", "LinearAlgebra", "Logging"]
+git-tree-sha1 = "9b9b347613394885fd1c8c7729bfc60528faa436"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.5.4"
+
+[[Arpack_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "5ba6c757e8feccf03a1554dfaf3e26b3cfc7fd5e"
+uuid = "68821587-b530-5797-8361-c406ea357684"
+version = "3.5.1+1"
+
 [[Artifacts]]
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
@@ -10,10 +25,40 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.4.1"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
@@ -27,18 +72,45 @@ git-tree-sha1 = "d1b46faefb7c2f48fdec69e6f3cc34857769bc15"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
 version = "3.8.0"
 
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
 [[Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
 [[OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "7eb1686b4f04b82f96ed7a4ea5890a4f0c7a09f1"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.4.0"
+
 [[Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
 [[Random]]
 deps = ["SHA", "Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -49,6 +121,9 @@ uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
 [[SparseArrays]]
 deps = ["LinearAlgebra", "Random"]
 uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -57,9 +132,33 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 deps = ["LinearAlgebra", "SparseArrays"]
 uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
 [[libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/Project.toml b/Project.toml
index 30a9014..9ed987f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,6 +4,7 @@ authors = ["Nick Mayhall <nmayhall@vt.edu> and contributors"]
 version = "0.1.0"
 
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
diff --git a/README.md b/README.md
index 8fa6fde..4b4c245 100644
--- a/README.md
+++ b/README.md
@@ -2,3 +2,38 @@
 
 [![Build Status](https://github.com/nmayhall-vt/BlockDavidson.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/nmayhall-vt/BlockDavidson.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/nmayhall-vt/BlockDavidson.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/nmayhall-vt/BlockDavidson.jl)
+
+---
+Simple code to solve for the eigenvectors/eigenvalues of a matrix or `LinearMap`. Preconditioner has not yet been implemented, so this is currently just a simple block lanczos code. 
+
+## Example usage 
+
+### Explicit matrix diagonalization
+
+```julia
+dav = Davidson(A)
+e,v = eigs(dav)
+```
+to use the diagonal preconditioner specify as argument to eigs
+```julia
+e,v = eigs(dav, Adiag=diag(A))
+```
+
+### Matrix-free diagonalization 
+
+We can also diagonalize an implicit matrix by defining a function `matvec` that algorithmically implments the action of `A` onto a vector or set of vectors. We can also specify several settings, including providing an initial guess `v_guess`:
+```julia
+using LinearMaps
+function matvec(v)
+    return A*v. # implement this however is appropriate
+end
+
+
+lmap = LinearMap(matvec)
+dav = Davidson(lmap; max_iter=200, max_ss_vecs=8, tol=1e-6, nroots=6, v0=v_guess, lindep_thresh=1e-10)
+e,v = eigs(dav)
+```
+
+Convergence issues:
+- the implementation is still not as robust as it should be, with loss of orthogonality making it very difficult to converge tightly in some scenarios. If this is the case, try tightening `lindep_thresh`. Whenever the determinant of the subspace drops below this value, the algorithm is restarted. As such, if it's too loose, then the correction vectors aren't very accurate after orthogonalization. If it's too tight, then it restarts too frequently, slowing down convergence. 
+
diff --git a/src/BlockDavidson.jl b/src/BlockDavidson.jl
index c04ecc8..d457139 100644
--- a/src/BlockDavidson.jl
+++ b/src/BlockDavidson.jl
@@ -3,19 +3,11 @@ using LinearAlgebra
 using Printf
 using InteractiveUtils
 
-# solve is exported by several packages.
-# check if its defined, and if so, extend it
-# this doesn't work
-#if @isdefined(solve)
-#    println(" Solve already defined: extend it. ", @which solve)
-#    a = @which solve
-#    eval(Meta.parse("import $a.solve"))
-#end    
-export solve
+export LinOpMat 
 export Davidson
+export eigs
 
-
-mutable struct Davidson
+mutable struct Davidson{T}
     op 
     dim::Int
     nroots::Int
@@ -25,27 +17,49 @@ mutable struct Davidson
     converged::Bool
     status::Vector{Bool}        # converged status of each root 
     iter_curr::Int
-    vec_curr::Array{Float64,2}
-    sig_curr::Array{Float64,2}
-    vec_prev::Array{Float64,2}
-    sig_prev::Array{Float64,2}
-    ritz_v::Array{Float64,2}
-    ritz_e::Vector{Float64}
-    resid::Vector{Float64}
+    vec_curr::Array{T,2}
+    sig_curr::Array{T,2}
+    vec_prev::Array{T,2}
+    sig_prev::Array{T,2}
+    ritz_v::Array{T,2}
+    ritz_e::Vector{T}
+    resid::Vector{T}
     lindep::Float64
     lindep_thresh::Float64
 end
 
-function Davidson(op; max_iter=100, max_ss_vecs=8, tol=1e-8, nroots=1, v0=nothing, lindep_thresh=1e-10)
-    size(op)[1] == size(op)[2] || throw(DimensionError())
+"""
+    Davidson(op; 
+    max_iter=100, 
+    max_ss_vecs=8, 
+    tol=1e-8, 
+    nroots=1, 
+    v0=nothing, 
+    lindep_thresh=1e-10, 
+    T=Float64)
+
+TBW
+"""
+function Davidson(op; 
+    max_iter=100, 
+    max_ss_vecs=8, 
+    tol=1e-8, 
+    nroots=1, 
+    v0=nothing, 
+    lindep_thresh=1e-10, 
+    T=Float64)
+
+    size(op)[1] == size(op)[2] || throw(DimensionMismatch)
     dim = size(op)[1]
-    if v0==nothing
-        v0 = rand(dim,nroots)
-        S = v0'*v0
-        v0 = v0*inv(sqrt(S))
+    if v0===nothing
+        F = qr(rand(T, dim,nroots))
+        v0 = Matrix(F.Q)
+    else
+        size(v0,1) == size(op,1) || throw(DimensionMismatch)
+        size(v0,2) == nroots || throw(DimensionMismatch)
     end
     #display(v0'*v0)
-    return Davidson(op, 
+    return Davidson{T}(op, 
                     dim, 
                     nroots, 
                     max_iter,
@@ -56,11 +70,11 @@ function Davidson(op; max_iter=100, max_ss_vecs=8, tol=1e-8, nroots=1, v0=nothin
                     0,
                     v0,
                     v0,
-                    zeros(dim,0),
-                    zeros(dim,0),
-                    zeros(nroots,nroots),
-                    zeros(nroots),
-                    zeros(nroots),
+                    zeros(T,dim,0),
+                    zeros(T,dim,0),
+                    zeros(T,nroots,nroots),
+                    zeros(T,nroots),
+                    zeros(T,nroots),
                     1.0,
                     lindep_thresh)
 end
@@ -104,19 +118,31 @@ function print_iter(solver::Davidson)
 end
 
 
-function _apply_diagonal_precond!(res_s::Vector{Float64}, Adiag::Vector{Float64}, denom::Float64)
-    dim = size(Adiag)[1]
-    length(size(Adiag)) == 1 || throw(Exception)
-    length(size(res_s)) == 1 || throw(Exception)
-    @inbounds @simd for i in 1:dim
+"""
+    _apply_diagonal_precond!(res_s::Vector{T}, Adiag::Vector{T}, denom::T) where {T}
+
+TBW
+"""
+function apply_diagonal_precond!(res_s::Vector{T}, Adiag::Vector{T}, denom::T) where {T}
+    dim = length(Adiag)
+    length(res_s) == length(Adiag) || throw(DimensionMismatch)
+    # res_s .= -100 .* res_s ./ (Adiag .- denom)
+    
+    # @inbounds @simd 
+    for i in 1:dim
         res_s[i] = res_s[i] / (denom - Adiag[i])
     end
 end
 
-function iteration(solver::Davidson; Adiag=nothing, iprint=0)
+function iteration(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1.0)
+
+    # 
+    # project out v_prev from v_curr
+    solver.vec_curr = solver.vec_curr - solver.vec_prev * (solver.vec_prev' * solver.vec_curr)
+    solver.vec_curr = Matrix(qr(solver.vec_curr).Q) 
     
     #
-    # perform sig_curr = A*vec_curr 
+    # perform σ_curr = A*v_curr 
     solver.sig_curr = Matrix(solver.op * solver.vec_curr)
    
     #
@@ -133,27 +159,19 @@ function iteration(solver::Davidson; Adiag=nothing, iprint=0)
     # form H in current subspace
     Hss = solver.vec_prev' * solver.sig_prev
     F = eigen(Hss)
-    #ritz_e = F.values[1:solver.nroots]
-    #ritz_v = F.vectors[:,1:solver.nroots]
-    ritz_e = F.values
-    ritz_v = F.vectors
-    ss_size = length(F.values)
-    if ss_size >= solver.max_ss_vecs
-        ritz_v = ritz_v[:,sortperm(ritz_e)][:,1:solver.nroots]
-        ritz_e = ritz_e[sortperm(ritz_e)][1:1:solver.nroots]
-        #ritz_v = ritz_v[:,sortperm(ritz_e)][:,1:solver.max_ss_vecs]
-        #ritz_e = ritz_e[sortperm(ritz_e)][1:1:solver.max_ss_vecs]
-    else
-        ritz_v = ritz_v[:,sortperm(ritz_e)]
-        ritz_e = ritz_e[sortperm(ritz_e)]
-    end
+    idx = sortperm(F.values)
 
-    solver.ritz_e = ritz_e
-    solver.ritz_v = ritz_v
+    ss_size = min(size(solver.vec_prev,2), solver.max_ss_vecs)
+    
+    idx = idx[1:ss_size]
 
-    solver.sig_prev = solver.sig_prev * ritz_v 
-    solver.vec_prev = solver.vec_prev * ritz_v
-    Hss = ritz_v' * Hss * ritz_v
+    solver.ritz_e = F.values[idx]
+    solver.ritz_v = F.vectors[:,idx]
+    
+    
+    solver.sig_prev = solver.sig_prev * solver.ritz_v 
+    solver.vec_prev = solver.vec_prev * solver.ritz_v
+    Hss = solver.ritz_v' * Hss * solver.ritz_v
 
     res = deepcopy(solver.sig_prev[:,1:solver.nroots])
     for s in 1:solver.nroots
@@ -161,13 +179,10 @@ function iteration(solver::Davidson; Adiag=nothing, iprint=0)
     end
 
     
-    ss = deepcopy(solver.vec_prev)
+    
 
-    new = zeros(solver.dim, 0) 
     #solver.statusconv = [false for i in 1:solver.nroots]
     for s in 1:solver.nroots
-        # |v'> = (I-|SS><SS|)|v>
-        #      = |v> - |SS><SS|V>
         solver.resid[s] = norm(res[:,s])
         if norm(res[:,s]) <= solver.tol
             solver.status[s] = true
@@ -175,39 +190,48 @@ function iteration(solver::Davidson; Adiag=nothing, iprint=0)
         else
             solver.status[s] = false 
         end
-        if Adiag != nothing
-            level_shift = 1e-3
-            denom = (ritz_e[s]  + level_shift)
-            _apply_diagonal_precond!(res[:,s], Adiag, denom)
-        end
-        scr = ss' * res[:,s]
-        res[:,s] .= res[:,s] .- (ss * scr)
-        nres = norm(res[:,s])
-        if nres > solver.tol
-            ss = hcat(ss,res[:,s]/nres)
-            new = hcat(new,res[:,s]/nres)
+        if Adiag !== nothing && solver.status[s] == false && solver.resid[s] < precond_start_thresh 
+        # if Adiag != nothing && solver.status[s] == false
+            tmp = deepcopy(res)
+            for i in 1:length(Adiag)
+                res[i, s] = res[i, s] / (solver.ritz_e[s] - Adiag[i] + 1e-12)
+                # if abs(solver.ritz_e[s] - Adiag[i]) < solver.tol
+                #     res[i,s] = 0
+                # else
+                #     res[i,s] = res[i,s] / (solver.ritz_e[s] - Adiag[i])
+                # end
+            end
         end
     end
-    return new
+
+
+    res = res - solver.vec_prev * (solver.vec_prev' * res)
+    return Matrix(qr(res).Q) 
 end
 
-function solve(solver::Davidson; Adiag=nothing, iprint=0)
+"""
+    eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1e-1)
+
+TBW
+"""
+function eigs(solver::Davidson; Adiag=nothing, iprint=0, precond_start_thresh=1e-1)
 
     for iter = 1:solver.max_iter
         #@btime $solver.vec_curr = $iteration($solver)
-        solver.vec_curr = iteration(solver, Adiag=Adiag, iprint=iprint)
+        solver.vec_curr = iteration(solver, Adiag=Adiag, iprint=iprint, precond_start_thresh=precond_start_thresh)
         solver.iter_curr = iter
         print_iter(solver)
         if all(solver.status)
-            solver.converged == true
+            solver.converged = true
             break
         end
         if solver.lindep > solver.lindep_thresh && iter < solver.max_iter
-            @warn "Linear dependency detected. Restarting."
+            # @warn "Linear dependency detected. Restarting."
+            println(" Linear dependency detected. Restarting.")
             flush(stdout)
             F = qr(solver.vec_prev[:,1:solver.nroots])
-            solver.vec_curr = Array(F.Q)
-            solver.sig_curr = Array(F.Q)
+            solver.vec_curr = Matrix(F.Q)
+            solver.sig_curr = Matrix(F.Q)
             solver.vec_prev = zeros(solver.dim, 0) 
             solver.sig_prev =  zeros(solver.dim, 0)
             solver.ritz_v = zeros(solver.nroots,solver.nroots)
@@ -219,11 +243,4 @@ function solve(solver::Davidson; Adiag=nothing, iprint=0)
     #return solver.fritz_e[1:solver.nroots], solver.vec_prev*solver.ritz_v[:,1:solver.nroots]
 end
 
-#function orthogonalize!(solver)
-#
-#    # |vv_i>  = |v_i> - |w
-#    new_v = zeros(solver.dim, nroots)
-#    for r in 1:nroots
-#end
-
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index e729f2f..8e6ade1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,4 +3,5 @@ using Test
 
 @testset "BlockDavidson.jl" begin
     include("test01.jl")
+    include("test_precond_01.jl")
 end
diff --git a/test/test01.jl b/test/test01.jl
index 8e219cd..ed35d4e 100644
--- a/test/test01.jl
+++ b/test/test01.jl
@@ -14,15 +14,32 @@ using LinearMaps
     #display(F.values)
 
     dav = Davidson(A)
-    e,v = solve(dav)
+    e,v = eigs(dav)
     e_ref = -18.260037795157675
     @test isapprox(e[1], -18.260037795157675)
 
     lmap = LinearMap(A)
 
     dav = Davidson(lmap)
-    e,v = solve(dav)
+    e,v = eigs(dav)
     e_ref = -18.260037795157675
     @test isapprox(e[1], -18.260037795157675)
 
+    e_ref = [
+             -18.260037795157675
+             -17.9716411644818
+             -17.47598854674256
+             -17.184105784827796
+             -16.939563610543257
+             -16.83937885452674
+            ]
+    # now with more settings specified and roots
+    dav = Davidson(lmap; max_iter=200, max_ss_vecs=8, tol=1e-6, nroots=6)
+    e,v = eigs(dav)
+    @test all(isapprox.(e, e_ref, atol=1e-10))
+
+    display(v'*Matrix(lmap*v))
+    # now test with initial guess
+    e,v = eigs(Davidson(lmap; max_iter=2, max_ss_vecs=8, tol=1e-8, nroots=6, v0=v))
+    @test all(isapprox.(e, e_ref, atol=1e-10))
 end
diff --git a/test/test_precond_01.jl b/test/test_precond_01.jl
new file mode 100644
index 0000000..d528d3c
--- /dev/null
+++ b/test/test_precond_01.jl
@@ -0,0 +1,34 @@
+using BlockDavidson
+using LinearAlgebra
+using Random
+using Test
+
+@testset "BlockDavidson preconditioning" begin
+    Random.seed!(2)
+
+    N = 1_000
+    nr = 4
+    A = -10 .* diagm(rand(N)) + 0.001 * rand(N, N)
+    A += A'
+
+    v0 = qr(rand(N, 2)).Q[:, 1:nr]
+
+    idx = sortperm(diag(A))
+    v0 = Matrix(1.0I, N, N)[:, idx[1:nr]]
+
+
+    println(" No preconditioning")
+    dav1 = Davidson(A, max_iter=50, nroots=nr, tol=1e-8, v0=v0)
+    dav2 = Davidson(A, max_iter=50, nroots=nr, tol=1e-8, v0=v0)
+    e1, v1 = eigs(dav1)
+    println()
+    println(" Preconditioned")
+    e2, v2 = eigs(dav2, Adiag=diag(A))
+   
+    flush(stdout)
+    display(dav1.converged)
+    display(dav2.converged)
+    @test dav1.converged == false
+    @test dav2.converged == true
+
+end
\ No newline at end of file