[Toon-members] TooN TooN.h determinant.h

[Edward Rosten]

CVSROOT:        /cvsroot/toon
Module name:    TooN
Changes by:     Edward Rosten <edrosten>        09/06/24 13:55:40

Modified files:
        .              : TooN.h 
Added files:
        .              : determinant.h 

Log message:
        Added determinant function. Switches between obvious 2x2 method, 
        modified version of gaussian_elimination.h and LU.
        
        The crossover point to LU is determined by the TOON_DETERMINANT_LAPACK 
macro
        if this is not defined then it will never use LU. At some point this 
will
        be moved in to the configuration system, but at the moment it's 
#define'd
        in TooN.h

CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/TooN.h?cvsroot=toon&r1=1.45&r2=1.46
http://cvs.savannah.gnu.org/viewcvs/TooN/determinant.h?cvsroot=toon&rev=1.1

Patches:
Index: TooN.h
===================================================================
RCS file: /cvsroot/toon/TooN/TooN.h,v
retrieving revision 1.45
retrieving revision 1.46
diff -u -b -r1.45 -r1.46
--- TooN.h      18 Jun 2009 14:52:49 -0000      1.45
+++ TooN.h      24 Jun 2009 13:55:40 -0000      1.46
@@ -44,6 +44,8 @@
 #include <ctime>
 #endif
 
+#define TOON_DETERMINANT_LAPACK 30
+
 namespace TooN {
 
 #ifdef TOON_TEST_INTERNALS

Index: determinant.h
===================================================================
RCS file: determinant.h
diff -N determinant.h
--- /dev/null   1 Jan 1970 00:00:00 -0000
+++ determinant.h       24 Jun 2009 13:55:40 -0000      1.1
@@ -0,0 +1,121 @@
+#ifndef TOON_INCLUDE_DETERMINANT_H
+#define TOON_INCLUDE_DETERMINANT_H
+#ifdef TOON_DETERMINANT_LAPACK
+       #include <TooN/LU.h>
+#endif
+
+namespace TooN
+{
+       namespace Internal
+       {
+               template<int R, int C> struct Square
+               {
+               };
+
+               template<int R> struct Square<R, R>
+               {
+                       static const int Size = R;
+               };
+
+               template<int R> struct Square<R, Dynamic>
+               {
+                       static const int Size = R;
+               };
+               template<int C> struct Square<Dynamic, C>
+               {
+                       static const int Size = C;
+               };
+               template<> struct Square<Dynamic, Dynamic>
+               {
+                       static const int Size = Dynamic;
+               };
+       };
+
+
+       template<int R, int C, typename Precision, typename Base>
+       Precision determinant_gaussian_elimination(const Matrix<R, C, 
Precision, Base>& A_)
+       {
+               Matrix<Internal::Square<R,C>::Size, 
Internal::Square<R,C>::Size,Precision> A = A_;
+               TooN::SizeMismatch<R, C>::test(A.num_rows(), A.num_cols());
+               using std::swap;
+
+               int size=A.num_rows();
+               
+               //If row operations of the form row_a += alpha * row_b
+               //then the determinant is unaffected. However, if a row
+               //is scaled, then the determinant is scaled by the same 
+               //amount. The total scaling is held in divisor.
+               Precision determinant=1;
+               Precision divisor=1;
+
+               for (int i=0; i<size; ++i) {
+
+                       //Find the pivot
+                       int argmax = i;
+                       Precision maxval = abs(A[i][i]);
+                       
+                       for (int ii=i+1; ii<size; ++ii) {
+                               double v =  abs(A[ii][i]);
+                               if (v > maxval) {
+                                       maxval = v;
+                                       argmax = ii;
+                               }
+                       }
+                       Precision pivot = A[argmax][i];
+
+                       //assert(abs(pivot) > 1e-16);
+                       
+                       //Swap the current row with the pivot row if necessary.
+                       if (argmax != i) {
+                               for (int j=i; j<size; ++j)
+                                       swap(A[i][j], A[argmax][j]);
+                       }
+
+                       determinant *= A[i][i];
+
+                       if(determinant == 0)
+                               return 0;
+                       
+                       for (int u=i+1; u<size; ++u) {
+                               //Multiply out the usual 1/pivot term
+                               //to avoid division.
+                               double factor = A[u][i];
+                               //A[u][i] = 0;
+                               divisor*=pivot;
+                               for (int j=i+1; j<size; ++j)
+                                       A[u][j] = A[u][j]*pivot - factor * 
A[i][j];
+                       }
+               }
+
+               if(size %2 == 0) 
+                       return -determinant/divisor;
+               else
+                       return determinant/divisor;
+       }
+       
+       #ifdef TOON_DETERMINANT_LAPACK
+               template<int R, int C, class P, class B>
+               P determinant_LU(const Matrix<R, C, P, B>& A)
+               {
+                       TooN::SizeMismatch<R, C>::test(A.num_rows(), 
A.num_cols());
+                       LU<Internal::Square<R,C>::Size, P> lu(A);
+                       return lu.determinant();
+               }
+       #endif
+
+       template<int R, int C, class P, class B>
+       P determinant(const Matrix<R, C, P, B>& A)
+       {
+               TooN::SizeMismatch<R, C>::test(A.num_rows(), A.num_cols());
+               if(A.num_rows() == 2)
+                       return A[0][0]*A[1][1] - A[1][0]*A[0][1];
+               #ifdef TOON_DETERMINANT_LAPACK
+                       else if(A.num_rows() > TOON_DETERMINANT_LAPACK)
+                               return determinant_LU(A);
+               #endif
+               else
+                       return determinant_gaussian_elimination(A);
+       }
+}
+
+#endif