[ADD] WTensor and WTensorSym templates

src/common/math/WBinom.cpp

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+
+#include "WBinom.h"

src/common/math/WBinom.h

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+#ifndef WBINOM_H
+#define WBINOM_H
+
+#include <string>
+
+namespace wmath
+{
+/**
+ * Implements compile-time calculation of binomial coefficients.
+ *
+ *
+ * WBinom< n, k >::value = n! / ( k!(n-k)! ).
+ *
+ * \note For k > n, compilation fails.
+ */
+template< std::size_t n, std::size_t k >
+struct WBinom
+{
+    enum
+    {
+        value = WBinom< n - 1, k - 1 >::value + WBinom< n - 1, k >::value
+    };
+};
+
+/**
+ * Specialization for n = k.
+ */
+template< std::size_t n >
+struct WBinom< n, n >
+{
+    enum
+    {
+        value = 1
+    };
+};
+
+/**
+ * Specialization for k = 0.
+ */
+template< std::size_t n >
+struct WBinom< n, 0 >
+{
+    enum
+    {
+        value = 1
+    };
+};
+
+/**
+ * This specialization of the WBinom struct is needed to avoid
+ * infinite recursion in case of k > n. The compiler should abort
+ * compilation with an error message.
+ */
+template< std::size_t k >
+struct WBinom< 0, k >
+{
+};
+
+} // namespace wmath
+
+#endif  // WBINOM_H

src/common/math/WTensor.cpp

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+
+#include "WTensor.h"

src/common/math/WTensor.h

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+#ifndef WTENSOR_H
+#define WTENSOR_H
+
+#include <vector>
+#include <iostream>
+
+#include "WTensorBase.h"
+
+namespace wmath
+{
+// ############################# class WTensor<> #############################################
+
+/**
+ * Implements a tensor that has the same number of components in every
+ * direction.
+ *
+ * The first template parameter is the order of the tensor. 
+ * The second template parameter is the dimension of the tensor, i.e. the number of components
+ * in each direction.
+ * The third template parameter is the datatype of the components, which is double by default.
+ *
+ * \note The type Data_T may not throw exceptions on construction, destruction or
+ * during operator =.
+ *
+ * Access to specific elements of the tensor can be achieved in 2 ways:
+ *
+ * - operator (), whose number of parameters matches the template parameter order.
+ * - operator [], whose parameter is either an array or a std::vector of appropriate size.
+ *
+ * \note The datatype of the array or std::vector can be any type castable to std::size_t.
+ * \note There is no bounds checking for the array version of operator [].
+ * \note Operator () is not supported for orders larger than 6.
+ *
+ * Examples:
+ *
+ * - Construct a tensor of order 2 and dimension 3 (i.e. a 3x3-Matrix):
+ *
+ *     wmath::WTensor< 2, 3 > w;
+ *
+ * - Change Element (2,0) to 4.0:
+ *
+ *     w( 2, 0 ) = 4.0;
+ *
+ * - Construct a 4D-vector:
+ *
+ *     wmath::WTensor< 1, 4 > v;
+ *
+ * - Access v at position 2:
+ *
+ *     v( 2 ) = ...;
+ *
+ *     std::vector< int > i( 1, 2 );
+ *     v[ i ] = ...; 
+ */
+template< std::size_t order, std::size_t dim, typename Data_T = double >
+class WTensor : public WTensorFunc< WTensorBase, order, dim, Data_T >
+{
+};
+
+/**
+ * Specialization for dim = 0. This essentially prohibits
+ * instantiation of a 0-dimensional tensor.
+ */
+template< std::size_t order, typename Data_T >
+class WTensor< order, 0, Data_T >
+{
+};
+
+// ###################### special case WTensor<> of order 0 ##############################
+
+/**
+ * The order 0 version of WTensor. This essentialy encapsulates a scalar.
+ */
+template< std::size_t dim, typename Data_T >
+class WTensor< 0, dim, Data_T >
+{
+public:
+    /**
+     * Standard constructor.
+     */
+    WTensor();
+
+    /**
+     * Get the dimension of this tensor.
+     *
+     * \return The dimension of this tensor.
+     */
+    std::size_t getDimension() const;
+
+    /**
+     * Get the order of this tensor.
+     *
+     * \return The order of this tensor.
+     */
+    std::size_t getOrder() const;
+
+    /**
+     * Access via std::vector of indices. Only for compatibility.
+     *
+     * \param indices A std::vector of indices.
+     *
+     * \return A reference to the value of this tensor.
+     */
+    template< typename Index_T >
+    Data_T& operator[] ( std::vector< Index_T > const& indices );
+
+    /**
+     * Access via std::vector of indices. Only for compatibility.
+     *
+     * \param indices A std::vector of indices.
+     *
+     * \return A reference to the value of this tensor.
+     */
+    template< typename Index_T >
+    Data_T const& operator[] ( std::vector< Index_T > const& indices ) const;
+
+    /**
+     * Access operator.
+     *
+     * \return A reference to the value of this tensor.
+     */
+    Data_T& operator() ();
+
+    /**
+     * Access operator.
+     *
+     * \return A reference to the value of this tensor.
+     */
+    Data_T const& operator() () const;
+
+private:
+    /**
+     * The data element.
+     */
+    Data_T m_data;
+};
+
+template< std::size_t dim, typename Data_T >
+WTensor< 0, dim, Data_T >::WTensor()
+    : m_data( Data_T() )
+{
+}
+
+template< std::size_t dim, typename Data_T >
+std::size_t WTensor< 0, dim, Data_T >::getDimension() const
+{
+    return 0;
+}
+
+template< std::size_t dim, typename Data_T >
+std::size_t WTensor< 0, dim, Data_T >::getOrder() const
+{
+    return 0;
+}
+
+template< std::size_t dim, typename Data_T > template< typename Index_T >
+Data_T& WTensor< 0, dim, Data_T >::operator[] ( std::vector< Index_T > const& indices )
+{
+    return m_data;
+}
+
+template< std::size_t dim, typename Data_T > template< typename Index_T >
+Data_T const& WTensor< 0, dim, Data_T >::operator[] ( std::vector< Index_T > const& indices ) const
+{
+    return m_data;
+}
+
+template< std::size_t dim, typename Data_T >
+Data_T& WTensor< 0, dim, Data_T >::operator() ()
+{
+    return m_data;
+}
+
+template< std::size_t dim, typename Data_T >
+Data_T const& WTensor< 0, dim, Data_T >::operator() () const
+{
+    return m_data;
+}
+
+// ######################## stream output operators #################################
+
+/**
+ * Write a tensor of order 0 to an output stream.
+ *
+ * \param o An output stream.
+ * \param t A WTensor.
+ *
+ * \return The output stream.
+ */
+template< std::size_t dim, typename Data_T >
+std::ostream& operator << ( std::ostream& o, WTensor< 0, dim, Data_T > const& t )
+{
+    o << t() << std::endl;
+    return o;
+}
+
+/**
+ * Write a tensor of order 1 to an output stream.
+ *
+ * \param o An output stream.
+ * \param t A WTensor.
+ *
+ * \return The output stream.
+ */
+template< std::size_t dim, typename Data_T >
+std::ostream& operator << ( std::ostream& o, WTensor< 1, dim, Data_T > const& t )
+{
+    for( std::size_t k = 0; k < dim; ++k )
+    {
+        o << t( k ) << " ";
+    }
+    o << std::endl;
+    return o;
+}
+
+/**
+ * Write a tensor of order 2 to an output stream.
+ *
+ * \param o An output stream.
+ * \param t A WTensor.
+ *
+ * \return The output stream.
+ */
+template< std::size_t dim, typename Data_T >
+std::ostream& operator << ( std::ostream& o, WTensor< 2, dim, Data_T > const& t )
+{
+    for( std::size_t k = 0; k < dim; ++k )
+    {
+        for( std::size_t l = 0; l < dim; ++l )
+        {
+            o << t( k, l ) << " ";
+        }
+        o << std::endl;
+    }
+    return o;
+}
+
+}  // namespace wmath
+
+#endif  // WTENSOR_H

src/common/math/WTensorBase.cpp

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+
+#include "WTensorBase.h"

src/common/math/WTensorFunctions.cpp

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+#include <vector>
+#include <cmath>
+
+#include "../WAssert.h"
+#include "../WLimits.h"
+
+#include "WTensorFunctions.h"
+
+namespace wmath
+{
+template< typename Data_T >
+void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
+                          std::vector< Data_T >* eigenValues,
+                          std::vector< WTensor< 1, 3, Data_T > >* eigenVectors )
+{
+    WTensorSym< 2, 3, Data_T > in = WTensorSym< 2, 3, Data_T >( mat );
+    WTensor< 2, 3, Data_T > ev;
+
+    int iter = 50;
+    Data_T evp[ 3 ];
+    Data_T evq[ 3 ];
+
+    while( iter >= 0 )
+    {
+        int p = 1;
+        int q = 0;
+
+        for( int i = 0; i < 2; ++i )
+        {
+            if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
+            {
+                p = 2;
+                q = i;
+            }
+        }
+
+        if( fabs( in( p, q ) ) == 0.0 )
+        {
+            for( int i = 0; i < 3; ++i )
+            {
+                eigenValues->at( i ) = in( i, i );
+                for( int j = 0; j < 3; ++j )
+                {
+                    eigenVectors->at( i )( j ) = ev( i, j );
+                }
+            }
+            return;
+        }
+
+        Data_T r = in( q, q ) - in( p, p );
+        Data_T o = r / ( 2.0 * in( p, q ) );
+
+        Data_T t;
+        Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
+        if( o * o > wlimits::MAX_DOUBLE )
+        {
+            t = signofo / ( 2.0 * fabs( o ) );
+        }
+        else
+        {
+            t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
+        }
+
+        Data_T c;
+
+        if( t * t < wlimits::DBL_EPS )
+        {
+            c = 1.0;
+        }
+        else
+        {
+            c = 1.0 / sqrt( t * t + 1.0 );
+        }
+
+        Data_T s = c * t;
+
+        int k = 0;
+        while( k == q || k == p )
+        {
+            ++k;
+        }
+        WAssert( k < 3, "" );
+
+        Data_T u = ( 1.0 - c ) / s;
+
+        Data_T x = in( k, p );
+        Data_T y = in( k, q );
+        in( p, k ) = in( k, p ) = x - s * ( y + u * x );
+        in( q, k ) = in( k, q ) = y + s * ( x - u * y );
+        x = in( p, p );
+        y = in( q, q );
+        in( p, p ) = x - t * in( p, q );
+        in( q, q ) = y + t * in( p, q );
+        in( q, p ) = in( p, q ) = 0.0;
+
+        for( int l = 0; l < 3; ++l )
+        {
+            evp[ l ] = ev( l, p );
+            evq[ l ] = ev( l, q );
+        }
+        for( int l = 0; l < 3; ++l )
+        {
+            ev( l, p ) = c * evp[ l ] - s * evq[ l ];
+            ev( l, q ) = s * evp[ l ] + c * evq[ l ];
+        }
+
+        --iter;
+    }
+    WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
+}
+}  // namespace wmath

src/common/math/WTensorFunctions.h

+//---------------------------------------------------------------------------
+//
+// Project: OpenWalnut ( http://www.openwalnut.org )
+//
+// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
+// For more information see http://www.openwalnut.org/copying
+//
+// This file is part of OpenWalnut.
+//
+// OpenWalnut is free software: you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// OpenWalnut is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
+//
+//---------------------------------------------------------------------------
+
+#ifndef WTENSORFUNCTIONS_H
+#define WTENSORFUNCTIONS_H
+
+#include <vector>
+
+#include "WTensorSym.h"
+
+namespace wmath
+{
+/**
+ * Compute all eigenvalues as well as the corresponding eigenvectors of a
+ * symmetric real Matrix.
+ *
+ * \pre Data_T must be castable to double.
+ *
+ * \param[in] mat A real symmetric matrix.
+ * \param[out] eigenValues A pointer to a vector of eigenvalues.
+ * \param[out] eigenVectors A pointer to a vector of eigenvectors.
+ */
+template< typename Data_T >
+void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
+                          std::vector< Data_T >* eigenValues,
+                          std::vector< WTensor< 1, 3, Data_T > >* eigenVectors );
+
+} // namespace wmath
+
+#endif  // WTENSORFUNCTIONS_H

src/common/math/WTensorSym.cpp

+//---------------------------------------------------------------------------