From c04a93df31c32912c81704535fc10849ac11f076 Mon Sep 17 00:00:00 2001
From: Thomas Capricelli <orzel@freehackers.org>
Date: Tue, 26 Jan 2010 14:08:52 +0100
Subject: [PATCH] clean r1mpyq: remove fortran leftovers

---
 .../HybridNonLinearSolver.h                   |  8 ++--
 .../Eigen/src/NonLinearOptimization/r1mpyq.h  | 40 +++++++------------
 2 files changed, 18 insertions(+), 30 deletions(-)

diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
index cbf89b1c69..e1d10bdccc 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h
@@ -360,8 +360,8 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
 
         /* compute the qr factorization of the updated jacobian. */
         ei_r1updt<Scalar>(n, n, R, wa1.data(), wa2.data(), wa3.data(), &sing);
-        ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
-        ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+        ei_r1mpyq<Scalar>(n, n, fjac.data(), wa2.data(), wa3.data());
+        ei_r1mpyq<Scalar>(1, n, qtf.data(), wa2.data(), wa3.data());
 
         jeval = false;
     }
@@ -609,8 +609,8 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
 
         /* compute the qr factorization of the updated jacobian. */
         ei_r1updt<Scalar>(n, n, R, wa1.data(), wa2.data(), wa3.data(), &sing);
-        ei_r1mpyq<Scalar>(n, n, fjac.data(), fjac.rows(), wa2.data(), wa3.data());
-        ei_r1mpyq<Scalar>(1, n, qtf.data(), 1, wa2.data(), wa3.data());
+        ei_r1mpyq<Scalar>(n, n, fjac.data(), wa2.data(), wa3.data());
+        ei_r1mpyq<Scalar>(1, n, qtf.data(), wa2.data(), wa3.data());
 
         jeval = false;
     }
diff --git a/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h b/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
index 9134dc0b63..ac1feeffc1 100644
--- a/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
+++ b/unsupported/Eigen/src/NonLinearOptimization/r1mpyq.h
@@ -1,29 +1,19 @@
 
+// TODO : move this to GivensQR once there's such a thing in Eigen
+
 template <typename Scalar>
-void ei_r1mpyq(int m, int n, Scalar *a, int
-        lda, const Scalar *v, const Scalar *w)
+void ei_r1mpyq(int m, int n, Scalar *a, const Scalar *v, const Scalar *w)
 {
-    /* System generated locals */
-    int a_dim1, a_offset;
-
     /* Local variables */
-    int i, j, nmj;
+    int i, j;
     Scalar cos__=0., sin__=0., temp;
 
-    /* Parameter adjustments */
-    --w;
-    --v;
-    a_dim1 = lda;
-    a_offset = 1 + a_dim1 * 1;
-    a -= a_offset;
-
     /* Function Body */
     if (n<=1)
         return;
 
     /*     apply the first set of givens rotations to a. */
-    for (nmj = 1; nmj <= n-1; ++nmj) {
-        j = n - nmj;
+    for (j = n-2; j>=0; --j) {
         if (ei_abs(v[j]) > 1.) {
             cos__ = 1. / v[j];
             sin__ = ei_sqrt(1. - ei_abs2(cos__));
@@ -31,15 +21,14 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
             sin__ = v[j];
             cos__ = ei_sqrt(1. - ei_abs2(sin__));
         }
-        for (i = 1; i <= m; ++i) {
-            temp = cos__ * a[i + j * a_dim1] - sin__ * a[i + n * a_dim1];
-            a[i + n * a_dim1] = sin__ * a[i + j * a_dim1] + cos__ * a[
-                i + n * a_dim1];
-            a[i + j * a_dim1] = temp;
+        for (i = 0; i<m; ++i) {
+            temp = cos__ * a[i+m*j] - sin__ * a[i+m*(n-1)];
+            a[i+m*(n-1)] = sin__ * a[i+m*j] + cos__ * a[i+m*(n-1)];
+            a[i+m*j] = temp;
         }
     }
     /*     apply the second set of givens rotations to a. */
-    for (j = 1; j <= n-1; ++j) {
+    for (j = 0; j<n-1; ++j) {
         if (ei_abs(w[j]) > 1.) {
             cos__ = 1. / w[j];
             sin__ = ei_sqrt(1. - ei_abs2(cos__));
@@ -47,11 +36,10 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
             sin__ = w[j];
             cos__ = ei_sqrt(1. - ei_abs2(sin__));
         }
-        for (i = 1; i <= m; ++i) {
-            temp = cos__ * a[i + j * a_dim1] + sin__ * a[i + n * a_dim1];
-            a[i + n * a_dim1] = -sin__ * a[i + j * a_dim1] + cos__ * a[
-                i + n * a_dim1];
-            a[i + j * a_dim1] = temp;
+        for (i = 0; i<m; ++i) {
+            temp = cos__ * a[i+m*j] + sin__ * a[i+m*(n-1)];
+            a[i+m*(n-1)] = -sin__ * a[i+m*j] + cos__ * a[i+m*(n-1)];
+            a[i+m*j] = temp;
         }
     }
     return;