Merge branch 'master' into expression_template

.gitignore

@@ -0,0 +1,19 @@
+regressions/*.result
+regressions/*.test
+regressions/*.testout
+tmp
+a.out
+*.s
+html/
+autom4te.cache/
+*/a.out
+config.log
+config.status
+*.gcda
+Makefile
+.*.swp
+.*.swo
+.*.swn
+regressions/results
+*.o
+TooN-doxy.tag

Makefile.in

@@ -44,7 +44,7 @@ docs:
 	doxygen 
 
 
-TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant chol_toon chol_lapack simplex sym_eigen fill so3 complex qr gr_svd diagonal_matrix
+TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant chol_toon chol_lapack simplex sym_eigen fill so3 complex qr gr_svd diagonal_matrix gaussian_elimination
 
 
 TEST_RESULT=$(TESTS:%=regressions/%.result)

QR.h

@@ -35,7 +35,20 @@
 
 namespace TooN
 {
+/**
+Performs %QR decomposition.
 
+@warning this will only work if the number of columns is greater than 
+the number of rows!
+
+The QR decomposition operates on a matrix A. 
+In general:
+\f[
+A = QR
+\f]
+
+@ingroup gDecomps
+*/
 template<int Rows=Dynamic, int Cols=Rows, typename Precision=double> class QR
 {
 

doc/documentation.h

@@ -598,10 +598,12 @@ This section is arranged as a FAQ. Most answers include code fragments. Assume
 	
 	Similarly for the other \ref sDecompos "decomposition objects"
 
+	For 2x2 matrices, the TooN::inv function can ve used.
+
 	\subsection sDecompos  Which decomposisions are there?
 
 	For general size matrices (not necessarily square) there are:
-	@link TooN::LU LU @endlink, @link TooN::SVD SVD @endlink and gauss_jordan()
+	@link TooN::LU LU @endlink, @link TooN::SVD SVD @endlink, @link TooN::QR QR@endlink, @link TooN::QR_Lapack LAPACK's QR@endlink and gauss_jordan()
 
 	For square symmetric matrices there are:
 	@link TooN::SymEigen SymEigen @endlink and @link TooN::Cholesky Cholesky @endlink
@@ -950,7 +952,8 @@ around LAPACK.
 
 /// @defgroup gDecomps Matrix decompositions
 /// Classes to perform matrix decompositions, used to solve 
-/// linear equations and provide information about matrices. These are wrappers for functionality
+/// linear equations and provide information about matrices. 
+/// Some of these are wrappers around LAPACK, others are built in.
 /// provided by the LAPACK library.
 
 /// @defgroup gTransforms Transformation matrices

gaussian_elimination.h

@@ -94,7 +94,7 @@ namespace TooN {
 	{
 		template<int i, int j, int k> struct Size3
 		{
-			static const int s=(i!= -1)?i:(j!=-1?j:k);
+			static const int s = !IsStatic<i>::is?i: (!IsStatic<j>::is?j:k);
 		};
 
 	};

helpers.h

@@ -56,6 +56,24 @@ namespace std {
 #endif
 
 namespace TooN {
+
+	///Invert a matrix.
+	///
+	///For sizes other than 2x2, @link gDecomps decompositions @endlink provide a suitable solition.
+	///
+	///@param m Matrix to invert.
+	///@ingroup gDecomps
+	inline Matrix<2> inv(const Matrix<2>& m)
+	{
+		double d = 1./(m[0][0]*m[1][1] - m[1][0]*m[0][1]);
+
+		return Data(
+		     m[1][1]*d, -m[0][1]*d,
+		    -m[1][0]*d,  m[0][0]*d
+		);
+	}
+
+
 
 	///\deprecated
 	///@ingroup gLinAlg

internal/matrix.hh

@@ -63,9 +63,7 @@ Matrix<3> M;
 which just is a synonym for <code>Matrix<3,3></code>. Matrices can also be
 constructed from pointers or static 1D or 2D arrays of doubles:
 @code
-
-  double dvals1[9]={1,2,3,4,5,6};
-  Matrix<2,3, Reference::RowMajor> M2 (dvals1);
+  Matrix<2,3, Reference::RowMajor> M2 = Data(1,2,3,4,5,6);
 @endcode
 
 \par Dynamically-sized matrices
@@ -85,8 +83,7 @@ note that the static dimension must be provided, but it is ignored.
 
 @endcode
 
-<code>Matrix<></code> is a synonym for <code> Matrix<Dynamic, Dynamic> </code> which is
-<code>%Matrix<-1,-1></code>
+<code>Matrix<></code> is a synonym for <code> Matrix<Dynamic, Dynamic> </code>.
 
 \par Row-major and column-major
 
@@ -100,13 +97,13 @@ You can override the default for a specific matrix by specifying the layout when
 you construct it:
 @code
 Matrix<3,3,double,ColMajor> M1;
-Matrix<-1,-1,double,RowMajor> M2(nrows, ncols);
+Matrix<Dynamic,Dynamic,double,RowMajor> M2(nrows, ncols);
 @endcode
 In this case the precision template argument must be given as it precedes the layout argument
 
 @ingroup gLinAlg
 **/
-template <int Rows=-1, int Cols=Rows, class Precision=DefaultPrecision, class Layout = RowMajor>
+template <int Rows=Dynamic, int Cols=Rows, class Precision=DefaultPrecision, class Layout = RowMajor>
 struct Matrix : public Layout::template MLayout<Rows, Cols, Precision>
 {
 public:
@@ -305,8 +302,9 @@ public:
 		Access an element from the matrix.
 		The index starts at zero, i.e. the top-left element is m(0, 0).
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
-		Matrix<2,3> m(d);
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		double e = m(1,2);     // now e = 6.0;
 		@endcode
 		@internal
@@ -331,7 +329,9 @@ public:
 		This can be used as either an r-value or an l-value. The index starts at zero,
 		i.e. the top-left element is m(0, 0).
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		m(1,2) = 8;     // now d = [1 2 3]
 					  //         [4 5 8]
@@ -349,7 +349,9 @@ public:
 		or an l-value. The index starts at zero, i.e. the first row (or column) is
 		m[0].
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		Vector<3> v = m[1];       // now v = [4 5 6];
 		Vector<2> v2 = m.T()[0];  // now v2 = [1 4];
@@ -367,7 +369,9 @@ public:
 		or an l-value. The index starts at zero, i.e. the first row (or column) is
 		m[0].
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		Zero(m[0]);   // set the first row to zero
 		Vector<2> v = 8,9;
@@ -396,7 +400,9 @@ public:
 		reinterprets a row-major matrix as column-major or vice-versa. This can be
 		used as an l-value.
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		Zero(m[0]);   // set the first row to zero
 		Vector<2> v = 8,9;
@@ -413,7 +419,9 @@ public:
 		reinterprets a row-major  matrix as column-major or vice-versa. The result can
 		be used as an l-value.
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		Vector<2> v = 8,9;
 		// Set the first column to v
@@ -433,7 +441,10 @@ public:
 		Extract a sub-matrix. The matrix extracted will be begin at element
 		(Rstart, Cstart) and will contain the next Rsize by Csize elements.
 		@code
-		double d[3][3] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6
+			7,8,9));
 		Matrix<3> m(d);
 		Extract the top-left 2x2 matrix
 		Matrix<2> b = m.slice<0,0,2,2>();  // b = [1 2]
@@ -450,7 +461,9 @@ public:
 		Cstart) and will contain the next Rsize by Csize elements. This can be used as
 		either an r-value or an l-value.
 		@code
-		double d[2][3] = {{1, 2, 3}, {4, 5, 6}};
+		Matrix<2,3> m(Data(
+			1,2,3
+			4,5,6));
 		Matrix<2,3> m(d);
 		Zero(m.slice<0,2,2,1>());  // b = [1 2 0]
 								  //     [4 5 0]

regressions/complex.cc

@@ -1,5 +1,5 @@
 #include <complex>
-#include "TooN/regressions/regression.h"
+#include "regressions/regression.h"
 #include <TooN/internal/planar_complex.hh>
 
 int main()

regressions/gaussian_elimination.cc

@@ -0,0 +1,40 @@
+#include "regressions/regression.h"
+#include <TooN/gaussian_elimination.h>
+using namespace TooN;
+using namespace std;
+
+
+int main()
+{
+	double err=0;
+
+	//Test a bunch of matrices 
+
+	for(int i=0; i < 10000; i++)
+	{
+		int rows = xor128u() % 100 + 2;
+		int cols = rows  + xor128u() % 20;
+
+		Matrix<> b(rows, cols);
+
+		for(int r=0; r < rows; r++)
+			for(int c=0; c < cols; c++)
+				b[r][c] = xor128d();
+
+		Matrix<> A(rows, rows);
+
+		for(int r=0; r < rows; r++)
+			for(int c=0; c < rows; c++)
+				A[r][c] = xor128d();
+
+		Matrix<> x = gaussian_elimination(A, b);
+
+		Matrix<> e = A*x - b;
+
+		for(int r=0; r < rows; r++)
+			err = max(err, norm_inf(e[r]));
+	}
+
+	cout << err << endl;
+}
+

regressions/gaussian_elimination.txt

@@ -0,0 +1 @@
+0