From 67fb7fa9ff549a41172a9eee08d8ee7c8c03d405 Mon Sep 17 00:00:00 2001
From: Frank Dellaert <dellaert@gmail.com>
Date: Thu, 11 Mar 2010 15:04:31 +0000
Subject: [PATCH] in-place house vector (was about 10% of all mallocs in Urban)

---
 cpp/Matrix.cpp |  6 +++---
 cpp/Vector.cpp | 20 +++++++++++++-------
 cpp/Vector.h   |  3 +++
 3 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp
index 7d3d70388..3f60e4cad 100644
--- a/cpp/Matrix.cpp
+++ b/cpp/Matrix.cpp
@@ -641,9 +641,9 @@ inline void householder_manual(Matrix &A, size_t k) {
 		for(size_t r = j ; r < m; r++)
 			xjm(r-j) = A(r,j);
 
-		// calculate the Householder vector
-		double beta; Vector vjm;
-		boost::tie(beta,vjm) = house(xjm);
+		// calculate the Householder vector, in place
+		double beta = houseInPlace(xjm);
+		Vector& vjm = xjm;
 
 		// do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
 		householder_update(A, j, beta, vjm);
diff --git a/cpp/Vector.cpp b/cpp/Vector.cpp
index 98364ef18..2bcb6262c 100644
--- a/cpp/Vector.cpp
+++ b/cpp/Vector.cpp
@@ -327,18 +327,18 @@ namespace gtsam {
   }
 
   /* ************************************************************************* */
-  pair<double, Vector > house(const Vector &x) {
-	  const double x0 = x(0);
+  // imperative version, pass in x
+  double houseInPlace(Vector &v) {
+	  const double x0 = v(0);
 	  const double x02 = x0*x0;
 
 	  // old code - GSL verison was actually a bit slower
-	  const double sigma = inner_prod(trans(x),x) - x02;
-	  double beta = 0.0;
+	  const double sigma = inner_prod(v,v) - x02;
 
-	  Vector v(x); v(0) = 1.0;
+	  v(0) = 1.0;
 
 	  if( sigma == 0.0 )
-		  beta = 0.0;
+		  return 0.0;
 	  else {
 		  double mu = sqrt(x02 + sigma);
 		  if( x0 <= 0.0 )
@@ -347,9 +347,15 @@ namespace gtsam {
 			  v(0) = -sigma / (x0 + mu);
 
 		  const double v02 = v(0)*v(0);
-		  beta = 2.0 * v02 / (sigma + v02);
 		  v = v / v(0);
+		  return 2.0 * v02 / (sigma + v02);
 	  }
+  }
+
+  /* ************************************************************************* */
+  pair<double, Vector > house(const Vector &x) {
+  	Vector v(x);
+  	double beta = houseInPlace(v);
 	  return make_pair(beta, v);
   }
   
diff --git a/cpp/Vector.h b/cpp/Vector.h
index 8b6282a6b..23b354a98 100644
--- a/cpp/Vector.h
+++ b/cpp/Vector.h
@@ -250,6 +250,9 @@ Vector operator/(double s, const Vector& v);
  */
 std::pair<double,Vector> house(const Vector &x);
 
+/** beta = house(x) computes the HouseHolder vector in place */
+double houseInPlace(Vector &x);
+
 /**
  * Weighted Householder solution vector,
  * a.k.a., the pseudoinverse of the column