1 files changed, 107 insertions, 0 deletions
diff --git a/eigen/unsupported/Eigen/src/NonLinearOptimization/dogleg.h b/eigen/unsupported/Eigen/src/NonLinearOptimization/dogleg.h
new file mode 100644
index 0000000..80c5d27
--- /dev/null
+++ b/eigen/unsupported/Eigen/src/NonLinearOptimization/dogleg.h
@@ -0,0 +1,107 @@
+namespace Eigen { 
+
+namespace internal {
+
+template <typename Scalar>
+void dogleg(
+        const Matrix< Scalar, Dynamic, Dynamic >  &qrfac,
+        const Matrix< Scalar, Dynamic, 1 >  &diag,
+        const Matrix< Scalar, Dynamic, 1 >  &qtb,
+        Scalar delta,
+        Matrix< Scalar, Dynamic, 1 >  &x)
+{
+    using std::abs;
+    using std::sqrt;
+    
+    typedef DenseIndex Index;
+
+    /* Local variables */
+    Index i, j;
+    Scalar sum, temp, alpha, bnorm;
+    Scalar gnorm, qnorm;
+    Scalar sgnorm;
+
+    /* Function Body */
+    const Scalar epsmch = NumTraits<Scalar>::epsilon();
+    const Index n = qrfac.cols();
+    eigen_assert(n==qtb.size());
+    eigen_assert(n==x.size());
+    eigen_assert(n==diag.size());
+    Matrix< Scalar, Dynamic, 1 >  wa1(n), wa2(n);
+
+    /* first, calculate the gauss-newton direction. */
+    for (j = n-1; j >=0; --j) {
+        temp = qrfac(j,j);
+        if (temp == 0.) {
+            temp = epsmch * qrfac.col(j).head(j+1).maxCoeff();
+            if (temp == 0.)
+                temp = epsmch;
+        }
+        if (j==n-1)
+            x[j] = qtb[j] / temp;
+        else
+            x[j] = (qtb[j] - qrfac.row(j).tail(n-j-1).dot(x.tail(n-j-1))) / temp;
+    }
+
+    /* test whether the gauss-newton direction is acceptable. */
+    qnorm = diag.cwiseProduct(x).stableNorm();
+    if (qnorm <= delta)
+        return;
+
+    // TODO : this path is not tested by Eigen unit tests
+
+    /* the gauss-newton direction is not acceptable. */
+    /* next, calculate the scaled gradient direction. */
+
+    wa1.fill(0.);
+    for (j = 0; j < n; ++j) {
+        wa1.tail(n-j) += qrfac.row(j).tail(n-j) * qtb[j];
+        wa1[j] /= diag[j];
+    }
+
+    /* calculate the norm of the scaled gradient and test for */
+    /* the special case in which the scaled gradient is zero. */
+    gnorm = wa1.stableNorm();
+    sgnorm = 0.;
+    alpha = delta / qnorm;
+    if (gnorm == 0.)
+        goto algo_end;
+
+    /* calculate the point along the scaled gradient */
+    /* at which the quadratic is minimized. */
+    wa1.array() /= (diag*gnorm).array();
+    // TODO : once unit tests cover this part,:
+    // wa2 = qrfac.template triangularView<Upper>() * wa1;
+    for (j = 0; j < n; ++j) {
+        sum = 0.;
+        for (i = j; i < n; ++i) {
+            sum += qrfac(j,i) * wa1[i];
+        }
+        wa2[j] = sum;
+    }
+    temp = wa2.stableNorm();
+    sgnorm = gnorm / temp / temp;
+
+    /* test whether the scaled gradient direction is acceptable. */
+    alpha = 0.;
+    if (sgnorm >= delta)
+        goto algo_end;
+
+    /* the scaled gradient direction is not acceptable. */
+    /* finally, calculate the point along the dogleg */
+    /* at which the quadratic is minimized. */
+    bnorm = qtb.stableNorm();
+    temp = bnorm / gnorm * (bnorm / qnorm) * (sgnorm / delta);
+    temp = temp - delta / qnorm * numext::abs2(sgnorm / delta) + sqrt(numext::abs2(temp - delta / qnorm) + (1.-numext::abs2(delta / qnorm)) * (1.-numext::abs2(sgnorm / delta)));
+    alpha = delta / qnorm * (1. - numext::abs2(sgnorm / delta)) / temp;
+algo_end:
+
+    /* form appropriate convex combination of the gauss-newton */
+    /* direction and the scaled gradient direction. */
+    temp = (1.-alpha) * (std::min)(sgnorm,delta);
+    x = temp * wa1 + alpha * x;
+}
+
+} // end namespace internal
+
+} // end namespace Eigen