merge with tip

8 files changed, 134 insertions, 224 deletions

diff --git a/unsupported/Eigen/FFT b/unsupported/Eigen/FFT
index c8521bbf0..fc2efc1d6 100644
--- a/unsupported/Eigen/FFT
+++ b/unsupported/Eigen/FFT
@@ -44,11 +44,42 @@
   * The build-in implementation is based on kissfft. It is a small, free, and
   * reasonably efficient default.
   *
-  * Frontends are
+  * There are currently two frontends:
   *
   * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
-  * - MLK (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form
+  * - MLK (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
   *
+  * \section FFTDesign Design
+  *
+  * The following design decisions were made concerning scaling and
+  * half-spectrum for real FFT.
+  *
+  * The intent is to facilitate generic programming and ease migrating code
+  * from  Matlab/octave.
+  * We think the default behavior of Eigen/FFT should favor correctness and
+  * generality over speed. Of course, the caller should be able to "opt-out" from this
+  * behavior and get the speed increase if they want it.
+  *
+  * 1) %Scaling:
+  * Other libraries (FFTW,IMKL,KISSFFT)  do not perform scaling, so there
+  * is a constant gain incurred after the forward&inverse transforms , so 
+  * IFFT(FFT(x)) = Kx;  this is done to avoid a vector-by-value multiply.  
+  * The downside is that algorithms that worked correctly in Matlab/octave 
+  * don't behave the same way once implemented in C++.
+  *
+  * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x. 
+  *
+  * 2) Real FFT half-spectrum
+  * Other libraries use only half the frequency spectrum (plus one extra 
+  * sample for the Nyquist bin) for a real FFT, the other half is the 
+  * conjugate-symmetric of the first half.  This saves them a copy and some 
+  * memory.  The downside is the caller needs to have special logic for the 
+  * number of bins in complex vs real.
+  *
+  * How Eigen/FFT differs: The full spectrum is returned from the forward 
+  * transform.  This facilitates generic template programming by obviating 
+  * separate specializations for real vs complex.  On the inverse
+  * transform, only half the spectrum is actually used if the output type is real.
   */
  
 
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..5895fb578 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. */
@@ -319,7 +320,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
 
         /* determine the levenberg-marquardt parameter. */
 
-        ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1, wa2);
+        ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1);
 
         /* store the direction p and x + p. calculate the norm of p. */
 
@@ -432,8 +433,6 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
 {
     n = x.size();
     m = functor.values();
-    JacobianType fjac(m, n);
-    VectorXi ipvt;
 
     /* check the input parameters for errors. */
     if (n <= 0 || m < n || tol < 0.)
@@ -537,8 +536,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.;
@@ -603,7 +603,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
 
         /* determine the levenberg-marquardt parameter. */
 
-        ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1, wa2);
+        ei_lmpar<Scalar>(fjac, ipvt, diag, qtf, delta, par, wa1);
 
         /* store the direction p and x + p. calculate the norm of p. */
 
diff --git a/unsupported/Eigen/src/NonLinearOptimization/lmpar.h b/unsupported/Eigen/src/NonLinearOptimization/lmpar.h
index c62881c81..b723a7e0a 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/lmpar.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/lmpar.h
@@ -7,16 +7,14 @@ void ei_lmpar(
         const Matrix< Scalar, Dynamic, 1 >  &qtb,
         Scalar delta,
         Scalar &par,
-        Matrix< Scalar, Dynamic, 1 >  &x,
-        Matrix< Scalar, Dynamic, 1 >  &sdiag)
+        Matrix< Scalar, Dynamic, 1 >  &x)
 {
     /* Local variables */
-    int i, j, k, l;
+    int i, j, l;
     Scalar fp;
-    Scalar sum, parc, parl;
+    Scalar parc, parl;
     int iter;
     Scalar temp, paru;
-    int nsing;
     Scalar gnorm;
     Scalar dxnorm;
 
@@ -28,31 +26,28 @@ void ei_lmpar(
     assert(n==qtb.size());
     assert(n==x.size());
 
-    Matrix< Scalar, Dynamic, 1 >  wa1(n), wa2(n);
+    Matrix< Scalar, Dynamic, 1 >  wa1, wa2;
 
     /* compute and store in x the gauss-newton direction. if the */
     /* jacobian is rank-deficient, obtain a least squares solution. */
 
-    nsing = n-1;
+    int nsing = n-1;
+    wa1 = qtb;
     for (j = 0; j < n; ++j) {
-        wa1[j] = qtb[j];
         if (r(j,j) == 0. && nsing == n-1)
             nsing = j - 1;
         if (nsing < n-1) 
             wa1[j] = 0.;
     }
-    for (k = 0; k <= nsing; ++k) {
-        j = nsing - k;
+    for (j = nsing; j>=0; --j) {
         wa1[j] /= r(j,j);
         temp = wa1[j];
         for (i = 0; i < j ; ++i)
             wa1[i] -= r(i,j) * temp;
     }
 
-    for (j = 0; j < n; ++j) {
-        l = ipvt[j];
-        x[l] = wa1[j];
-    }
+    for (j = 0; j < n; ++j)
+        x[ipvt[j]] = wa1[j];
 
     /* initialize the iteration counter. */
     /* evaluate the function at the origin, and test */
@@ -77,8 +72,10 @@ void ei_lmpar(
             l = ipvt[j];
             wa1[j] = diag[l] * (wa2[l] / dxnorm);
         }
+        // it's actually a triangularView.solveInplace(), though in a weird
+        // way:
         for (j = 0; j < n; ++j) {
-            sum = 0.;
+            Scalar sum = 0.;
             for (i = 0; i < j; ++i) 
                 sum += r(i,j) * wa1[i];
             wa1[j] = (wa1[j] - sum) / r(j,j);
@@ -89,13 +86,9 @@ void ei_lmpar(
 
     /* calculate an upper bound, paru, for the zero of the function. */
 
-    for (j = 0; j < n; ++j) {
-        sum = 0.;
-        for (i = 0; i <= j; ++i)
-            sum += r(i,j) * qtb[i];
-        l = ipvt[j];
-        wa1[j] = sum / diag[l];
-    }
+    for (j = 0; j < n; ++j)
+        wa1[j] = r.col(j).start(j+1).dot(qtb.start(j+1)) / diag[ipvt[j]];
+
     gnorm = wa1.stableNorm();
     paru = gnorm / delta;
     if (paru == 0.)
@@ -119,12 +112,10 @@ void ei_lmpar(
         if (par == 0.)
             par = std::max(dwarf,Scalar(.001) * paru); /* Computing MAX */
 
-        temp = ei_sqrt(par);
-        wa1 = temp * diag;
+        wa1 = ei_sqrt(par)* diag;
 
-        ipvt.cwise()+=1; // qrsolv() expects the fortran convention (as qrfac provides)
-        ei_qrsolv<Scalar>(n, r.data(), r.rows(), ipvt.data(), wa1.data(), qtb.data(), x.data(), sdiag.data(), wa2.data());
-        ipvt.cwise()-=1;
+        Matrix< Scalar, Dynamic, 1 > sdiag(n);
+        ei_qrsolv<Scalar>(r, ipvt, wa1, qtb, x, sdiag);
 
         wa2 = diag.cwise() * x;
         dxnorm = wa2.blueNorm();
@@ -143,7 +134,6 @@ void ei_lmpar(
         for (j = 0; j < n; ++j) {
             l = ipvt[j];
             wa1[j] = diag[l] * (wa2[l] / dxnorm);
-            /* L180: */
         }
         for (j = 0; j < n; ++j) {
             wa1[j] /= sdiag[j];
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. */
diff --git a/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h b/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h
index 57884870c..1ee21d5ed 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/qrsolv.h
@@ -1,166 +1,92 @@
 
-    template <typename Scalar>
-void ei_qrsolv(int n, Scalar *r__, int ldr, 
-        const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x, 
-        Scalar *sdiag, Scalar *wa)
-{
-    /* System generated locals */
-    int r_dim1, r_offset;
+template <typename Scalar>
+void ei_qrsolv(
+        Matrix< Scalar, Dynamic, Dynamic > &r,
+        VectorXi &ipvt, // TODO : const once ipvt mess fixed
+        const Matrix< Scalar, Dynamic, 1 >  &diag,
+        const Matrix< Scalar, Dynamic, 1 >  &qtb,
+        Matrix< Scalar, Dynamic, 1 >  &x,
+        Matrix< Scalar, Dynamic, 1 >  &sdiag)
 
+{
     /* Local variables */
-    int i, j, k, l, jp1, kp1;
-    Scalar tan__, cos__, sin__, sum, temp, cotan;
-    int nsing;
-    Scalar qtbpj;
-
-    /* Parameter adjustments */
-    --wa;
-    --sdiag;
-    --x;
-    --qtb;
-    --diag;
-    --ipvt;
-    r_dim1 = ldr;
-    r_offset = 1 + r_dim1 * 1;
-    r__ -= r_offset;
+    int i, j, k, l;
+    Scalar sum, temp;
+    int n = r.cols();
+    Matrix< Scalar, Dynamic, 1 >  wa(n);
 
     /* Function Body */
 
     /*     copy r and (q transpose)*b to preserve input and initialize s. */
     /*     in particular, save the diagonal elements of r in x. */
 
-    for (j = 1; j <= n; ++j) {
-        for (i = j; i <= n; ++i) {
-            r__[i + j * r_dim1] = r__[j + i * r_dim1];
-            /* L10: */
-        }
-        x[j] = r__[j + j * r_dim1];
-        wa[j] = qtb[j];
-        /* L20: */
-    }
+    x = r.diagonal();
+    wa = qtb;
 
-    /*     eliminate the diagonal matrix d using a givens rotation. */
+    for (j = 0; j < n; ++j)
+        for (i = j+1; i < n; ++i)
+            r(i,j) = r(j,i);
 
-    for (j = 1; j <= n; ++j) {
+    /*     eliminate the diagonal matrix d using a givens rotation. */
+    for (j = 0; j < n; ++j) {
 
         /*        prepare the row of d to be eliminated, locating the */
         /*        diagonal element using p from the qr factorization. */
 
         l = ipvt[j];
-        if (diag[l] == 0.) {
-            goto L90;
-        }
-        for (k = j; k <= n; ++k) {
-            sdiag[k] = 0.;
-            /* L30: */
-        }
+        if (diag[l] == 0.)
+            break;
+        sdiag.segment(j,n-j).setZero();
         sdiag[j] = diag[l];
 
         /*        the transformations to eliminate the row of d */
         /*        modify only a single element of (q transpose)*b */
         /*        beyond the first n, which is initially zero. */
 
-        qtbpj = 0.;
-        for (k = j; k <= n; ++k) {
-
+        Scalar qtbpj = 0.;
+        for (k = j; k < n; ++k) {
             /*           determine a givens rotation which eliminates the */
             /*           appropriate element in the current row of d. */
-
-            if (sdiag[k] == 0.)
-                goto L70;
-            if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k]))
-                goto L40;
-            cotan = r__[k + k * r_dim1] / sdiag[k];
-            /* Computing 2nd power */
-            sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
-            cos__ = sin__ * cotan;
-            goto L50;
-L40:
-            tan__ = sdiag[k] / r__[k + k * r_dim1];
-            /* Computing 2nd power */
-            cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
-            sin__ = cos__ * tan__;
-L50:
+            PlanarRotation<Scalar> givens;
+            givens.makeGivens(-r(k,k), sdiag[k]);
 
             /*           compute the modified diagonal element of r and */
             /*           the modified element of ((q transpose)*b,0). */
 
-            r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[
-                k];
-            temp = cos__ * wa[k] + sin__ * qtbpj;
-            qtbpj = -sin__ * wa[k] + cos__ * qtbpj;
+            r(k,k) = givens.c() * r(k,k) + givens.s() * sdiag[k];
+            temp = givens.c() * wa[k] + givens.s() * qtbpj;
+            qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj;
             wa[k] = temp;
 
             /*           accumulate the tranformation in the row of s. */
-
-            kp1 = k + 1;
-            if (n < kp1) {
-                goto L70;
+            for (i = k+1; i<n; ++i) {
+                temp = givens.c() * r(i,k) + givens.s() * sdiag[i];
+                sdiag[i] = -givens.s() * r(i,k) + givens.c() * sdiag[i];
+                r(i,k) = temp;
             }
-            for (i = kp1; i <= n; ++i) {
-                temp = cos__ * r__[i + k * r_dim1] + sin__ * sdiag[i];
-                sdiag[i] = -sin__ * r__[i + k * r_dim1] + cos__ * sdiag[
-                    i];
-                r__[i + k * r_dim1] = temp;
-                /* L60: */
-            }
-L70:
-            /* L80: */
-            ;
         }
-L90:
-
-        /*        store the diagonal element of s and restore */
-        /*        the corresponding diagonal element of r. */
-
-        sdiag[j] = r__[j + j * r_dim1];
-        r__[j + j * r_dim1] = x[j];
-        /* L100: */
     }
 
+    // restore
+    sdiag = r.diagonal();
+    r.diagonal() = x;
+
     /*     solve the triangular system for z. if the system is */
     /*     singular, then obtain a least squares solution. */
 
-    nsing = n;
-    for (j = 1; j <= n; ++j) {
-        if (sdiag[j] == 0. && nsing == n) {
-            nsing = j - 1;
-        }
-        if (nsing < n) {
-            wa[j] = 0.;
-        }
-        /* L110: */
-    }
-    if (nsing < 1) {
-        goto L150;
-    }
-    for (k = 1; k <= nsing; ++k) {
-        j = nsing - k + 1;
+    int nsing;
+    for (nsing=0; nsing<n && sdiag[nsing]!=0; nsing++);
+    wa.segment(nsing,n-nsing).setZero();
+    nsing--; // nsing is the last nonsingular index
+
+    for (j = nsing; j>=0; j--) {
         sum = 0.;
-        jp1 = j + 1;
-        if (nsing < jp1) {
-            goto L130;
-        }
-        for (i = jp1; i <= nsing; ++i) {
-            sum += r__[i + j * r_dim1] * wa[i];
-            /* L120: */
-        }
-L130:
+        for (i = j+1; i <= nsing; ++i)
+            sum += r(i,j) * wa[i];
         wa[j] = (wa[j] - sum) / sdiag[j];
-        /* L140: */
     }
-L150:
 
     /*     permute the components of z back to components of x. */
-
-    for (j = 1; j <= n; ++j) {
-        l = ipvt[j];
-        x[l] = wa[j];
-        /* L160: */
-    }
-    return;
-
-    /*     last card of subroutine qrsolv. */
-
-} /* qrsolv_ */
+    for (j = 0; j < n; ++j) x[ipvt[j]] = wa[j];
+}
 
diff --git a/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h b/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
index a06072a54..70a6d30c3 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
@@ -18,28 +18,18 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
     a -= a_offset;
 
     /* Function Body */
-
-    /*     apply the first set of givens rotations to a. */
-
     nm1 = n - 1;
-    if (nm1 < 1) {
-        /* goto L50; */
+    if (nm1 < 1)
         return;
-    }
+
+    /*     apply the first set of givens rotations to a. */
     for (nmj = 1; nmj <= nm1; ++nmj) {
         j = n - nmj;
         if (ei_abs(v[j]) > 1.) {
             cos__ = 1. / v[j];
-        }
-        if (ei_abs(v[j]) > 1.) {
-            /* Computing 2nd power */
             sin__ = ei_sqrt(1. - ei_abs2(cos__));
-        }
-        if (ei_abs(v[j]) <= 1.) {
+        } else  {
             sin__ = v[j];
-        }
-        if (ei_abs(v[j]) <= 1.) {
-            /* Computing 2nd power */
             cos__ = ei_sqrt(1. - ei_abs2(sin__));
         }
         for (i = 1; i <= m; ++i) {
@@ -47,26 +37,15 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
             a[i + n * a_dim1] = sin__ * a[i + j * a_dim1] + cos__ * a[
                 i + n * a_dim1];
             a[i + j * a_dim1] = temp;
-            /* L10: */
         }
-        /* L20: */
     }
-
     /*     apply the second set of givens rotations to a. */
-
     for (j = 1; j <= nm1; ++j) {
         if (ei_abs(w[j]) > 1.) {
             cos__ = 1. / w[j];
-        }
-        if (ei_abs(w[j]) > 1.) {
-            /* Computing 2nd power */
             sin__ = ei_sqrt(1. - ei_abs2(cos__));
-        }
-        if (ei_abs(w[j]) <= 1.) {
+        } else  {
             sin__ = w[j];
-        }
-        if (ei_abs(w[j]) <= 1.) {
-            /* Computing 2nd power */
             cos__ = ei_sqrt(1. - ei_abs2(sin__));
         }
         for (i = 1; i <= m; ++i) {
@@ -74,14 +53,8 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
             a[i + n * a_dim1] = -sin__ * a[i + j * a_dim1] + cos__ * a[
                 i + n * a_dim1];
             a[i + j * a_dim1] = temp;
-            /* L30: */
         }
-        /* L40: */
     }
-    /* L50: */
     return;
-
-    /*     last card of subroutine r1mpyq. */
-
 } /* r1mpyq_ */
 
diff --git a/unsupported/doc/Doxyfile.in b/unsupported/doc/Doxyfile.in
index a9c7d6004..e55d53f7e 100644
--- a/unsupported/doc/Doxyfile.in
+++ b/unsupported/doc/Doxyfile.in
@@ -601,6 +601,7 @@ EXCLUDE_PATTERNS       = CMake* \
                          *.txt \
                          *.sh \
                          *.diff \
+                         *.orig \
                          diff
                          *~