[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN TooN.h determinant.h
From: |
Edward Rosten |
Subject: |
[Toon-members] TooN TooN.h determinant.h |
Date: |
Wed, 24 Jun 2009 13:55:40 +0000 |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN TooN.h determinant.h,
Edward Rosten <=