From 4209dd43f332b10a151f860d0759cd9e7c1d2d70 Mon Sep 17 00:00:00 2001
From: thorek1 <thorek1@users.noreply.github.com>
Date: Tue, 3 Oct 2023 22:54:59 +0200
Subject: [PATCH] optimised matrix equation solvers for large probs

---
 src/MacroModelling.jl | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/src/MacroModelling.jl b/src/MacroModelling.jl
index 21e5e5f5..a7b91554 100644
--- a/src/MacroModelling.jl
+++ b/src/MacroModelling.jl
@@ -23,7 +23,7 @@ import AbstractDifferentiation as AD
 import SpeedMapping: speedmapping
 import REPL
 import Unicode
-import MatrixEquations
+import MatrixEquations # good overview: https://cscproxy.mpi-magdeburg.mpg.de/mpcsc/benner/talks/Benner-Melbourne2019.pdf
 # import NLboxsolve: nlboxsolve
 # using NamedArrays
 # using AxisKeys
@@ -3762,10 +3762,9 @@ function calculate_second_order_solution(∇₁::AbstractMatrix{<: Real}, #first
     push!(dimensions,size(C))
     push!(dimensions,size(X))
 
-    𝐒₂, solved = solve_matrix_equation_forward(values, coords = coordinates, dims = dimensions, solver = :iterative, sparse_output = true)
-    # 𝐒₂, solved = solve_matrix_equation_forward([vec(B) ;vec(C) ;vec(X)], dims = [size(B) ;size(C) ;size(X)], tol = tol)
-    # 𝐒₂, solved = solve_matrix_equation_AD([vec(B) ;vec(C) ;vec(X)], dims = [size(B) ;size(C) ;size(X)], sparse_output = true)
-    # 𝐒₂, solved = solve_matrix_equation_forward([vec(B) ;vec(C) ;vec(X)], dims = [size(B) ;size(C) ;size(X)], sparse_output = true)
+    solver = length(X.nzval) / length(X) < .1 ? :sylvester : :gmres
+
+    𝐒₂, solved = solve_matrix_equation_forward(values, coords = coordinates, dims = dimensions, solver = solver, sparse_output = true)
 
     𝐒₂ *= M₂.𝐔₂
 
@@ -3870,9 +3869,8 @@ function calculate_third_order_solution(∇₁::AbstractMatrix{<: Real}, #first
     push!(dimensions,size(B))
     push!(dimensions,size(C))
     push!(dimensions,size(X))
-    
 
-    𝐒₃, solved = solve_matrix_equation_forward(values, coords = coordinates, dims = dimensions, solver = :iterative, sparse_output = true)
+    𝐒₃, solved = solve_matrix_equation_forward(values, coords = coordinates, dims = dimensions, solver = :gmres, sparse_output = true)
 
     𝐒₃ *= M₃.𝐔₃
 
@@ -4447,7 +4445,7 @@ function solve_matrix_equation_forward(ABC::Vector{Float64};
         𝐂¹ = C
         while change > eps(Float32) && iter < 10000
             𝐂¹ = A * 𝐂 * B - C
-            if !(A isa DenseMatrix)
+            if !(𝐂¹ isa DenseMatrix)
                 droptol!(𝐂¹, eps())
             end
             if iter > 500
@@ -4476,10 +4474,10 @@ function solve_matrix_equation_forward(ABC::Vector{Float64};
         end
         solved = change < eps(Float32)
     elseif solver == :sylvester
-        𝐂 = sylvd(collect(-A),collect(B),-C)
+        𝐂 = MatrixEquations.sylvd(collect(-A),collect(B),-C)
         solved = isapprox(𝐂, A * 𝐂 * B - C, rtol = eps(Float32))
     elseif solver == :lyapunov
-        𝐂 = MatrixEquations.lyapd(A,-C)
+        𝐂 = MatrixEquations.lyapd(collect(A),-C)
         solved = isapprox(𝐂, A * 𝐂 * A' - C, rtol = eps(Float32))
     elseif solver == :speedmapping
         soll = speedmapping(collect(-C); m! = (X, x) -> X .= A * x * B - C, stabilize = true)