Compute determinant directly

tkf · Jul 19, 2018 · 2814580 · 2814580
1 parent 6dfc0e6
commit 2814580
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 5 deletions.
diff --git a/src/codim2/analysis.jl b/src/codim2/analysis.jl
@@ -48,11 +48,7 @@ function guess_point_type(::HopfCont, ::Continuous, cache, opts)
     end
 
     if length(cache.eigvals) >= 3
-        # TODO: directly compute determinant and avoid computing
-        # eigenvalues all the time
-        d0 = prod(cache.prev_eigvals)
-        d1 = prod(cache.eigvals)
-        if d0 * d1 <= 0
+        if cache.prev_det * cache.det <= 0
             ev0 = sort_by_abs_real!(copy(cache.eigvals))[3]
             ev1 = sort_by_abs_real!(copy(cache.prev_eigvals))[3]
             if mean(imag.((ev0, ev1))) <= opts.atol

diff --git a/src/codim2/solver.jl b/src/codim2/solver.jl
@@ -82,6 +82,8 @@ mutable struct Codim2Cache{P, C <: ContinuationCache{P},
                            } <: BifurcationCache{P}
     super::C
     J::JType
+    det::Float64
+    prev_det::Float64
     eigvals::eType
     prev_eigvals::eType
     point_type::PointType
@@ -109,6 +111,7 @@ function Codim2Cache(super::C,
                      ) where {P, C <: ContinuationCache{P},
                               JType, eType}
     return Codim2Cache{P, C, JType, eType}(super, J,
+                                           NaN, NaN,
                                            eigvals, copy(eigvals),
                                            point_type,
                                            NaN, NaN,
@@ -158,6 +161,8 @@ end
 function analyze!(cache::Codim2Cache, opts)
     cnt = contkind(cache)
     cache.J = J = ds_jacobian(cache)
+    cache.prev_det = cache.det
+    cache.det = det(J)
     cache.prev_eigvals = cache.eigvals
     cache.eigvals = ds_eigvals(timekind(cache), J)
     set_quadratic_coefficient!(cnt, cache)