From 547a912d5e9a8d9228fb4eb5b827199ebfff6e10 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien.villemot@ens.fr>
Date: Fri, 19 Feb 2010 16:55:03 +0100
Subject: [PATCH] Estimation DLL, matrix library: * fixed various bugs * added
 utility for computing the infinite norm * changed the direction of
 mat::transpose(A, B) for consistency with GSL

---
 mex/sources/estimation/DecisionRules.cc |  3 +-
 mex/sources/estimation/libmat/Matrix.hh | 42 +++++++++++++++++++------
 mex/sources/estimation/libmat/Vector.hh | 22 +++++++++++--
 3 files changed, 54 insertions(+), 13 deletions(-)

diff --git a/mex/sources/estimation/DecisionRules.cc b/mex/sources/estimation/DecisionRules.cc
index 72f9afc38..36d5d1903 100644
--- a/mex/sources/estimation/DecisionRules.cc
+++ b/mex/sources/estimation/DecisionRules.cc
@@ -116,7 +116,8 @@ DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (
   for (size_t j = 0; j < n_fwrd; j++)
     mat::col_copy(A, n_back_mixed + zeta_fwrd_mixed[pi_fwrd[j]], n_static, n - n_static,
                   E, n_back_mixed + pi_fwrd[j], 0);
-  mat::negate(MatrixView(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed)); // Here we take the opposite of some of the zeros initialized in the constructor, but it is not a problem
+  MatrixView E_tmp(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed);
+  mat::negate(E_tmp); // Here we take the opposite of some of the zeros initialized in the constructor, but it is not a problem
 
   // Perform the generalized Schur
   size_t sdim;
diff --git a/mex/sources/estimation/libmat/Matrix.hh b/mex/sources/estimation/libmat/Matrix.hh
index d0571bd23..2175b8b53 100644
--- a/mex/sources/estimation/libmat/Matrix.hh
+++ b/mex/sources/estimation/libmat/Matrix.hh
@@ -27,6 +27,7 @@
 #include <cstdlib>
 #include <cassert>
 #include <cstring>
+#include <cmath>
 
 #include "Vector.hh"
 
@@ -285,28 +286,29 @@ namespace mat
         std::swap(M(i,j), M(j,i));
   }
 
-  //! Computes M2 = M1' (even for rectangular matrices)
+  //! Computes M1 = M2' (even for rectangular matrices)
   template<class Mat1, class Mat2>
   inline void
-  transpose(Mat1 &M1, Mat2 &M2)
+  transpose(Mat1 &M1, const Mat2 &M2)
   {
     assert(M1.getRows() == M2.getCols() && M1.getCols() == M2.getRows());
     for (size_t i = 0; i < M1.getRows(); i++)
       for (size_t j = 0; j < M1.getCols(); j++)
-        M2(j, i) = M1(i, j);
+        M1(i, j) = M2(j, i);
   }
 
   //! Computes m1 = m1 + m2
   template<class Mat1, class Mat2>
   void
-  add(Mat1 m1, const Mat2 m2)
+  add(Mat1 &m1, const Mat2 &m2)
   {
     assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
     double *p1 = m1.getData();
     const double *p2 = m2.getData();
     while (p1 < m1.getData() + m1.getCols() * m1.getLd())
       {
-        double *pp1 = p1, *pp2 = p2;
+        double *pp1 = p1;
+        const double *pp2 = p2;
         while (pp1 < p1 + m1.getRows())
           *pp1++ += *pp2++;
 
@@ -315,17 +317,18 @@ namespace mat
       }
   }
 
-  //! Computes m1 = m1  m2
+  //! Computes m1 = m1 - m2
   template<class Mat1, class Mat2>
   void
-  sub(Mat1 m1, const Mat2 m2)
+  sub(Mat1 &m1, const Mat2 &m2)
   {
     assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
     double *p1 = m1.getData();
     const double *p2 = m2.getData();
     while (p1 < m1.getData() + m1.getCols() * m1.getLd())
       {
-        double *pp1 = p1, *pp2 = p2;
+        double *pp1 = p1;
+        const double *pp2 = p2;
         while (pp1 < p1 + m1.getRows())
           *pp1++ -= *pp2++;
 
@@ -337,7 +340,7 @@ namespace mat
   //! Does m = -m
   template<class Mat>
   void
-  negate(Mat m)
+  negate(Mat &m)
   {
     double *p = m.getData();
     while (p < m.getData() + m.getCols() * m.getLd())
@@ -353,6 +356,27 @@ namespace mat
       }
   }
 
+  // Computes the infinite norm of a matrix
+  template<class Mat>
+  double
+  nrminf(const Mat &m)
+  {
+    double nrm = 0;
+    const double *p = m.getData();
+    while (p < m.getData() + m.getCols() * m.getLd())
+      {
+        const double *pp = p;
+        while (pp < p + m.getRows())
+          {
+            if (fabs(*pp) > nrm)
+              nrm = fabs(*pp);
+            pp++;
+          }
+
+        p += m.getLd();
+      }
+    return nrm;
+  }
 } // End of namespace
 
 #endif
diff --git a/mex/sources/estimation/libmat/Vector.hh b/mex/sources/estimation/libmat/Vector.hh
index 43efcf0bb..23d96997d 100644
--- a/mex/sources/estimation/libmat/Vector.hh
+++ b/mex/sources/estimation/libmat/Vector.hh
@@ -25,6 +25,7 @@
 
 #include <cstdlib>
 #include <cassert>
+#include <cmath>
 
 /*
   This header defines three vector classes, which implement a "vector concept"
@@ -159,7 +160,7 @@ namespace vec
   //! Computes v1 = v1 + v2
   template<class Vec1, class Vec2>
   void
-  add(Vec1 v1, const Vec2 v2)
+  add(Vec1 &v1, const Vec2 &v2)
   {
     assert(v1.getSize() == v2.getSize());
     double *p1 = v1.getData();
@@ -175,7 +176,7 @@ namespace vec
   //! Computes v1 = v1 - v2
   template<class Vec1, class Vec2>
   void
-  sub(Vec1 v1, const Vec2 v2)
+  sub(Vec1 &v1, const Vec2 &v2)
   {
     assert(v1.getSize() == v2.getSize());
     double *p1 = v1.getData();
@@ -191,7 +192,7 @@ namespace vec
   //! Does v = -v
   template<class Vec>
   void
-  negate(Vec v)
+  negate(Vec &v)
   {
     double *p = v.getData();
     while (p < v.getData() + v.getSize() * v.getStride())
@@ -201,6 +202,21 @@ namespace vec
       }
   }
 
+  // Computes the infinite norm of a vector
+  template<class Vec>
+  double
+  nrminf(const Vec &v)
+  {
+    double nrm = 0;
+    const double *p = v.getData();
+    while (p < v.getData() + v.getSize() * v.getStride())
+      {
+        if (fabs(*p) > nrm)
+          nrm = fabs(*p);
+        p += v.getStride();
+      }
+    return nrm;
+  }
 } // End of namespace
 
 #endif