---

more/num/tensin.h

Go to the documentation of this file.

//  more/tensin.h -- inlined vectors, matrices and tensors
//  Copyright (C) 1998--2001  Petter Urkedal (petter.urkedal@matfys.lth.se)

//  This file 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.

//  This file 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 this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.

//  $Id: tensin.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $



#ifndef MORE_TENSIN_H
#define MORE_TENSIN_H

#include <iosfwd>
#include <more/math/math.h>
#include <more/num/vectin.h>

#ifdef MORE_CHECK_RANGE
#  include <assert.h>
#  define MORE_CHECKRNG(i, N) assert(0 <= i && i < N)
#else
#  define MORE_CHECKRNG(i, N)
#endif

namespace more {
namespace num {

  template <int I> struct index {};

  template <typename T, int N, int M>
    class matrin {
        T data[M][N]; // use fortran-compatible storage
      public:
        static const int CTC_size = M*N, CTC_size0 = N, CTC_size1 = M;
        typedef T scalar_type;
        T& operator()(int i, int j) { return data[j][i]; }
        T operator()(int i, int j) const { return data[j][i]; }
        template<typename U1, typename U2>
          T& operator()(U1 const& i, U2 const& j)
        { return (*this)(enumerator(i), enumerator(j)); }
        template<typename U1, typename U2>
          T operator()(U1 const& i, U2 const& j) const
        { return (*this)(enumerator(i), enumerator(j)); }
        static int size(int i) { return i? M : N; }
        static int size(index<0>) { return N; }
        static int size(index<1>) { return M; }

        typedef T* storage_iterator;
        typedef const T* const_storage_iterator;
        storage_iterator storage_begin() { return &data[0][0]; }
        const_storage_iterator storage_begin() const { return &data[0][0]; }
        storage_iterator storage_end() { return &data[M-1][N-1]+1; }
        const_storage_iterator storage_end() const
        { return &data[M-1][N-1]+1; }

        matrin& operator+=(const matrin& rhs) {
            for(int i = 0; i < N; ++i) for(int j = 0; j < M; ++j)
                data[i][j] += rhs(i, j);
            return *this;
        }
        matrin& operator-=(const matrin& rhs) {
            for(int i = 0; i < N; ++i) for(int j = 0; j < M; ++j)
                data[i][j] += rhs(i, j);
            return *this;
        }
    };

  template<typename T, int N, int M>
    inline matrin<T, N, M>
    operator-(const matrin<T, N, M>& A) {
        matrin<T, N, M> C;
        for(int i = 0; i < N; ++i) for(int j = 0; j < M; ++j)
            C(i,j) = -A(i,j);
        return C;
    }

  template<typename T, int N, int M>
    inline matrin<T, N, M>
    operator+(const matrin<T, N, M>& A, const matrin<T, N, M>& B) {
        matrin<T, N, M> C;
        for(int i = 0; i < N; ++i) for(int j = 0; j < M; ++j)
            C(i,j) = A(i,j) + B(i,j);
        return C;
    }

  template<typename T, int N, int M>
    inline matrin<T, N, M>
    operator-(const matrin<T, N, M>& A, const matrin<T, N, M>& B) {
        matrin<T, N, M> C;
        for(int i = 0; i < N; ++i) for(int j = 0; j < M; ++j)
            C(i,j) = A(i,j) - B(i,j);
        return C;
    }

  template<typename T, int N, int M, int L>
    inline matrin<T, N, L>
    operator*(const matrin<T, N, M>& A, const matrin<T, M, L>& B) {
        matrin<T, N, L> C;
        for(int i = 0; i < N; ++i)
            for(int j = 0; j < L; ++j) {
                T sum = 0;
                for(int k = 0; k < M; ++k) sum += A(i,k)*B(k,j);
                C(i,j) = sum;
            }
        return C;
    }


  ///\if bits
  namespace bits {
    template<int M0, int M1, int M2, int M3, int M4, int M5>
      struct rank_finder
      {
          static const int rank =
          1+rank_finder<M1, M2, M3, M4, M5, 1>::rank;
      };
    template<>
      struct rank_finder<1, 1, 1, 1, 1, 1>
      {
          static const int rank = 0;
      };
  }
  ///\endif

  template<typename T, int N0 = 1, int N1 = 1, int N2 = 1, int N3 = 1,
           int N4 = 1, int N5 = 1>
    class tensin {
        T data[N5][N4][N3][N2][N1][N0];

      public:
        T& operator()(int i0 = 0, int i1 = 0, int i2 = 0, int i3 = 0,
                      int i4 = 0, int i5 = 0)
        { return data[i5][i4][i3][i2][i1][i0]; }
        const T& operator()(int i0 = 0, int i1 = 0, int i2 = 0, int i3 = 0,
                            int i4 = 0, int i5 = 0) const
        { return data[i5][i4][i3][i2][i1][i0]; }
        static int size() { return N0*N1*N2*N3*N4*N5; }
        static int size(int i) {
            return (i==0? N0 : (i==1? N1 :
                                (i==2? N2 : (i==3? N3 :
                                             (i==4? N4 : N5)))));
        }
        static int size(index<0>) { return N0; }
        static int size(index<1>) { return N1; }
        static int size(index<2>) { return N2; }
        static int size(index<3>) { return N3; }
        static int size(index<4>) { return N4; }
        static int size(index<5>) { return N5; }

        typedef T* iterator;
        typedef const T* const_iterator;
        typedef T scalar_type;
        iterator begin() { return &data[0][0][0][0][0][0]; }
        iterator end() {
            return &data[N5-1][N4-1][N3-1][N2-1][N1-1][N0-1]+1;
        }
        const_iterator begin() const { return &data[0][0][0][0][0][0]; }
        const_iterator end() const {
            return &data[N5-1][N4-1][N3-1][N2-1][N1-1][N0-1]+1;
        }
        static int rank() {
            return bits::rank_finder<N0, N1, N2, N3, N4, N5>::rank;
        }
    };

  template<typename T, int N, int M>
    inline vectin<T, N>
    operator*(const matrin<T, N, M>& A, const vectin<T, M>& x) {
        vectin<T, N> y;
        for(int i = 0; i < N; ++i) {
            T sum = 0;
            for(int j = 0; j < M; ++j) sum += A(i,j)*x(j);
            y(i) = sum;
        }
        return y;
    }
}} // more::num

#undef MORE_CHECKRNG

#endif

---