[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN gaussian_elimination.h helpers.h lapack.h ...
|
From: |
Edward Rosten |
|
Subject: |
[Toon-members] TooN gaussian_elimination.h helpers.h lapack.h ... |
|
Date: |
Fri, 20 Mar 2009 16:49:26 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 09/03/20 16:49:26
Modified files:
. : gaussian_elimination.h helpers.h lapack.h
internal : matrix.hh vector.hh
Added files:
benchmark : solve_ax_equals_b.cc
test : gaussian_elimination_test.cc
Log message:
- Make gaussian_elimination work on Matrices
- Added test program for gaussian_elimination on matrices
- Added benchmark code ao test the speed of solve x for Ax=b
- Fixed all bugs which cropped up as a result
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/gaussian_elimination.h?cvsroot=toon&r1=1.3&r2=1.4
http://cvs.savannah.gnu.org/viewcvs/TooN/helpers.h?cvsroot=toon&r1=1.26&r2=1.27
http://cvs.savannah.gnu.org/viewcvs/TooN/lapack.h?cvsroot=toon&r1=1.8&r2=1.9
http://cvs.savannah.gnu.org/viewcvs/TooN/benchmark/solve_ax_equals_b.cc?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/internal/matrix.hh?cvsroot=toon&r1=1.15&r2=1.16
http://cvs.savannah.gnu.org/viewcvs/TooN/internal/vector.hh?cvsroot=toon&r1=1.23&r2=1.24
http://cvs.savannah.gnu.org/viewcvs/TooN/test/gaussian_elimination_test.cc?cvsroot=toon&rev=1.1
Patches:
Index: gaussian_elimination.h
===================================================================
RCS file: /cvsroot/toon/TooN/gaussian_elimination.h,v
retrieving revision 1.3
retrieving revision 1.4
diff -u -b -r1.3 -r1.4
--- gaussian_elimination.h 10 Mar 2009 16:05:26 -0000 1.3
+++ gaussian_elimination.h 20 Mar 2009 16:49:25 -0000 1.4
@@ -6,7 +6,7 @@
namespace TooN {
template<int N, typename Precision>
- inline Vector<N, Precision>::type> gaussian_elimination
(Matrix<N,N,Precision> A, Vector<N, Precision> b) {
+ inline Vector<N, Precision> gaussian_elimination (Matrix<N,N,Precision>
A, Vector<N, Precision> b) {
using std::swap;
int size=b.size();
@@ -44,7 +44,7 @@
}
}
- Vector<N,Precision>::type> x(size);
+ Vector<N,Precision> x(size);
for (int i=size-1; i>=0; --i) {
x[i] = b[i];
for (int j=i+1; j<size; ++j)
@@ -52,5 +52,69 @@
}
return x;
}
+
+ namespace Internal
+ {
+ template<int i, int j, int k> struct Size3
+ {
+ static const int s=(i!= -1)?i:(j!=-1?j:k);
+ };
+
+ };
+
+ template<int R1, int C1, int R2, int C2, typename Precision>
+ inline Matrix<Internal::Size3<R1, C1, R2>::s, C2, Precision>
gaussian_elimination (Matrix<R1,C1,Precision> A, Matrix<R2, C2, Precision> b) {
+ using std::swap;
+ SizeMismatch<R1, C1>::test(A.num_rows(), A.num_cols());
+ SizeMismatch<R1, R2>::test(A.num_rows(), b.num_rows());
+
+ int size=A.num_rows();
+
+ for (int i=0; i<size; ++i) {
+ 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);
+ Precision inv_pivot = static_cast<Precision>(1)/pivot;
+ if (argmax != i) {
+ for (int j=i; j<size; ++j)
+ swap(A[i][j], A[argmax][j]);
+
+ for(int j=0; j < b.num_cols(); j++)
+ swap(b[i][j], b[argmax][j]);
+ }
+ //A[i][i] = 1;
+ for (int j=i+1; j<size; ++j)
+ A[i][j] *= inv_pivot;
+ b[i] *= inv_pivot;
+
+ for (int u=i+1; u<size; ++u) {
+ double factor = A[u][i];
+ //A[u][i] = 0;
+ for (int j=i+1; j<size; ++j)
+ A[u][j] -= factor * A[i][j];
+ b[u] -= factor * b[i];
+ }
+ }
+
+ Matrix<Internal::Size3<R1, C1, R2>::s,C2,Precision>
x(b.num_rows(), b.num_cols());
+ for (int i=size-1; i>=0; --i) {
+ for(int k=0; k <b.num_cols(); k++)
+ {
+ x[i][k] = b[i][k];
+ for (int j=i+1; j<size; ++j)
+ x[i][k] -= A[i][j] * x[j][k];
+ }
+ }
+ return x;
+ }
}
#endif
Index: helpers.h
===================================================================
RCS file: /cvsroot/toon/TooN/helpers.h,v
retrieving revision 1.26
retrieving revision 1.27
diff -u -b -r1.26 -r1.27
--- helpers.h 18 Feb 2009 18:10:37 -0000 1.26
+++ helpers.h 20 Mar 2009 16:49:25 -0000 1.27
@@ -40,6 +40,16 @@
}
+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++)
Index: lapack.h
===================================================================
RCS file: /cvsroot/toon/TooN/lapack.h,v
retrieving revision 1.8
retrieving revision 1.9
diff -u -b -r1.8 -r1.9
--- lapack.h 26 Feb 2009 21:50:19 -0000 1.8
+++ lapack.h 20 Mar 2009 16:49:25 -0000 1.9
@@ -66,12 +66,12 @@
dgetrf_(M, N, A, lda, IPIV, INFO);
}
-inline void trsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M,
int* N, float* alpha, float* A, int* lda, float* B, int* ldb)
-{ strsm_(SIDE, UPLO, TRANSA, DIAG, M, N, alpha, A, lda, B, ldb);
+inline void trsm_(const char* SIDE, const char* UPLO, const char* TRANSA,
const char* DIAG, int* M, int* N, float* alpha, float* A, int* lda, float* B,
int* ldb) {
+ strsm_(const_cast<char*>(SIDE), const_cast<char*>(UPLO),
const_cast<char*>(TRANSA), const_cast<char*>(DIAG), M, N, alpha, A, lda, B,
ldb);
}
-inline void trsm_(char* SIDE, char* UPLO, char* TRANSA, char* DIAG, int* M,
int* N, double* alpha, double* A, int* lda, double* B, int* ldb) {
- dtrsm_(SIDE, UPLO, TRANSA, DIAG, M, N, alpha, A, lda, B, ldb);
+inline void trsm_(const char* SIDE, const char* UPLO, const char* TRANSA,
const char* DIAG, int* M, int* N, double* alpha, double* A, int* lda, double*
B, int* ldb) {
+ dtrsm_(const_cast<char*>(SIDE), const_cast<char*>(UPLO),
const_cast<char*>(TRANSA), const_cast<char*>(DIAG), M, N, alpha, A, lda, B,
ldb);
}
void getri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork,
int* INFO){
Index: internal/matrix.hh
===================================================================
RCS file: /cvsroot/toon/TooN/internal/matrix.hh,v
retrieving revision 1.15
retrieving revision 1.16
diff -u -b -r1.15 -r1.16
--- internal/matrix.hh 10 Mar 2009 15:30:22 -0000 1.15
+++ internal/matrix.hh 20 Mar 2009 16:49:25 -0000 1.16
@@ -60,7 +60,7 @@
SizeMismatch<Cols, Cols>::test(num_cols(), from.num_cols());
for(int r=0; r < num_rows(); r++)
- for(int c=0; c < num_rows(); c++)
+ for(int c=0; c < num_cols(); c++)
(*this)[r][c] = from[r][c];
return *this;
Index: internal/vector.hh
===================================================================
RCS file: /cvsroot/toon/TooN/internal/vector.hh,v
retrieving revision 1.23
retrieving revision 1.24
diff -u -b -r1.23 -r1.24
--- internal/vector.hh 20 Mar 2009 12:24:19 -0000 1.23
+++ internal/vector.hh 20 Mar 2009 16:49:25 -0000 1.24
@@ -103,7 +103,7 @@
Vector& operator-=(const Vector<Size2, Precision2, Base2>& rhs) {
SizeMismatch<Size,Size2>::test(size(),rhs.size());
for(int i=0; i<size(); i++)
- (*this)[i]-=rhs;
+ (*this)[i]-=rhs[i];
return *this;
}
Index: benchmark/solve_ax_equals_b.cc
===================================================================
RCS file: benchmark/solve_ax_equals_b.cc
diff -N benchmark/solve_ax_equals_b.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ benchmark/solve_ax_equals_b.cc 20 Mar 2009 16:49:25 -0000 1.1
@@ -0,0 +1,164 @@
+#include <TooN/TooN.h>
+#include <TooN/LU.h>
+#include <TooN/gaussian_elimination.h>
+#include <cvd/timer.h>
+#include <tr1/random>
+
+using namespace TooN;
+using namespace std;
+using namespace tr1;
+using namespace CVD;
+
+std::tr1::mt19937 eng;
+std::tr1::uniform_real<double> rnd;
+double global_sum;
+
+struct UseLU
+{
+ template<int R, int C> static void solve(const Matrix<R, R>& a, const
Matrix<R, C>& b, Matrix<R, C>& x)
+ {
+ LU<R> lu(a);
+
+ x = lu.backsub(b);
+ }
+
+ static string name()
+ {
+ return "LU";
+ }
+};
+
+struct UseLUInv
+{
+ template<int R, int C> static void solve(const Matrix<R, R>& a, const
Matrix<R, C>& b, Matrix<R, C>& x)
+ {
+ LU<R> lu(a);
+
+ x = lu.get_inverse() * b;
+ }
+
+ static string name()
+ {
+ return "LI";
+ }
+};
+
+struct UseGaussianElimination
+{
+ template<int R, int C> static void solve(const Matrix<R, R>& a, const
Matrix<R, C>& b, Matrix<R, C>& x)
+ {
+ x = gaussian_elimination(a, b);
+ }
+
+ static string name()
+ {
+ return "GE";
+ }
+};
+
+template<int Size, int Cols, class Solver> void benchmark_ax_eq_b()
+{
+ cvd_timer t;
+ double time=0;
+ double sum=0;
+ int n=0;
+
+ while(time < 1)
+ {
+ Matrix<Size> a;
+ for(int r=0; r < Size; r++)
+ for(int c=0; c < Size; c++)
+ a[r][c] = rnd(eng);
+
+
+ Matrix<Size, Cols> b, x;
+
+ for(int r=0; r < Size; r++)
+ for(int c=0; c < Cols; c++)
+ b[r][c] = rnd(eng);
+
+ t.reset();
+ Solver::template solve<Size, Cols>(a, b, x);
+ time += t.get_time();
+
+
+ for(int r=0; r < Size; r++)
+ for(int c=0; c < Cols; c++)
+ sum += x[r][c];
+
+ n++;
+ }
+
+ cout << Solver::name() << " " << Size << " " << Cols << " " << n / time
<< endl;
+
+ global_sum += sum;
+}
+
+template<class A, class B> struct TypeList
+{
+ typedef A a;
+ typedef B b;
+};
+
+struct Null{};
+
+template<int R, int C, class L> struct benchmark_iter
+{
+ static void iter()
+ {
+ benchmark_ax_eq_b<R, C, typename L::a>();
+ benchmark_iter<R, C, typename L::b>::iter();
+ }
+};
+
+
+template<int R, int C> struct benchmark_iter<R, C, Null>
+{
+ static void iter()
+ {
+ }
+};
+
+
+
+
+template<int Size, int Cols, class Test> struct ColIter
+{
+ static void iter()
+ {
+ benchmark_iter<Size, Cols, Test>::iter();
+ ColIter<Size, Cols-1, Test>::iter();
+ }
+};
+
+template<int Size,class Test> struct ColIter<Size, 1, Test>
+{
+ static void iter()
+ {
+ }
+};
+
+
+template<int Size, class Test> struct SizeIter
+{
+ static void iter()
+ {
+ ColIter<Size, Size*8+1, Test>::iter();
+ SizeIter<Size-1, Test>::iter();
+ }
+};
+
+template<class Test> struct SizeIter<0, Test>
+{
+ static void iter()
+ {
+ }
+};
+
+
+int main()
+{
+ SizeIter<5, TypeList<UseGaussianElimination, TypeList<UseLU, Null> >
>::iter();
+
+ return global_sum != 123456789.0;
+}
Index: test/gaussian_elimination_test.cc
===================================================================
RCS file: test/gaussian_elimination_test.cc
diff -N test/gaussian_elimination_test.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ test/gaussian_elimination_test.cc 20 Mar 2009 16:49:25 -0000 1.1
@@ -0,0 +1,50 @@
+#include <cstdlib>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <TooN/TooN.h>
+#include <TooN/helpers.h>
+#include <TooN/gaussian_elimination.h>
+#include <tr1/random>
+
+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,5> m(5,5);
+
+ for(int i=0; i< m.num_rows(); i++)
+ for(int j=0; j< m.num_rows(); j++)
+ m[i][j] = rnd(eng);
+
+ cout << m << endl;
+
+ Matrix<5,5> i;
+ Identity(i);
+
+ Matrix<5,5> inv = gaussian_elimination(m, i);
+
+ Matrix<5,5> a = m*inv;
+
+ for(int i=0; i< m.num_rows(); i++)
+ for(int j=0; j< m.num_rows(); j++)
+ {
+ if(round(a[i][j]) < 1e-10)
+ a[i][j] = 0;
+ }
+
+ cout << a;
+
+
+
+}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN gaussian_elimination.h helpers.h lapack.h ...,
Edward Rosten <=