librheolef-dev 6.5-1build1 / usr / include / rheolef / tiny_matvec.h

# ifndef _RHEO_TINY_MATVEC_H
# define _RHEO_TINY_MATVEC_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//
// very small matrix - vector
//
// authors: Pierre.Saramito@imag.fr
//
// date: 7 july 1997
//
#include "rheolef/compiler.h"

namespace rheolef { 
// take a 2^n since a(i,j) -> table (tiny_size_max*j + i)
// and tiny_size_max*j == j << log2(tiny_size_max*j) is fast

const unsigned int tiny_size_max = 32;

template <class T>
class tiny_vector {
public:
    typedef typename std::vector<int>::size_type size_type;

    tiny_vector (size_type n = 0);
    tiny_vector (size_type n, const T& value);
    size_type size() const { return size_; }
    void resize(size_type n);
    const T& operator() (size_type i) const { return t_[i+i0_]; }
    const T& operator[] (size_type i) const { return t_[i+i0_]; }
    T& operator() (size_type i) { return t_[i+i0_]; }
    T& operator[] (size_type i) { return t_[i+i0_]; }
    void set_origin(size_type i) { i0_ = i; }
    size_type get_origin() const { return i0_; }
    void fill(const T& val) {
	for (size_type i = i0_; i < i0_ + size_; i++) t_ [i] = val; }
    void reset() { fill(T()); }
protected:
    T t_ [tiny_size_max];
    size_type size_;
    size_type i0_;
};
template <class T>
class tiny_matrix {
public:
    typedef typename tiny_vector<T>::size_type size_type;
    tiny_matrix (size_type nr = 0, size_type nc = 0);
    size_type nrow() const { return nrow_; }
    size_type ncol() const { return ncol_; }
    T& operator() (size_type i, size_type j) { return t_[i+i0_][j+j0_]; }
    const T& operator() (size_type i, size_type j) const { return t_[i+i0_][j+j0_]; }
    T& operator() (size_type i) { return t_[i+i0_][i+j0_]; }
    const T& operator() (size_type i) const { return t_[i+i0_][i+j0_]; }
    void set_origin(size_type i, size_type j) { i0_ = i; j0_ = j; }
    void resize(size_type nr, size_type nc);
    size_type get_row_origin() const { return i0_; }
    size_type get_col_origin() const { return j0_; }
    void fill(const T& val);
    void reset() { fill(T()); }
private:
    T t_ [tiny_size_max][tiny_size_max];
    size_type nrow_;
    size_type ncol_;
    size_type i0_;
    size_type j0_;
};
// =====================================================================
// inlined
// =====================================================================

template <class T>
inline
tiny_vector<T>::tiny_vector (size_type n)
  : size_(n), i0_(0)
{ 
    check_macro (n <= tiny_size_max, "invalid size"); 
#ifdef _RHEOLEF_PARANO
    std::fill (t_, t_+tiny_size_max, std::numeric_limits<T>::max());
#endif // _RHEOLEF_PARANO
}
template <class T>
inline
tiny_vector<T>::tiny_vector (size_type n, const T& value)
  : size_(n), i0_(0)
{ 
    check_macro (n <= tiny_size_max, "invalid size"); 
    fill (t_, t_+tiny_size_max, value);
}
template <class T>
inline
void
tiny_vector<T>::resize(size_type n)
{
    size_ = n;
    check_macro (n <= tiny_size_max, "invalid size"); 
#ifdef _RHEOLEF_PARANO
    std::fill (t_, t_+tiny_size_max, std::numeric_limits<T>::max());
#endif // _RHEOLEF_PARANO
}
template <class T>
inline
tiny_matrix<T>::tiny_matrix (size_type nr, size_type nc)
  : nrow_(nr), ncol_(nc), i0_(0), j0_(0)
{
    check_macro (nr <= tiny_size_max && nc <= tiny_size_max, "invalid sizes");
#ifdef _RHEOLEF_PARANO
    for (size_type i = 0; i < tiny_size_max; i++) 
      for (size_type j = 0; j < tiny_size_max; j++) 
	t_[i][j] = std::numeric_limits<T>::max();
#endif // _RHEOLEF_PARANO

}
template <class T>
inline
void
tiny_matrix<T>::resize(size_type nr, size_type nc)
{
    nrow_ = nr;
    ncol_ = nc;
    check_macro (nr <= tiny_size_max && nc <= tiny_size_max, "invalid sizes");
#ifdef _RHEOLEF_PARANO
    for (size_type i = 0; i < tiny_size_max; i++) 
      for (size_type j = 0; j < tiny_size_max; j++) 
	t_[i][j] = std::numeric_limits<T>::max();
#endif // _RHEOLEF_PARANO
}
template <class T>
inline
void
tiny_matrix<T>::fill(const T& val)
{
    for (size_type i = i0_; i < i0_ + nrow_; i++) 
      for (size_type j = j0_; j < j0_ + ncol_; j++) 
	t_ [i][j] = val;
}
template <class T>
void
trans(const tiny_matrix<T>& a, tiny_matrix<T>& b)
{
    typedef typename tiny_matrix<T>::size_type size_type;
    b.resize (a.ncol(), a.nrow());
    for (size_type i = 0; i < a.nrow(); i++) 
      for (size_type j = 0; j < a.ncol(); j++)
	b(j,i) = a(i,j);
}
template<class T>
tiny_matrix<T>
operator* (const tiny_matrix<T>& a, const tiny_matrix<T>& b)
 {
    check_macro(a.ncol()==b.nrow(),"Error in matrices sizes for multiplication, "
   	<< a.nrow()<<"x"<<a.ncol() <<" and "<< b.nrow()<<"x"<<b.ncol());
    typedef typename tiny_matrix<T>::size_type size_type;
    tiny_matrix<T> c(a.nrow(),b.ncol());
    c.fill(0);
    for (size_type i=0; i<a.nrow(); i++)
     for (size_type j=0; j<b.ncol(); j++)
      for (size_type k=0; k<b.nrow(); k++)
       c(i,j)+=a(i,k)*b(k,j);
    return c;
 }
template<class T>
tiny_vector<T>
operator* (const tiny_matrix<T>& a, const tiny_vector<T>& u)
 {
    check_macro(a.ncol()==u.size(),"Error in matrice-vector sizes for multiplication, "
   	<< a.nrow()<<"x"<<a.ncol() <<" and "<< u.size());
    typedef typename tiny_matrix<T>::size_type size_type;
    tiny_vector<T> v(a.nrow());
    v.fill(0);
    for (size_type i=0; i<a.nrow(); i++)
     for (size_type j=0; j<u.size(); j++)
      v(i)+=a(i,j)*u(j);
    return v;
 }
}// namespace rheolef
# endif /* _RHEO_TINY_MATVEC_H */