define mul! for quasi-Newton operators

.cirrus.yml

@@ -4,7 +4,6 @@ task:
   name: FreeBSD
   env:
     matrix:
-      - JULIA_VERSION: 1.0
       - JULIA_VERSION: 1
       - JULIA_VERSION: nightly
   allow_failures: $JULIA_VERSION == 'nightly'

Project.toml

@@ -12,11 +12,12 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 [compat]
 FastClosures = "0.2, 0.3"
 TimerOutputs = "^0.5"
-julia = "^1.0.0"
+julia = "^1.3.0"
 
 [extras]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
 
 [targets]
-test = ["Test", "TestSetExtensions"]
+test = ["Arpack", "Test", "TestSetExtensions"]

src/lbfgs.jl

@@ -292,3 +292,10 @@ function reset!(op :: LBFGSOperator)
   op.nctprod = 0
   return op
 end
+
+# define mul! so we can call, e.g., Arpack
+function mul!(y :: AbstractVector, op :: LBFGSOperator, x :: AbstractVector)
+  op.prod(x)
+  y .= op.data.Ax
+  return y
+end

src/lsr1.jl

@@ -214,3 +214,10 @@ function reset!(op :: LSR1Operator)
   op.nctprod = 0
   return op
 end
+
+# define mul! so we can call, e.g., Arpack
+function mul!(y :: AbstractVector, op :: LSR1Operator, x :: AbstractVector)
+  op.prod(x)
+  y .= op.data.Ax
+  return y
+end

test/runtests.jl

@@ -1,4 +1,4 @@
-using Test, TestSetExtensions, LinearOperators
+using Arpack, Test, TestSetExtensions, LinearOperators
 
 include("test_aux.jl")
 

test/test_lbfgs.jl

@@ -182,6 +182,23 @@ function test_lbfgs()
     @test nallocs == 0
   end
 
+  @testset "LBFGS eigenvalues" begin
+    n = 100
+    mem = 20
+    B = LBFGSOperator(n, mem=mem)
+    H = InverseLBFGSOperator(n, mem=mem)
+    for _ = 1:2:n
+      s = rand(n)
+      y = rand(n)
+      push!(B, s, y)
+      push!(H, s, y)
+    end
+    λ = eigs(B, nev=n)[1]  # call Arpack
+    @test minimum(λ) > 0
+    λ = eigs(H)[1]
+    @test minimum(λ) > 0
+  end
+
 end
 
 test_lbfgs()

test/test_lsr1.jl

@@ -94,6 +94,20 @@ function test_lsr1()
     nallocs = @allocated diag!(B, x)
     @test nallocs == 0
   end
+
+  @testset "LSR1 eigenvalues" begin
+    n = 100
+    mem = 20
+    B = LSR1Operator(n, mem=mem)
+    for _ = 1:2:n
+      s = rand(n)
+      y = rand(n)
+      push!(B, s, y)
+    end
+    vals = eigs(B, nev=n)
+    resid = vals[end]
+    @test norm(resid) ≤ sqrt(eps(eltype(B)))
+  end
 end
 
 test_lsr1()