computes column norms outside of ei_qrfac()

3 files changed, 27 insertions, 36 deletions

diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
index 6269a3d89..b2e297741 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
@@ -226,15 +226,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
         return UserAksed;
     ++njev;
 
-    /* compute the qr factorization of the jacobian. */
-
-    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data(), wa2.data());
-
+    wa2 = fjac.colwise().blueNorm();
 
+    /* on the first iteration and if mode is 1, scale according */
+    /* to the norms of the columns of the initial jacobian. */
     if (iter == 1) {
-
-        /* on the first iteration and if mode is 1, scale according */
-        /* to the norms of the columns of the initial jacobian. */
         if (mode != 2)
             for (j = 0; j < n; ++j) {
                 diag[j] = wa2[j];
@@ -251,6 +247,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
             delta = parameters.factor;
     }
 
+    /* compute the qr factorization of the jacobian. */
+    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data());
+
     /* form (q transpose)*fvec and store in qtf. */
 
     qtf = fvec;
@@ -269,18 +268,16 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
     sing = false;
     for (j = 0; j < n; ++j) {
         l = j;
-        if (j)
-            for (i = 0; i < j; ++i) {
-                R[l] = fjac(i,j);
-                l = l + n - i -1;
-            }
+        for (i = 0; i < j; ++i) {
+            R[l] = fjac(i,j);
+            l = l + n - i -1;
+        }
         R[l] = wa1[j];
         if (wa1[j] == 0.)
             sing = true;
     }
 
     /* accumulate the orthogonal factor in fjac. */
-
     ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
     /* rescale if necessary. */
@@ -543,13 +540,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
         return UserAksed;
     nfev += std::min(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1, n);
 
-    /* compute the qr factorization of the jacobian. */
-
-    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data(), wa2.data());
+    wa2 = fjac.colwise().blueNorm();
 
     /* on the first iteration and if mode is 1, scale according */
     /* to the norms of the columns of the initial jacobian. */
-
     if (iter == 1) {
         if (mode != 2)
             for (j = 0; j < n; ++j) {
@@ -560,7 +554,6 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
 
         /* on the first iteration, calculate the norm of the scaled x */
         /* and initialize the step bound delta. */
-
         wa3 = diag.cwise() * x;
         xnorm = wa3.stableNorm();
         delta = parameters.factor * xnorm;
@@ -568,6 +561,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
             delta = parameters.factor;
     }
 
+    /* compute the qr factorization of the jacobian. */
+    ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), false, iwa, wa1.data());
+
     /* form (q transpose)*fvec and store in qtf. */
 
     qtf = fvec;
@@ -586,18 +582,16 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
     sing = false;
     for (j = 0; j < n; ++j) {
         l = j;
-        if (j)
-            for (i = 0; i < j; ++i) {
-                R[l] = fjac(i,j);
-                l = l + n - i -1;
-            }
+        for (i = 0; i < j; ++i) {
+            R[l] = fjac(i,j);
+            l = l + n - i -1;
+        }
         R[l] = wa1[j];
         if (wa1[j] == 0.)
             sing = true;
     }
 
     /* accumulate the orthogonal factor in fjac. */
-
     ei_qform<Scalar>(n, n, fjac.data(), fjac.rows(), wa1.data());
 
     /* rescale if necessary. */
diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
index 1a2a69561..c611ec595 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h
@@ -248,8 +248,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
 
     /* compute the qr factorization of the jacobian. */
 
-    ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data(), wa2.data());
-    ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
+    wa2 = fjac.colwise().blueNorm();
+    ei_qrfac<Scalar>(m, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data());
+    ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convention (1->n), convert it to c (0->n-1)
 
     /* on the first iteration and if mode is 1, scale according */
     /* to the norms of the columns of the initial jacobian. */
@@ -537,8 +538,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
     }
     if (sing) {
         ipvt.cwise()+=1;
-        ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data(), wa2.data());
-        ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
+        wa2 = fjac.colwise().blueNorm();
+        ei_qrfac<Scalar>(n, n, fjac.data(), fjac.rows(), true, ipvt.data(), wa1.data());
+        ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convention (1->n), convert it to c (0->n-1)
         for (j = 0; j < n; ++j) {
             if (fjac(j,j) != 0.) {
                 sum = 0.;
diff --git a/unsupported/Eigen/src/NonLinearOptimization/qrfac.h b/unsupported/Eigen/src/NonLinearOptimization/qrfac.h
index 481fe57d8..0c1ecf394 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/qrfac.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/qrfac.h
@@ -1,8 +1,7 @@
 
 template <typename Scalar>
 void ei_qrfac(int m, int n, Scalar *a, int
-        lda, int pivot, int *ipvt, Scalar *rdiag,
-        Scalar *acnorm)
+        lda, int pivot, int *ipvt, Scalar *rdiag)
 {
     /* System generated locals */
     int a_dim1, a_offset;
@@ -18,7 +17,6 @@ void ei_qrfac(int m, int n, Scalar *a, int
     Matrix< Scalar, Dynamic, 1 > wa(n+1);
 
     /* Parameter adjustments */
-    --acnorm;
     --rdiag;
     a_dim1 = lda;
     a_offset = 1 + a_dim1 * 1;
@@ -31,13 +29,10 @@ void ei_qrfac(int m, int n, Scalar *a, int
     /*     compute the initial column norms and initialize several arrays. */
 
     for (j = 1; j <= n; ++j) {
-        acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
-        rdiag[j] = acnorm[j];
+        rdiag[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
         wa[j] = rdiag[j];
-        if (pivot) {
+        if (pivot)
             ipvt[j] = j;
-        }
-        /* L10: */
     }
 
     /*     reduce a to r with householder transformations. */