[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN gauss_jordan.h helpers.h benchmark/solve_a...
From: |
Edward Rosten |
Subject: |
[Toon-members] TooN gauss_jordan.h helpers.h benchmark/solve_a... |
Date: |
Fri, 20 Mar 2009 18:24:29 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 09/03/20 18:24:29
Modified files:
. : gauss_jordan.h helpers.h
benchmark : solve_ax_equals_b.cc
internal : matrix.hh
Added files:
test : gauss_jordan.cc
Log message:
Added 0-ary operators for Matrix.
You can now do:
mat.slice<0,0,2,2>() = Identity;
The old version:
Identity(mat.slice<0,0,2,2>());
failed due to slice objects being temporaries.
Added Gauss-Jordan test which relies on this, and fixed
gauss-jordan to work with TooN 2.
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/gauss_jordan.h?cvsroot=toon&r1=1.2&r2=1.3
http://cvs.savannah.gnu.org/viewcvs/TooN/helpers.h?cvsroot=toon&r1=1.27&r2=1.28
http://cvs.savannah.gnu.org/viewcvs/TooN/benchmark/solve_ax_equals_b.cc?cvsroot=toon&r1=1.2&r2=1.3
http://cvs.savannah.gnu.org/viewcvs/TooN/internal/matrix.hh?cvsroot=toon&r1=1.16&r2=1.17
http://cvs.savannah.gnu.org/viewcvs/TooN/test/gauss_jordan.cc?cvsroot=toon&rev=1.1
Patches:
Index: gauss_jordan.h
===================================================================
RCS file: /cvsroot/toon/TooN/gauss_jordan.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -b -r1.2 -r1.3
--- gauss_jordan.h 20 Mar 2009 17:35:31 -0000 1.2
+++ gauss_jordan.h 20 Mar 2009 18:24:27 -0000 1.3
@@ -15,12 +15,12 @@
/// partial pivoting.
///
/// @param m The matrix to be reduced.
-template<int R, int C, class X> void gauss_jordan(FixedMatrix<R, C, X>& m)
+template<int R, int C, class Precision, class Base> void
gauss_jordan(Matrix<R, C, Precision, Base>& m)
{
using std::swap;
//Loop over columns to reduce.
- for(int col=0; col < R; col++)
+ for(int col=0; col < m.num_rows(); col++)
{
//Reduce the current column to a single element
@@ -34,19 +34,21 @@
{
int pivotpos = col;
double pivotval = abs(m[pivotpos][col]);
- for(int p=col+1; p <R; p++)
+ for(int p=col+1; p <m.num_rows(); p++)
if(abs(m[p][col]) > pivotval)
{
pivotpos = p;
pivotval = abs(m[pivotpos][col]);
}
- swap(m[col], m[pivotpos]);
+ if(col != pivotpos)
+ for(int c=0; c < m.num_cols(); c++)
+ swap(m[col][c], m[pivotpos][c]);
}
//Reduce the current column in every row to zero, excluding
elements on
//the leading diagonal.
- for(int row = 0; row < R; row++)
+ for(int row = 0; row < m.num_rows(); row++)
{
if(row != col)
{
@@ -59,7 +61,7 @@
//Subtract the pivot row from all other rows,
to make
//column col zero.
m[row][col] = 0;
- for(int c=col+1; c < C; c++)
+ for(int c=col+1; c < m.num_cols(); c++)
m[row][c] = m[row][c] - m[col][c] *
multiple;
}
}
@@ -68,13 +70,13 @@
//Final pass to make diagonal elements one. Performing this in a final
//pass allows us to avoid any significant computations on the left-hand
//square matrix, since it is diagonal, and ends up as the identity.
- for(int row=0;row < R; row++)
+ for(int row=0;row < m.num_rows(); row++)
{
double mul = 1/m[row][row];
m[row][row] = 1;
- for(int col=R; col < C; col++)
+ for(int col=m.num_rows(); col < m.num_cols(); col++)
m[row][col] *= mul;
}
}
Index: helpers.h
===================================================================
RCS file: /cvsroot/toon/TooN/helpers.h,v
retrieving revision 1.27
retrieving revision 1.28
diff -u -b -r1.27 -r1.28
--- helpers.h 20 Mar 2009 16:49:25 -0000 1.27
+++ helpers.h 20 Mar 2009 18:24:27 -0000 1.28
@@ -40,16 +40,6 @@
}
-template<int Rows, int Cols, class Precision, class Base> void
Identity(Matrix<Rows, Cols, Precision, Base>& m)
-{
- SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());
-
- Zero(m);
- for(int i=0; i < m.num_rows(); i++)
- m[i][i] = 1;
-}
-
-
template<int Size, class Precision, class Base> void Fill(Vector<Size,
Precision, Base>& v, const Precision& p)
{
for(int i=0; i < v.size(); i++)
@@ -65,5 +55,24 @@
+namespace Internal{
+
+struct Identity
+{
+ template<int R, int C, class P, class B> static void eval(Matrix<R, C,
P, B>& m)
+ {
+ SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
+
+ for(int r=0; r < m.num_rows(); r++)
+ for(int c=0; c < m.num_rows(); c++)
+ m[r][c] = 0;
+
+ for(int i=0; i < m.num_rows(); i++)
+ m[i][i] = 1;
+ }
+};
+}
+static Operator<Internal::Identity> Identity;
+
}
#endif
Index: benchmark/solve_ax_equals_b.cc
===================================================================
RCS file: /cvsroot/toon/TooN/benchmark/solve_ax_equals_b.cc,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -b -r1.2 -r1.3
--- benchmark/solve_ax_equals_b.cc 20 Mar 2009 17:21:30 -0000 1.2
+++ benchmark/solve_ax_equals_b.cc 20 Mar 2009 18:24:27 -0000 1.3
@@ -82,11 +82,11 @@
template<int Size, int Cols, class Solver> void benchmark_ax_eq_b()
{
- double time=0, t_tmp;
+ double time=0, t_tmp, start = get_time_of_day();
double sum=0;
int n=0;
- while(time < 1)
+ while(get_time_of_day() - start < .1)
{
Matrix<Size> a;
for(int r=0; r < Size; r++)
@@ -112,7 +112,7 @@
n++;
}
- cout << Solver::name() << " " << Size << " " << Cols << " " << n / time
<< endl;
+ cout << Solver::name() << "\t" << n / time << "\t";
global_sum += sum;
}
@@ -149,8 +149,10 @@
{
static void iter()
{
+ cout << Size << "\t" << Cols << "\t";
benchmark_iter<Size, Cols, Test>::iter();
- ColIter<Size, Cols-50, Test>::iter();
+ cout << endl;
+ ColIter<Size, Cols-5, Test>::iter();
}
};
Index: internal/matrix.hh
===================================================================
RCS file: /cvsroot/toon/TooN/internal/matrix.hh,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -b -r1.16 -r1.17
--- internal/matrix.hh 20 Mar 2009 16:49:25 -0000 1.16
+++ internal/matrix.hh 20 Mar 2009 18:24:27 -0000 1.17
@@ -24,7 +24,6 @@
:Layout::template Layout<Rows,Cols,Precision>(rows, cols)
{}
-
//See vector.hh and allocator.hh for details about why the
//copy constructor should be default.
@@ -66,6 +65,12 @@
return *this;
}
+ // operator = 0-ary operator
+ template<class Op> inline Matrix& operator= (const Operator<Op>&)
+ {
+ Op::eval(*this);
+ }
+
// operator =
template<int Rows2, int Cols2, typename Precision2, typename Base2>
Matrix& operator= (const Matrix<Rows2, Cols2, Precision2, Base2>& from)
Index: test/gauss_jordan.cc
===================================================================
RCS file: test/gauss_jordan.cc
diff -N test/gauss_jordan.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ test/gauss_jordan.cc 20 Mar 2009 18:24:29 -0000 1.1
@@ -0,0 +1,34 @@
+#include <TooN/gauss_jordan.h>
+#include <TooN/helpers.h>
+#include <tr1/random>
+#include <fstream>
+
+using namespace TooN;
+using namespace std;
+using namespace tr1;
+
+int main()
+{
+ unsigned int s;
+ ifstream("/dev/urandom").read((char*)&s, 4);
+
+ std::tr1::mt19937 eng;
+ std::tr1::uniform_real<double> rnd;
+ eng.seed(s);
+ Matrix<5,10> m, reduced;
+
+ for(int i=0; i< m.num_rows(); i++)
+ for(int j=0; j< m.num_rows(); j++)
+ m[i][j] = rnd(eng);
+
+ m.slice<0, 5, 5, 5>() = Identity;
+
+ cout << m << endl;
+
+
+ reduced = m;
+ gauss_jordan(reduced);
+
+ cout << reduced.slice<0,5,5,5>() * m.slice<0,0,5,5>() << endl;
+
+}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN gauss_jordan.h helpers.h benchmark/solve_a...,
Edward Rosten <=