[Toon-members] TooN Cholesky.h gaussian_elimination.h

[Tom Drummond]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Tom Drummond <twd20>    09/03/10 16:05:27

Modified files:
        .              : Cholesky.h gaussian_elimination.h 

Log message:
        fixed gaussian_elimination.h to work with new TooN
        lots of work to do on Cholesky.h

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.18&r2=1.19
http://cvs.savannah.gnu.org/viewcvs/TooN/gaussian_elimination.h?cvsroot=toon&r1=1.2&r2=1.3

Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.18
retrieving revision 1.19
diff -u -b -r1.18 -r1.19
--- Cholesky.h  1 Aug 2008 20:37:01 -0000       1.18
+++ Cholesky.h  10 Mar 2009 16:05:26 -0000      1.19
@@ -433,6 +433,8 @@
        int rank;
     };
   
+
+
     template <>
     class Cholesky<-1> {
     public:

Index: gaussian_elimination.h
===================================================================
RCS file: /cvsroot/toon/TooN/gaussian_elimination.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -b -r1.2 -r1.3
--- gaussian_elimination.h      29 May 2008 16:33:08 -0000      1.2
+++ gaussian_elimination.h      10 Mar 2009 16:05:26 -0000      1.3
@@ -5,62 +5,52 @@
 #include <TooN/TooN.h>
 
 namespace TooN {
-    template <class R, int N, class B> inline
-    R gaussian_elimination(Matrix<N> A, B b)
-    {
+    template<int N, typename Precision>
+       inline Vector<N, Precision>::type>  gaussian_elimination 
(Matrix<N,N,Precision> A, Vector<N, Precision> b) {
        using std::swap;
 
-       for (int i=0; i<N; ++i) {
+               int size=b.size();
+
+               for (int i=0; i<size; ++i) {
            int argmax = i;
-           double maxval = abs(A[i][i]);
+                       Precision maxval = abs(A[i][i]);
        
-           for (int ii=i+1; ii<N; ++ii) {
+                       for (int ii=i+1; ii<size; ++ii) {
                double v =  abs(A[ii][i]);
                if (v > maxval) {
                    maxval = v;
                    argmax = ii;
                }
            }
-           double pivot = A[argmax][i];
+                       Precision pivot = A[argmax][i];
            //assert(abs(pivot) > 1e-16);
-           double inv_pivot = 1.0/pivot;
+                       Precision inv_pivot = static_cast<Precision>(1)/pivot;
            if (argmax != i) {
-               for (int j=i; j<N; ++j)
+                               for (int j=i; j<size; ++j)
                    swap(A[i][j], A[argmax][j]);
                swap(b[i], b[argmax]);
            }
            //A[i][i] = 1;
-           for (int j=i+1; j<N; ++j)
+                       for (int j=i+1; j<size; ++j)
                A[i][j] *= inv_pivot;
            b[i] *= inv_pivot;
 
-           for (int u=i+1; u<N; ++u) {
+                       for (int u=i+1; u<size; ++u) {
                double factor = A[u][i];
                //A[u][i] = 0;
-               for (int j=i+1; j<N; ++j)
+                               for (int j=i+1; j<size; ++j)
                    A[u][j] -= factor * A[i][j];
                b[u] -= factor * b[i];
            }
        }
     
-       R x;
-       for (int i=N-1; i>=0; --i) {
+               Vector<N,Precision>::type> x(size);
+               for (int i=size-1; i>=0; --i) {
            x[i] = b[i];
-           for (int j=i+1; j<N; ++j)
+                       for (int j=i+1; j<size; ++j)
                x[i] -= A[i][j] * x[j];
        }
        return x;
     }
-
-    template <int N, class A1, class A2>
-    inline Vector<N> gaussian_elimination(const FixedMatrix<N,N,A1>& A, const 
FixedVector<N,A2>& b) {
-       return gaussian_elimination<Vector<N>, N, Vector<N> >(Matrix<N>(A), 
Vector<N>(b));
-    }
-
-    template <int N, class A1, class A2, int M>
-    inline Matrix<N> gaussian_elimination(const FixedMatrix<N,N,A1>& A, const 
FixedMatrix<N,M,A2>& b) {
-       return gaussian_elimination<Matrix<N,M>,N, Matrix<N,M> >(Matrix<N>(A), 
Matrix<N,M>(b));
-    }
-
 }
 #endif