TinyMatrix.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program 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, or (at your option)
//  any later version.

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

//  $Id: TinyMatrix.hpp,v 1.15 2007/11/02 20:23:16 delpinux Exp $

#ifndef TINY_MATRIX_HPP
#define TINY_MATRIX_HPP

#include <TinyVector.hpp>

template <size_t M, size_t N,
          typename T=real_t>
class TinyMatrix
{
public:
  typedef TinyVector<N>
  AssociatedVector;             
  typedef TinyVector<M>
  AssociatedTransposedVector;   
  typedef T ValueType;          
  enum {
    NumberOfLines = M,          
    NumberOfRow   = N           
  };

private:
  T __x[M][N];                  
public:
  
  inline T& operator()(const size_t& i, const size_t& j)
  {
    ASSERT ((i<M) && (j<N));
    return __x[i][j];
  }

  inline const T& operator()(const size_t& i, const size_t& j) const
  {
    ASSERT ((i<M) && (j<N));
    return __x[i][j];
  }

  inline const TinyMatrix<M,N,T>& operator=(const TinyMatrix<M,N,T>& B)
  {
    T* y = __x[0];
    const T* B_y = B.__x[0];
    for (size_t i=0; i<M*N; i++)
      y[i] = B_y[i];
    return (*this);
  }

  inline const TinyMatrix<M,N,T>& operator += (const TinyMatrix<M,N,T>& B)
  {
    T* y = __x[0];
    const T* B_y = B.__x[0];
    for (size_t i=0; i<M*N; i++)
      y[i] += B_y[i];
    return (*this);
  }

  inline const TinyMatrix<M,N,T>& operator -= (const TinyMatrix<M,N,T>& B)
  {
    T* y = __x[0];
    const T* B_y = B.__x[0];
    for (size_t i=0; i<M*N; i++)
      y[i] -= B_y[i];
    return (*this);
  }

  template <size_t P>
  inline TinyMatrix<M,P> operator * (const TinyMatrix<N,P>& B) const
  {
    TinyMatrix<M,P> C(0);
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<P; j++)
        for (size_t k=0; k<N; k++)
          C(i,j) += __x[i][k] * B(k,j);
    return C;
  }

  inline TinyVector<M,T> operator * (const TinyVector<N,T>& v) const
  {
    TinyVector<M,T> u = 0;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
          u[i] += __x[i][j] * v[j];
    return u;
  }

  inline TinyMatrix<M,N,T> operator * (const T& t) const
  {
    TinyMatrix<M,N,T> B;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
          B.__x[i][j] = __x[i][j] * t;
    return B;
  }

  inline TinyMatrix<M,N,T> operator / (const T& t) const
  {
    ASSERT(t != 0);
    TinyMatrix<M,N,T> B;
    const T t_1 = 1./t;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
          B.__x[i][j] = __x[i][j] * t_1;
    return B;
  }

  inline const TinyMatrix<M,N,T>& operator *= (const T& t)
  {
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] *= t;
    return (*this);
  }

  inline const TinyMatrix<M,N,T>& operator /= (const T& t)
  {
    const T t_1 = 1./t;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] *= t_1;
    return (*this);
  }

  inline friend TinyMatrix<M,N,T> operator * (const T& t,
                                              const TinyMatrix<M,N,T>& A)
  {
    return (A*t);
  }

  inline TinyMatrix<M,N,T> operator + (const TinyMatrix<M,N,T>& B) const
  {
    TinyMatrix<M,N,T> C;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
        C.__x[i][j] = __x[i][j] + B.__x[i][j];
    return C;
  }

  inline TinyMatrix<M,N,T> operator - (const TinyMatrix<M,N,T>& B) const
  {
    TinyMatrix<M,N,T> C;
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
          C.__x[i][j] = __x[i][j] - B.__x[i][j];
    return C;
  }

  friend std::ostream& operator << (std::ostream& os,
                                    const TinyMatrix<M,N,T>& A)
  {
    for (size_t i=0; i<M; i++) {
      for (size_t j=0; j<N; j++)
        os << A(i,j) << " ";
      os << '\n';
    }
    return os;
  }

  TinyMatrix()
  {
    T* y = __x[0];
    for (size_t i=0; i<M*N; i++) {
      y[i] = 0;
    }
  }

  TinyMatrix(const TinyMatrix<M,N,T>& B)
  {
    T* y = __x[0];
    const T* B_y = B.__x[0];
    for (size_t i=0; i<M*N; i++) {
      y[i] = B_y[i];
    }
  }

  TinyMatrix(const T& t)
  {
    for (size_t i=0; i<M; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = (i==j)?t:0;
  }

  ~TinyMatrix()
  {
    ;
  }
};

template <size_t N,
          typename T>
class TinyMatrix<N,N,T>
{
public:
  typedef TinyVector<N>
  AssociatedVector;             
  typedef TinyVector<N> 
  AssociatedTransposedVector;   
  typedef T ValueType;          
  enum {
    NumberOfLines = N,          
    NumberOfRow   = N           
  };

private:
  T __x[N][N];                  
public:

  size_t size() const
  {
    return N;
  }

  inline T& operator()(const size_t& i, const size_t& j)
  {
    ASSERT ((i<N) && (j<N));
    return __x[i][j];
  }

  inline const T& operator()(const size_t& i, const size_t& j) const
  {
    ASSERT ((i<N) && (j<N));
    return __x[i][j];
  }

  inline TinyMatrix<N,N,T> invert() const
  {
    TinyMatrix<N,N,T> M;
    TinyMatrix<N,N,T> lu;
    TinyMatrix<N,N,T> P;
    LU(*this, lu, P);

    TinyVector<N,T> v;
    TinyVector<N,T> u;

    for (size_t i=0; i<N; ++i) {
      for (size_t j=0; j<N; ++j)
        v[j] = (i==j)?1:0;

      u = lu.applyLU(P*v);
      for (size_t j=0; j<N; ++j)
        M(j,i) = u[j];
    }

    return M;
  }

  inline TinyMatrix<N,N,T>& operator = (const TinyMatrix<N,N,T>& B)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = B.__x[i][j];
    return (*this);
  }

  inline const TinyMatrix<N,N,T>& operator += (const TinyMatrix<N,N,T>& B)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] += B.__x[i][j];
    return (*this);
  }

  inline const TinyMatrix<N,N,T>& operator -= (const TinyMatrix<N,N,T>& B)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] -= B.__x[i][j];
    return (*this);
  }

  inline const TinyMatrix<N,N,T>& operator = (const T& t)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = (i==j)?t:0;
    return (*this);
  }

  template <size_t P>
  inline TinyMatrix<N,P> operator * (const TinyMatrix<N,P>& B) const
  {
    TinyMatrix<N,P> C(0);
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<P; j++)
        for (size_t k=0; k<N; k++)
          C(i,j) += __x[i][k] * B(k,j);
    return C;
  }

  inline TinyVector<N,T> operator * (const TinyVector<N,T>& v) const
  {
    TinyVector<N,T> u = 0;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
          u[i] += __x[i][j] * v[j];
    return u;
  }

  inline TinyMatrix<N,N,T> operator * (const T& t) const
  {
    TinyMatrix<N,N,T> B;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
          B.__x[i][j] = __x[i][j] * t;
    return B;
  }

  inline TinyMatrix<N,N,T> operator / (const T& t) const
  {
    ASSERT (t!=0);
    const T t_1 = 1./t;
    TinyMatrix<N,N,T> B;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
          B.__x[i][j] = __x[i][j] * t_1;
    return B;
  }

  inline const TinyMatrix<N,N,T>& operator *= (const T& t)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] *= t;
    return (*this);
  }

  inline TinyMatrix<N,N,T>& operator /= (const T& t) {
    const T t_1 = 1./t;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] *= t_1;
    return (*this);
  }

  inline TinyMatrix<N,N,T> operator + (const TinyMatrix<N,N,T>& B) const
  {
    TinyMatrix<N,N,T> C;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
          C.__x[i][j] = __x[i][j] + B.__x[i][j];
    return C;
  }

  inline TinyMatrix<N,N,T> operator - (const TinyMatrix<N,N,T>& B) const
  {
    TinyMatrix<N,N,T> C;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
          C.__x[i][j] = __x[i][j] - B.__x[i][j];
    return C;
  }

  inline friend TinyMatrix<N,N,T> operator * (const T& t,
                                              const TinyMatrix<N,N,T>& A)
  {
    return (A*t);
  }

  inline TinyMatrix<N,N,T> operator - () const
  {
    TinyMatrix<N,N,T> B;
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        B.__x[i][j] = - __x[i][j];
    return B;
  }

  friend std::ostream& operator << (std::ostream& os,
                                    const TinyMatrix<N,N,T>& A)
  {
    for (size_t i=0; i<N; i++) {
      for (size_t j=0; j<N; j++)
        os << A(i,j) << " ";
      os << '\n';
    }
    return os;
  }

  friend T det(const TinyMatrix<N,N,T>& A)
  {
    T d = 1.;
    TinyMatrix<N,N,T> B = A;

    for (size_t k=0; k<N-1; k++) {
      size_t maxLine = k;
      for (size_t i=k+1; i<N; i++) {
        maxLine = (std::abs(B(i,k))>std::abs(B(maxLine,k)))?i:maxLine;
      }
      if (k != maxLine) {
        B.swapLines(k,maxLine);
        d *= -1;
      }

      for(size_t i=k+1; i<N; ++i) {
        for (size_t j=k+1; j<N; j++) {
          B(i,j) -= B(i,k)*B(k,j)/B(k,k);
        }
        B(i,k) = 0;
      }
    }
    for (size_t i=0; i<N; ++i)
      d *= B(i,i);
    return d;
  }

  template <typename T2>
  friend T gaussPivot(const TinyMatrix<N,N,T>& A,
                      const TinyVector<N,T2>& rightHandSide,
                      TinyVector<N,T2>& X)
  {
    T t = 1.;
    TinyMatrix<N,N,T> M = A;
    TinyVector<N,T2> b = rightHandSide;

    for (size_t k=0; k<N-1; k++) {
      size_t maxLine = k;
      for (size_t i=k+1; i<N; i++) {
        maxLine = (std::abs(M(i,k))>std::abs(M(maxLine,k)))?i:maxLine;
      }
      if (k != maxLine) {
        M.swapLines(k,maxLine);
        std::swap(b[k], b[maxLine]);
        t *= -1;
      }

      for(size_t i=k+1; i<N; ++i) {
        const real_t coef = M(i,k)/M(k,k);
        for (size_t j=k; j<N; j++) {
          M(i,j) -= coef * M(k,j);
        }
        b[i] -= coef*b[k];
      }
    }
    for (size_t i=0; i<N; ++i)
      t *= M(i,i);

    for (size_t i = N-1; i < N; --i) {
      X[i] = 0;
      for (size_t j = i+1; j<N; j++) {
        X[i] -= M(i,j)*X[j];
      }
      X[i] += b[i];
      X[i] /= M(i,i);
    }
    return t;
  }

  void swapLines(const size_t& n, const size_t& m)
  {
    if (n==m)
      return;

    T temp;
    for (size_t i=0; i<N; ++i) {
      temp = __x[n][i];
      __x[n][i] = __x[m][i];
      __x[m][i] = temp;
    }
  }

  friend void LU(const TinyMatrix<N,N,T>& A,
                 TinyMatrix<N,N,T>& lu,
                 TinyMatrix<N,N,T>& P)
  {
    // Reinitialize P to Identity
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++) {
        P(i,j) = (i==j);
      }

    lu = A;
    
    for (size_t k=0; k<N-1; k++) {
      // Checkes for Pivot.
      size_t maxLine = k;
      for (size_t i=k+1; i<N; i++) {
        maxLine = (std::abs(lu(i,k))<std::abs(lu(maxLine,k)))?maxLine:i;
      }
      lu.swapLines(k,maxLine);
      P.swapLines(k,maxLine);

      for (size_t i=k+1; i<N; i++) {
        lu(i,k) /= lu(k,k);
        for (size_t j=k+1; j<N; j++)
          lu(i,j) -=lu(i,k)*lu(k,j);
      }
    }
  }

  inline TinyVector<N,T> applyLU(const TinyVector<N,T>& b) const
  {
    const TinyMatrix<N,N,T>& lu = (*this);

    TinyVector<N,T> y = b;
    // Go down
    for (size_t i=0; i<N; i++) {
      for (size_t j=0; j<i; j++)
        y[i] -= y[j]*lu(i,j);
    }

    // go up
    for (size_t i=N-1; i<N; i--) {
      for (size_t j=i+1; j<N; j++)
        y[i] -= y[j]*lu(i,j);

      y[i] /= lu(i,i);
    }
        
    return y;
  }

  friend TinyVector<N,T> operator/(const TinyVector<N,T>& b,
                                   const TinyMatrix<N,N,T>& A)
  {
    TinyMatrix<N,N,T> lu;
    TinyMatrix<N,N,T> P;
    LU(A,lu,P);

    return lu.applyLU(P*b);
  }

  TinyMatrix()
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = 0;
  }

  TinyMatrix(const TinyMatrix<N,N,T>& B)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = B.__x[i][j];
  }

  TinyMatrix(const T& t)
  {
    for (size_t i=0; i<N; i++)
      for (size_t j=0; j<N; j++)
        __x[i][j] = (i==j)?t:0;
  }

  ~TinyMatrix()
  {
    ;
  }
};


template<typename T>
class TinyMatrix<1,1,T>
{
public:
  typedef T
  AssociatedVector;             
  typedef T
  AssociatedTransposedVector;   
  typedef T ValueType;          
  enum {
    NumberOfLines = 1,          
    NumberOfRow   = 1           
  };

private:
  T __x;                        
public:
  inline operator T&()
  {
    return __x;
  }

  inline operator const T&() const
  {
    return __x;
  }

  T& operator()(const size_t& i, const size_t& j)
  {
    ASSERT ((i==0)&&(j==0));
    return __x;
  }

  const T& operator()(const size_t& i, const size_t& j) const
  {
    ASSERT ((i==0)&&(j==0));
    return __x;
  }

  friend const T& det(const TinyMatrix<1,1,T>& A)
  {
    return A.__x;
  }

  TinyMatrix()
    : __x(0)
  {
    ;
  }

  TinyMatrix(const T& t)
    : __x(t)
  {
    ;
  }

  TinyMatrix(const TinyMatrix<1,1,T>& B)
    : __x(B.__x)
  {
    ;
  }

  ~TinyMatrix()
  {
    ;
  }  
};

#endif // TINY_MATRIX_HPP

---