Added in-place insertion functions to Matrix, as well as a #define flag to enable or disable GSL code without a reconfigure.

cpp/Matrix.cpp

@@ -28,6 +28,9 @@
 #include "Vector.h"
 #include "svdcmp.h"
 
+// use for switching quickly between GSL and nonGSL versions without reconfigure
+//#define REVERTGSL
+
 
 using namespace std;
 namespace ublas = boost::numeric::ublas;
@@ -199,13 +202,6 @@ Vector row_(const Matrix& A, size_t i) {
 
 	const double * Aptr = A.data().begin() + A.size2() * i;
 	return Vector_(A.size2(), Aptr);
-
-	// TODO: improve this
-//	size_t n = A.size2();
-//	Vector a(n);
-//	for (size_t j=0; j<n; ++j)
-//		a(j) = A(i,j);
-//	return a;
 }
 
 /* ************************************************************************* */
@@ -235,6 +231,33 @@ Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2) {
   return B;
 }
 
+/* ************************************************************************* */
+void insertSub(Matrix& big, const Matrix& small, size_t i, size_t j) {
+	// direct pointer method
+	size_t ib = big.size1(), jb = big.size2(),
+		   is = small.size1(), js = small.size2();
+
+	// pointer to start of window in big
+	double * bigptr = big.data().begin() + i*jb + j;
+	const double * smallptr = small.data().begin();
+	for (size_t row=0; row<is; ++row)
+		memcpy(bigptr+row*jb, smallptr+row*js, js*sizeof(double));
+}
+
+/* ************************************************************************* */
+void insertColumn(Matrix& A, const Vector& col, size_t j) {
+	ublas::matrix_column<Matrix> colproxy(A, j);
+	colproxy = col;
+}
+
+/* ************************************************************************* */
+void insertColumn(Matrix& A, const Vector& col, size_t i, size_t j) {
+	ublas::matrix_column<Matrix> colproxy(A, j);
+	ublas::vector_range<ublas::matrix_column<Matrix> > colsubproxy(colproxy,
+			ublas::range (i, i+col.size()));
+	colsubproxy = col;
+}
+
 /* ************************************************************************* */
 void solve(Matrix& A, Matrix& B)
 {
@@ -310,21 +333,8 @@ pair<Matrix,Matrix> qr(const Matrix& A) {
  * on a number of different matrices for which all columns change.
  */
 /* ************************************************************************* */
-void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
+inline void householder_update_manual(Matrix &A, int j, double beta, const Vector& vjm) {
-
 	const size_t m = A.size1(), n = A.size2();
-
-#ifdef GSL
-  // use GSL version
-  gsl_vector_const_view v = gsl_vector_const_view_array(vjm.data().begin(), m-j);
-  gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
-  gsl_matrix_view Ag_view = gsl_matrix_submatrix (&(Ag.matrix), j, 0, m-j, n);
-  gsl_linalg_householder_hm (beta, &(v.vector), &(Ag_view.matrix));
-
-#else
-  // elegant but slow: A -= Matrix(outer_prod(v,beta*trans(A)*v));
-  // faster code below
-
 	// w = beta*transpose(A(j:m,:))*v(j:m)
 	Vector w(n);
 	for( size_t c = 0; c < n; c++) {
@@ -345,6 +355,23 @@ void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
 			// A(r,c) -= vjm(r-j) * wjn(c-j);
 			(*a) -= v[s] * wc;
 	}
+}
+
+void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
+#ifdef GSL
+#ifndef REVERTGSL
+	const size_t m = A.size1(), n = A.size2();
+	// use GSL version
+	gsl_vector_const_view v = gsl_vector_const_view_array(vjm.data().begin(), m-j);
+	gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
+	gsl_matrix_view Ag_view = gsl_matrix_submatrix (&(Ag.matrix), j, 0, m-j, n);
+	gsl_linalg_householder_hm (beta, &(v.vector), &(Ag_view.matrix));
+#else
+	householder_update_manual(A,j,beta,vjm);
+#endif
+
+#else
+	householder_update_manual(A,j,beta,vjm);
 #endif
 }
 
@@ -352,52 +379,51 @@ void householder_update(Matrix &A, int j, double beta, const Vector& vjm) {
 // update A, b
 // A' \define A_{S}-ar and b'\define b-ad
 // __attribute__ ((noinline))	// uncomment to prevent inlining when profiling
-static void updateAb(Matrix& A, Vector& b, int j, const Vector& a,
+inline void updateAB_manual(Matrix& A, Vector& b, int j, const Vector& a,
 		const Vector& r, double d) {
 	const size_t m = A.size1(), n = A.size2();
-#ifdef GSL
+	for (size_t i = 0; i < m; i++) { // update all rows
+		double ai = a(i);
+		b(i) -= ai * d;
+		double *Aij = A.data().begin() + i * n + j + 1;
+		const double *rptr = r.data().begin() + j + 1;
+		// A(i,j+1:end) -= ai*r(j+1:end)
+		for (size_t j2 = j + 1; j2 < n; j2++, Aij++, rptr++)
+			*Aij -= ai * (*rptr);
+	}
+}
 
+static void updateAb(Matrix& A, Vector& b, int j, const Vector& a,
+		const Vector& r, double d) {
+#ifdef GSL
+#ifndef REVERTGSL
+	const size_t m = A.size1(), n = A.size2();
 	// update A
 	// A(0:m,j+1:end) = A(0:m,j+1:end) - a(0:m)*r(j+1:end)'
 	// get a view for A
-//	gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
+	gsl_matrix_view Ag = gsl_matrix_view_array(A.data().begin(), m, n);
-//	gsl_matrix_view Ag_view = gsl_matrix_submatrix (&(Ag.matrix), 0, j+1, m, n-j-1);
+	gsl_matrix_view Ag_view = gsl_matrix_submatrix (&(Ag.matrix), 0, j+1, m, n-j-1);
-//	// get a view for r
+	// get a view for r
-//	gsl_vector_const_view rg = gsl_vector_const_view_array(r.data().begin()+j+1, n-j-1);
+	gsl_vector_const_view rg = gsl_vector_const_view_array(r.data().begin()+j+1, n-j-1);
-//	// get a view for a
+	// get a view for a
-//	gsl_vector_const_view ag = gsl_vector_const_view_array(a.data().begin(), m);
+	gsl_vector_const_view ag = gsl_vector_const_view_array(a.data().begin(), m);
-//
-//	// rank one update
-//	gsl_blas_dger (-1.0, &(ag.vector), &(rg.vector), &(Ag_view.matrix));
-//
-//	// update b
-//	double * bptr = b.data().begin();
-//	const double * aptr = a.data().begin();
-//	for (size_t i = 0; i < m; i++) {
-//		*(bptr+i) -= d* *(aptr+i);
-//	}
 
-	// original
+	// rank one update
-	for (size_t i = 0; i < m; i++) { // update all rows
+	gsl_blas_dger (-1.0, &(ag.vector), &(rg.vector), &(Ag_view.matrix));
-		double ai = a(i);
+
-		b(i) -= ai * d;
+	// update b
-		double *Aij = A.data().begin() + i * n + j + 1;
+	double * bptr = b.data().begin();
-		const double *rptr = r.data().begin() + j + 1;
+	const double * aptr = a.data().begin();
-		// A(i,j+1:end) -= ai*r(j+1:end)
+	for (size_t i = 0; i < m; i++) {
-		for (size_t j2 = j + 1; j2 < n; j2++, Aij++, rptr++)
+		*(bptr+i) -= d* *(aptr+i);
-			*Aij -= ai * (*rptr);
 	}
 
 #else
-	for (size_t i = 0; i < m; i++) { // update all rows
+	updateAB_manual(A, b,j,a,r,d);
-		double ai = a(i);
+#endif
-		b(i) -= ai * d;
+
-		double *Aij = A.data().begin() + i * n + j + 1;
+#else
-		const double *rptr = r.data().begin() + j + 1;
+	updateAb_manual(A,b,j,a,r,d);
-		// A(i,j+1:end) -= ai*r(j+1:end)
-		for (size_t j2 = j + 1; j2 < n; j2++, Aij++, rptr++)
-			*Aij -= ai * (*rptr);
-	}
 #endif
 }
 
@@ -458,11 +484,36 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
  * version with Householder vectors below diagonal, as in GVL
  */
 /* ************************************************************************* */
+inline void householder_manual(Matrix &A, size_t k) {
+	const size_t m = A.size1(), n = A.size2(), kprime = min(k,min(m,n));
+	// loop over the kprime first columns
+	for(size_t j=0; j < kprime; j++){
+		// below, the indices r,c always refer to original A
+
+		// copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
+		Vector xjm(m-j);
+		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);
+
+		// do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
+		householder_update(A, j, beta, vjm);
+
+		// the Householder vector is copied in the zeroed out part
+		for( size_t r = j+1 ; r < m ; r++ )
+			A(r,j) = vjm(r-j);
+
+	} // column j
+}
+
 void householder_(Matrix &A, size_t k) 
 {
-	const size_t m = A.size1(), n = A.size2(), kprime = min(k,min(m,n));
-
 #ifdef GSL
+#ifndef REVERTGSL
+	const size_t m = A.size1(), n = A.size2(), kprime = min(k,min(m,n));
 	// loop over the kprime first columns
 	for(size_t j=0; j < kprime; j++){
 		// below, the indices r,c always refer to original A
@@ -510,27 +561,11 @@ void householder_(Matrix &A, size_t k)
 	} // column j
 
 #else
-	// loop over the kprime first columns
+	householder_manual(A, k);
-	for(size_t j=0; j < kprime; j++){
+#endif
-		// below, the indices r,c always refer to original A
 
-		// copy column from matrix to xjm, i.e. x(j:m) = A(j:m,j)
+#else
-		Vector xjm(m-j);
+	householder_manual(A, 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);
-
-		// do outer product update A = (I-beta vv')*A = A - v*(beta*A'*v)' = A - v*w'
-		householder_update(A, j, beta, vjm);
-
-		// the Householder vector is copied in the zeroed out part
-		for( size_t r = j+1 ; r < m ; r++ )
-			A(r,j) = vjm(r-j);
-
-	} // column j
 #endif
 }
 

cpp/Matrix.h

@@ -133,6 +133,16 @@ void print(const Matrix& A, const std::string& s = "");
  */
 Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2);
 
+/**
+ * insert a submatrix IN PLACE at a specified location in a larger matrix
+ * NOTE: there is no size checking
+ * @param large matrix to be updated
+ * @param small matrix to be inserted
+ * @param i is the row of the upper left corner insert location
+ * @param j is the column of the upper left corner insert location
+ */
+void insertSub(Matrix& big, const Matrix& small, size_t i, size_t j);
+
 /**
  * extracts a column from a matrix
  * NOTE: using this without the underscore is the ublas version!
@@ -142,6 +152,17 @@ Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2);
  */
 Vector column_(const Matrix& A, size_t j);
 
+/**
+ * inserts a column into a matrix IN PLACE
+ * NOTE: there is no size checking
+ * Alternate form allows for vectors smaller than the whole column to be inserted
+ * @param A matrix to be modified in place
+ * @param col is the vector to be inserted
+ * @param j is the index to insert the column
+ */
+void insertColumn(Matrix& A, const Vector& col, size_t j);
+void insertColumn(Matrix& A, const Vector& col, size_t i, size_t j);
+
 /**
  * extracts a row from a matrix
  * @param matrix to extract row from
@@ -164,6 +185,12 @@ void solve(Matrix& A, Matrix& B);
  */
 Matrix inverse(const Matrix& A);
 
+/**
+ * Perform updates of system matrices
+ */
+static void updateAb(Matrix& A, Vector& b, int j, const Vector& a,
+		const Vector& r, double d);
+
 /**
  * QR factorization, inefficient, best use imperative householder below
  * m*n matrix -> m*m Q, m*n R

cpp/testMatrix.cpp

@@ -10,6 +10,7 @@
 #include <boost/tuple/tuple.hpp>
 #include <boost/foreach.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/io.hpp>
 #include "Matrix.h"
 
 using namespace std;
@@ -165,6 +166,45 @@ TEST( matrix, column )
 	CHECK(assert_equal(a3, exp3));
 }
 
+/* ************************************************************************* */
+TEST( matrix, insert_column )
+{
+	Matrix big = zeros(5,6);
+	Vector col = ones(5);
+	size_t j = 3;
+
+	insertColumn(big, col, j);
+
+	Matrix expected = Matrix_(5,6,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0);
+
+	CHECK(assert_equal(expected, big));
+}
+
+/* ************************************************************************* */
+TEST( matrix, insert_subcolumn )
+{
+	Matrix big = zeros(5,6);
+	Vector col = ones(2);
+	size_t i = 1;
+	size_t j = 3;
+
+	insertColumn(big, col, i, j);
+
+	Matrix expected = Matrix_(5,6,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+	CHECK(assert_equal(expected, big));
+}
+
 /* ************************************************************************* */
 TEST( matrix, row )
 {
@@ -198,6 +238,24 @@ TEST( matrix, zeros )
   EQUALITY(A , zero);
 }
 
+/* ************************************************************************* */
+TEST( matrix, insert_sub )
+{
+	Matrix big = zeros(5,6),
+		   small = Matrix_(2,3, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0);
+
+	insertSub(big, small, 1, 2);
+
+	Matrix expected = Matrix_(5,6,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+			0.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+			0.0, 0.0, 1.0, 1.0, 1.0, 0.0,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+			0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+	CHECK(assert_equal(expected, big));
+}
+
 /* ************************************************************************* */
 TEST( matrix, scale_columns )
 {

cpp/timeMatrix.cpp

@@ -198,6 +198,31 @@ double timeHouseholder(size_t reps) {
 	}
 	return elapsed;
 }
+/**
+ * Results: (Alex's machine)
+ * reps: 200000
+ *
+ * Initial (boost matrix proxies) - 12.08
+ * Direct pointer method          - 4.62
+ */
+double timeMatrixInsert(size_t reps) {
+	// create a matrix
+	Matrix bigBase = zeros(100, 100);
+	Matrix small = eye(5,5);
+
+	// perform timing
+	double elapsed;
+	{
+		boost::timer t;
+		Matrix big = bigBase;
+		for (size_t rep=0; rep<reps; ++rep)
+			for (size_t i=0; i<100; i += 5)
+				for (size_t j=0; j<100; j += 5)
+					insertSub(big, small, i,j);
+		elapsed = t.elapsed();
+	}
+	return elapsed;
+}
 
 int main(int argc, char ** argv) {
 
@@ -233,5 +258,11 @@ int main(int argc, char ** argv) {
 	double house_time = timeHouseholder(reps_house);
 	cout << "Elapsed time for householder_() : " << house_time << " sec" << endl;
 
+	// Time matrix insertion
+	cout << "Starting insertSub() Timing" << endl;
+	size_t reps_insert = 200000;
+	double insert_time = timeMatrixInsert(reps_insert);
+	cout << "Elapsed time for insertSub() : " << insert_time << " sec" << endl;
+
 	return 0;
 }