news.swissonline.ch
2003-11-19 21:57:22 UTC
Hello
Since a week or two I'm learning functional languages. I realy like the
features like lazy evaluation and I think they are much cleaner then
imperativ languages. Today I did a little benchmark between my matrix
lib in C++ and Haskell. I was shocked ! The C++ versions needs 25
Seconds to calc the determinant of a 12x12 identity matrix. The Haskell
versions needs 1min15sec for a 11x11 identity matrix. Maybe you have
some tricks for me to speed up the haskell version. The Haskell version
is compiled with GHC (ghc -o).
Cheers Tobias
----------
C++ source
----------
#include <iostream>
namespace math
{
template<typename T, unsigned int SIZE>
class vec;
template<typename T, unsigned int ROWS, unsigned int COLS>
class matrix {
private:
T m_fields[ROWS][COLS];
void checkindex(unsigned int row, unsigned int col) const {
if(row >= ROWS) {
throw std::out_of_range("row index is out of range");
}
if(col >= COLS) {
throw std::out_of_range("col index is out of range");
}
}
public:
matrix() {
}
matrix(const matrix<T, ROWS+1, COLS+1>& c, unsigned int ai, unsigned
int ak) {
for(unsigned int i = 0, thisi = 0; i < ROWS+1; ++i) {
if(i != ai) {
for(unsigned int k = 0, thisk = 0; k < COLS+1; ++k) {
if(k != ak) {
m_fields[thisi][thisk] = c.at(i,k);
++thisk;
}
}
++thisi;
}
}
}
T& at(unsigned int row, unsigned int col) {
#ifdef _DEBUG
checkindex(row, col);
#endif
return m_fields[row][col];
}
const T& at(unsigned int row, unsigned int col) const {
#ifdef _DEBUG
checkindex(row, col);
#endif
return m_fields[row][col];
}
};
template<typename T, unsigned int ORDER>
T det(const matrix<T, ORDER, ORDER>& a)
{
T d = T(0.0);
for(unsigned int icol = 0; icol < ORDER; ++icol) {
matrix<T, ORDER-1, ORDER-1> sub(a, 0, icol);
T pre = (icol % 2) == 0 ? T(1.0) : T(-1.0);
d += pre * a.at(0,icol) * det(sub);
}
return d;
}
template<typename T>*/
T det(const matrix<T, 1, 1>& a)
{
return a.at(0,0);
}
template<typename T, unsigned int ORDER>
matrix<T, ORDER, ORDER> identity()
{
matrix<T, ORDER, ORDER> r;
for(unsigned int irow = 0; irow < ORDER; ++irow) {
for(unsigned int icol = 0; icol < ORDER; ++icol) {
r.at(irow,icol) = (icol == irow) ? T(1.0) : T(0.0);
}
}
return r;
}
int main()
{
matrix<float,12,12> mi(identity<float,12>());
int dummy;
cin >> dummy;
if(dummy == 1) {
mi.at(1,1) = 0.9f;
mi.at(1,5) = 0.9f;
mi.at(1,3) = 1.9f;
}
cout << det(mi);
}
--------------
Haskell source
--------------
import List
identity :: Int -> [[Float]]
identity 0 = []
identity (n+1) = [1.0: replicate n 0.0] ++ map (0:) (identity n)
sub :: [a] -> [[a]]
sub [x] = [[]]
sub (x:xs) = xs: (map (x:) (sub xs))
submat :: [[Float]] -> [[[[Float]]]]
submat [[x]] = [[[[]]]]
submat x = [transpose (map sub subm) | subm<-(sub x)]
det :: [[Float]] -> Float
det [[x]] = x
det x = sum (zipWith (*) (zipWith (*) detmatrix sign) (head x))
where detmatrix = [det subm | subm<-(head (submat x))]
sign = [1.0,-1.0] ++ sign
main :: IO ()
main = do {d <- return (det (identity 11)); putStr (show d)}
Since a week or two I'm learning functional languages. I realy like the
features like lazy evaluation and I think they are much cleaner then
imperativ languages. Today I did a little benchmark between my matrix
lib in C++ and Haskell. I was shocked ! The C++ versions needs 25
Seconds to calc the determinant of a 12x12 identity matrix. The Haskell
versions needs 1min15sec for a 11x11 identity matrix. Maybe you have
some tricks for me to speed up the haskell version. The Haskell version
is compiled with GHC (ghc -o).
Cheers Tobias
----------
C++ source
----------
#include <iostream>
namespace math
{
template<typename T, unsigned int SIZE>
class vec;
template<typename T, unsigned int ROWS, unsigned int COLS>
class matrix {
private:
T m_fields[ROWS][COLS];
void checkindex(unsigned int row, unsigned int col) const {
if(row >= ROWS) {
throw std::out_of_range("row index is out of range");
}
if(col >= COLS) {
throw std::out_of_range("col index is out of range");
}
}
public:
matrix() {
}
matrix(const matrix<T, ROWS+1, COLS+1>& c, unsigned int ai, unsigned
int ak) {
for(unsigned int i = 0, thisi = 0; i < ROWS+1; ++i) {
if(i != ai) {
for(unsigned int k = 0, thisk = 0; k < COLS+1; ++k) {
if(k != ak) {
m_fields[thisi][thisk] = c.at(i,k);
++thisk;
}
}
++thisi;
}
}
}
T& at(unsigned int row, unsigned int col) {
#ifdef _DEBUG
checkindex(row, col);
#endif
return m_fields[row][col];
}
const T& at(unsigned int row, unsigned int col) const {
#ifdef _DEBUG
checkindex(row, col);
#endif
return m_fields[row][col];
}
};
template<typename T, unsigned int ORDER>
T det(const matrix<T, ORDER, ORDER>& a)
{
T d = T(0.0);
for(unsigned int icol = 0; icol < ORDER; ++icol) {
matrix<T, ORDER-1, ORDER-1> sub(a, 0, icol);
T pre = (icol % 2) == 0 ? T(1.0) : T(-1.0);
d += pre * a.at(0,icol) * det(sub);
}
return d;
}
template<typename T>*/
T det(const matrix<T, 1, 1>& a)
{
return a.at(0,0);
}
template<typename T, unsigned int ORDER>
matrix<T, ORDER, ORDER> identity()
{
matrix<T, ORDER, ORDER> r;
for(unsigned int irow = 0; irow < ORDER; ++irow) {
for(unsigned int icol = 0; icol < ORDER; ++icol) {
r.at(irow,icol) = (icol == irow) ? T(1.0) : T(0.0);
}
}
return r;
}
int main()
{
matrix<float,12,12> mi(identity<float,12>());
int dummy;
cin >> dummy;
if(dummy == 1) {
mi.at(1,1) = 0.9f;
mi.at(1,5) = 0.9f;
mi.at(1,3) = 1.9f;
}
cout << det(mi);
}
--------------
Haskell source
--------------
import List
identity :: Int -> [[Float]]
identity 0 = []
identity (n+1) = [1.0: replicate n 0.0] ++ map (0:) (identity n)
sub :: [a] -> [[a]]
sub [x] = [[]]
sub (x:xs) = xs: (map (x:) (sub xs))
submat :: [[Float]] -> [[[[Float]]]]
submat [[x]] = [[[[]]]]
submat x = [transpose (map sub subm) | subm<-(sub x)]
det :: [[Float]] -> Float
det [[x]] = x
det x = sum (zipWith (*) (zipWith (*) detmatrix sign) (head x))
where detmatrix = [det subm | subm<-(head (submat x))]
sign = [1.0,-1.0] ++ sign
main :: IO ()
main = do {d <- return (det (identity 11)); putStr (show d)}