Add a test with a variable preconditioner

src/fgmres.jl

@@ -29,6 +29,7 @@ This implementation allows a left preconditioner M and a flexible right precondi
 A situation in which the preconditioner is "not constant" is when a relaxation-type method,
 a Chebyshev iteration or another Krylov subspace method is used as a preconditioner. 
 Compared to GMRES, there is no additional cost incurred in the arithmetic but the memory requirement almost doubles.
+Thus, GMRES is recommended if the right preconditioner N is identical as each iteration.
 
 Full reorthogonalization is available with the `reorthogonalization` option.
 
@@ -93,9 +94,8 @@ function fgmres!(solver :: FgmresSolver{T,FC,S}, A, b :: AbstractVector{FC};
   length(b) == m || error("Inconsistent problem size")
   (verbose > 0) && @printf("FGMRES: system of size %d\n", n)
 
-  # Check M = Iₙ and N = Iₙ
+  # Check M = Iₙ
   MisI = (M === I)
-  NisI = (N === I)
 
   # Check type consistency
   eltype(A) == FC || error("eltype(A) ≠ $FC")

src/precompile.jl

@@ -0,0 +1,60 @@
+# This precompile statements are a mix of hand-generated ones (e.g., to ensure
+# coverage across the breadth of common constraint types, as well as some
+# SnoopCompile ones that are expensive.)
+
+function _precompile_()
+    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
+    T = Float64
+    scalar_sets = (
+        MOI.LessThan{T},
+        MOI.GreaterThan{T},
+        MOI.EqualTo{T},
+        MOI.Interval{T},
+        MOI.Integer,
+        MOI.ZeroOne,
+        MOI.Semiinteger{T},
+        MOI.Semicontinuous{T},
+    )
+    scalar_functions = (
+        MOI.VariableIndex,
+        MOI.ScalarAffineFunction{T},
+        MOI.ScalarQuadraticFunction{T},
+    )
+    vector_sets = (
+        MOI.Nonnegatives,
+        MOI.Nonpositives,
+        MOI.Zeros,
+        MOI.Reals,
+        MOI.SecondOrderCone,
+        MOI.RotatedSecondOrderCone,
+        MOI.ExponentialCone,
+        MOI.DualExponentialCone,
+        MOI.PositiveSemidefiniteConeSquare,
+        MOI.PositiveSemidefiniteConeTriangle,
+    )
+    vector_functions = (
+        MOI.VectorOfVariables,
+        MOI.VectorAffineFunction{T},
+        MOI.VectorQuadraticFunction{T},
+    )
+    constraints = vcat(
+        [(F, S) for F in scalar_functions for S in scalar_sets],
+        [(F, S) for F in vector_functions for S in vector_sets],
+    )
+    for (F, S) in constraints
+        Base.precompile(
+            _moi_add_constraint,
+            (
+                MOIU.CachingOptimizer{
+                    MOI.AbstractOptimizer,
+                    MOIU.UniversalFallback{MOIU.Model{Float64}},
+                },
+                F,
+                S,
+            ),
+        )
+    end
+
+    Base.precompile(Tuple{typeof(model_string),Type,Model})   # time: 0.25714603
+    return
+end

test/test_fgmres.jl

@@ -1,3 +1,16 @@
+import LinearAlgebra.mul!
+
+mutable struct FlexiblePreconditioner{T,S}
+  D::Diagonal{T, S}
+  ω::T
+end
+
+function mul!(y::Vector, P::FlexiblePreconditioner, x::Vector)
+  P.ω = -P.ω
+  mul!(y, P.D, x)
+  y .*= P.ω
+end
+
 @testset "fgmres" begin
   fgmres_tol = 1.0e-6
 
@@ -128,18 +141,14 @@
         @test(stats.solved)
       end
 
-      # test callback function
-      # A, b = sparse_laplacian(FC=FC)
-      # solver = FgmresSolver(A, b)
-      # tol = 1.0e-1
-      # N = Diagonal(1 ./ diag(A))
-      # stor = StorageGetxRestartedGmres(solver, N = N)
-      # storage_vec = similar(b)
-      # fgmres!(solver, A, b, N = N, atol = 0.0, rtol = 0.0, restart = true, callback = solver -> restarted_fgmres_callback_n2(solver, A, b, stor, N, storage_vec, tol))
-      # @test solver.stats.status == "user-requested exit"
-      # @test norm(A * x - b) ≤ tol
-      #
-      # @test_throws TypeError fgmres(A, b, callback = solver -> "string", history = true)
+      A, b = polar_poisson(FC=FC)
+      J = inv(Diagonal(A))  # Jacobi preconditioner
+      N = FlexiblePreconditioner(J, 1.0)
+      (x, stats) = fgmres(A, b, N=N)
+      r = b - A * x
+      resid = norm(r) / norm(b)
+      @test(resid ≤ fgmres_tol)
+      @test(stats.solved)
     end
   end
 end

test/test_utils.jl

@@ -470,24 +470,3 @@ function restarted_gmres_callback_n2(solver::GmresSolver, A, b, stor, N, storage
   storage_vec .-= b
   return (norm(storage_vec) ≤ tol)
 end
-
-# Successive over-relaxation (SOR) method
-function sor!(x, A, b, ω, k)
-  x .= 0
-  n = length(x)
-  for iter = 1:k
-    for i = 1:n
-      sum1 = sum(A[i,j] * x[j] for j = 1:i-1; init = 0)
-      sum2 = sum(A[i,j] * x[j] for j = i+1:n; init = 0)
-      x[i] = (1 - ω) * x[i] + (ω / A[i,i]) * (b[i] - sum1 - sum2)
-    end
-  end
-  return x
-end
-
-function test_sor()
-  A = [4 -1 -6 0; -5 -4 10 8; 0 9 4 -2; 1 0 -7 5]
-  b = [2; 21; -12; -6]
-  ω = 0.5
-  return A, b, ω
-end