[Toon-members] TooN doc/Doxyfile gauss_jordan.h [Maintenance_Branch_1_x]

[Edward Rosten]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Branch:         Maintenance_Branch_1_x
Changes by:     Edward Rosten <edrosten>        08/12/17 18:59:48

Modified files:
        doc            : Doxyfile 
Added files:
        .              : gauss_jordan.h 

Log message:
        Gauss-Jordan reduction. The input matrix [A|B] is modified in place to 
        become [I|X], such that AX=B.
        
        On my machine:
        
        Size | Speed relative to LU::get_inverse()
        -----+------------------------------------
        3    |    8.6
        4    |    6.5
        10   |    2.1
        20   |    0.94

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/gauss_jordan.h?cvsroot=toon&only_with_tag=Maintenance_Branch_1_x&rev=1.1.2.1
http://cvs.savannah.gnu.org/viewcvs/TooN/doc/Doxyfile?cvsroot=toon&only_with_tag=Maintenance_Branch_1_x&r1=1.3&r2=1.3.2.1

Patches:
Index: doc/Doxyfile
===================================================================
RCS file: /cvsroot/toon/TooN/doc/Doxyfile,v
retrieving revision 1.3
retrieving revision 1.3.2.1
diff -u -b -r1.3 -r1.3.2.1
--- doc/Doxyfile        23 Oct 2007 21:35:46 -0000      1.3
+++ doc/Doxyfile        17 Dec 2008 18:59:48 -0000      1.3.2.1
@@ -400,7 +400,7 @@
 # *.c *.cc *.cxx *.cpp *.c++ *.java *.ii *.ixx *.ipp *.i++ *.inl *.h *.hh 
*.hxx *.hpp 
 # *.h++ *.idl *.odl *.cs *.php *.php3 *.inc
 
-FILE_PATTERNS          = documentation.h *doc.h irls.h wls.h wls_cholesky.h 
downhill_simplex.h
+FILE_PATTERNS          = documentation.h *doc.h irls.h wls.h wls_cholesky.h 
downhill_simplex.h gauss_jordan.h
 
 # The RECURSIVE tag can be used to turn specify whether or not subdirectories 
 # should be searched for input files as well. Possible values are YES and NO. 

Index: gauss_jordan.h
===================================================================
RCS file: gauss_jordan.h
diff -N gauss_jordan.h
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ gauss_jordan.h      17 Dec 2008 18:59:48 -0000      1.1.2.1
@@ -0,0 +1,83 @@
+#ifndef TOON_INC_GAUSS_JORDAN_H
+#define TOON_INC_GAUSS_JORDAN_H
+
+#include <utility>
+#include <TooN/TooN.h>
+
+namespace TooN
+{
+/// Perform Gauss-Jordan reduction on m
+///
+/// If m is of the form \f$[A | I ]\f$, then after reduction, m
+/// will be \f$[ I | A^{-1}]\f$. There is no restriction on the input, 
+/// in that the matrix augmenting A does not need to be I, or square.
+/// The reduction is performed using elementary row operations and 
+/// partial pivoting.
+///
+/// @param m The matrix to be reduced.
+template<int R, int C, class X> void gauss_jordan(FixedMatrix<R, C, X>& m)
+{
+       using std::swap;
+
+       //Loop over columns to reduce.
+       for(int col=0; col < R; col++)
+       {
+               //Reduce the current column to a single element
+
+
+               //Search down the current column in the lower triangle for the 
largest
+               //absolute element (pivot).  Then swap the pivot row, so that 
the pivot
+               //element is on the diagonal. The benchmarks show that it is 
actually
+               //faster to swap whole rows than it is to access the rows via 
indirection 
+               //and swap the indirection element. This holds for both pointer 
indirection
+               //and using a permutation vector over rows.
+               {
+                       int pivotpos = col;
+                       double pivotval = abs(m[pivotpos][col]);
+                       for(int p=col+1; p <R; p++)
+                               if(abs(m[p][col]) > pivotval)
+                               {
+                                       pivotpos = p;
+                                       pivotval = abs(m[pivotpos][col]);
+                               }
+
+                       swap(m[col], m[pivotpos]);
+               }
+
+               //Reduce the current column in every row to zero, excluding 
elements on
+               //the leading diagonal.
+               for(int row = 0; row < R; row++)
+               {
+                       if(row != col)
+                       {
+                               double multiple = m[row][col] / m[col][col];
+               
+                               //We could eliminate some of the computations 
in the augmented
+                               //matrix, if the augmented half is the 
identity. In general, it
+                               //is not. 
+
+                               //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++)
+                                       m[row][c] = m[row][c] - m[col][c] * 
multiple;
+                       }
+               }
+       }
+       
+       //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++)
+       {
+               double mul = 1/m[row][row];
+
+               m[row][row] = 1;
+
+               for(int col=R; col < C; col++)
+                       m[row][col] *= mul;
+       }
+}
+
+}
+#endif