[Toon-members] TooN gauss_jordan.h helpers.h benchmark/solve_a...

[Edward Rosten]

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;
+
+}