From 3e91e51b6d0c6e028f4e9d7a2bd38ca150f169ab Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Thu, 30 Apr 2020 09:08:35 +0200
Subject: [PATCH] Root Finding: Brent behavior more consistent between managed
 and native providers

#690
---
 src/Numerics/RootFinding/Broyden.cs | 66 ++++++++++++++++-------------
 1 file changed, 37 insertions(+), 29 deletions(-)

diff --git a/src/Numerics/RootFinding/Broyden.cs b/src/Numerics/RootFinding/Broyden.cs
index 84b1d464..e4cd2a83 100644
--- a/src/Numerics/RootFinding/Broyden.cs
+++ b/src/Numerics/RootFinding/Broyden.cs
@@ -82,42 +82,50 @@ namespace MathNet.Numerics.RootFinding
 
             Matrix<double> B = CalculateApproximateJacobian(f, initialGuess, y0, jacobianStepSize);
 
-            for (int i = 0; i <= maxIterations; i++)
+            try
             {
-                var dx = (DenseVector) (-B.LU().Solve(y));
-                var xnew = x + dx;
-                var ynew = new DenseVector(f(xnew.Values));
-                double gnew = ynew.L2Norm();
-
-                if (gnew > g)
+                for (int i = 0; i <= maxIterations; i++)
                 {
-                    double g2 = g*g;
-                    double scale = g2/(g2 + gnew*gnew);
-                    if (scale == 0.0)
+                    var dx = (DenseVector) (-B.LU().Solve(y));
+                    var xnew = x + dx;
+                    var ynew = new DenseVector(f(xnew.Values));
+                    double gnew = ynew.L2Norm();
+
+                    if (gnew > g)
                     {
-                        scale = 1.0e-4;
+                        double g2 = g*g;
+                        double scale = g2/(g2 + gnew*gnew);
+                        if (scale == 0.0)
+                        {
+                            scale = 1.0e-4;
+                        }
+
+                        dx = scale*dx;
+                        xnew = x + dx;
+                        ynew = new DenseVector(f(xnew.Values));
+                        gnew = ynew.L2Norm();
                     }
 
-                    dx = scale*dx;
-                    xnew = x + dx;
-                    ynew = new DenseVector(f(xnew.Values));
-                    gnew = ynew.L2Norm();
-                }
-
-                if (gnew < accuracy)
-                {
-                    root = xnew.Values;
-                    return true;
-                }
+                    if (gnew < accuracy)
+                    {
+                        root = xnew.Values;
+                        return true;
+                    }
 
-                // update Jacobian B
-                DenseVector dF = ynew - y;
-                Matrix<double> dB = (dF - B.Multiply(dx)).ToColumnMatrix() * dx.Multiply(1.0 / Math.Pow(dx.L2Norm(), 2)).ToRowMatrix();
-                B = B + dB;
+                    // update Jacobian B
+                    DenseVector dF = ynew - y;
+                    Matrix<double> dB = (dF - B.Multiply(dx)).ToColumnMatrix() * dx.Multiply(1.0 / Math.Pow(dx.L2Norm(), 2)).ToRowMatrix();
+                    B = B + dB;
 
-                x = xnew;
-                y = ynew;
-                g = gnew;
+                    x = xnew;
+                    y = ynew;
+                    g = gnew;
+                }
+            }
+            catch (InvalidParameterException)
+            {
+                root = null;
+                return false;
             }
 
             root = null;