SquareMatrix.h

Go to the documentation of this file.

// -*- C++ -*-

#if !defined(__SquareMatrix_h__)
#define __SquareMatrix_h__

#if defined(DEBUG_ads) && !defined(DEBUG_SquareMatrix)
#define DEBUG_SquareMatrix
#endif

#include "TensorTypes.h"

#include "../array/FixedArray.h"

#include "../../third-party/loki/static_check.h"

#include <iosfwd>
#include <algorithm>
#include <functional>
#include <numeric>

#include <cassert>
#include <cmath>

BEGIN_NAMESPACE_ADS

//
// NxN matrices.
//


template<int N, typename T = double>
class SquareMatrix :
  public TensorTypes<T> {
private:

  //
  // Private types.
  //

  typedef TensorTypes<T> base_type;

public:

  //
  // public types
  //

  typedef typename base_type::value_type value_type;
  typedef typename base_type::pointer pointer;
  typedef typename base_type::const_pointer const_pointer;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::index_type index_type;

protected:

  //
  // Data
  //

  value_type _elem[N][N];

public:

  //
  // Constructors
  //

  SquareMatrix() {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  }

  ~SquareMatrix() {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  }

  SquareMatrix(const SquareMatrix& x);

  template<typename T2>
  SquareMatrix(const SquareMatrix<N,T2>& x);

  SquareMatrix(const_pointer x);

  SquareMatrix(const value_type x);

  //
  // Assignment operators
  //

  SquareMatrix& 
  operator=(const SquareMatrix& x);

  SquareMatrix& 
  operator=(const value_type x);

  template<typename T2> 
  SquareMatrix&
  operator=(const SquareMatrix<N,T2>& x);

  //
  // Accessors and Manipulators
  //

  iterator
  begin() { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return &(_elem[0][0]);
  }

  iterator 
  end() { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return begin() + N * N;
  }

  const_iterator 
  begin() const { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return &(_elem[0][0]);
  }

  const_iterator 
  end() const { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return begin() + N * N;
  }

  pointer
  data() { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return &(_elem[0][0]);
  }

  const_pointer
  data() const { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return &(_elem[0][0]);
  }

  size_type
  size() const { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
    return N * N;
  }

  value_type
  operator()(const index_type i) const {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < size());
#endif
    return *(begin() + i);
  }

  reference
  operator()(const index_type i) { 
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < size());
#endif
    return *(begin() + i);
  }

  value_type 
  operator[](const index_type i) const {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < size());
#endif
    return *(begin() + i);
  }

  value_type& 
  operator[](const index_type i) {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < size());
#endif
    return *(begin() + i);
  }

  value_type
  operator()(const index_type i, const index_type j) const {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < N && 0 <= j && j < N);
#endif
    return _elem[i][j];
  }

  reference
  operator()(const index_type i, const index_type j) {
    LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < N && 0 <= j && j < N);
#endif
    return _elem[i][j];
  }

  void
  get(pointer x) const;

  void
  set(const_pointer x);


  void
  getMinor(const int i, const int j, SquareMatrix<N-1,T>& minor) const;

  void
  negate();

  void
  transpose();

  //
  // Assignment operators with scalar operand.
  //

  SquareMatrix& 
  operator+=(const value_type x);

  SquareMatrix& 
  operator-=(const value_type x);

  SquareMatrix& 
  operator*=(const value_type x);

  SquareMatrix& 
  operator/=(const value_type x);

  SquareMatrix& 
  operator%=(const value_type x);

  //
  // Assignment operators with SquareMatrix operand
  //

  template<typename T2>
  SquareMatrix& 
  operator+=(const SquareMatrix<N,T2> & x);

  template<typename T2>
  SquareMatrix& 
  operator-=(const SquareMatrix<N,T2> & x);

  template<typename T2>
  SquareMatrix& 
  operator*=(const SquareMatrix<N,T2> & x);

  //
  // Unary operators
  //

  SquareMatrix
  operator+() {
    return *this;
  }
    
  SquareMatrix
  operator-();

};

//
// Binary operators
//
    
template<typename T, int N>
inline
SquareMatrix<N,T>
operator+(const SquareMatrix<N,T>& m, const T x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(m);
  result += x;
  return result;
}

template<typename T, int N>
inline
SquareMatrix<N,T>
operator+(const T x, const SquareMatrix<N,T>& m) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return m + x;
}

template<typename T, int N>
inline
SquareMatrix<N,T>
operator+(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(x);
  result += y;
  return result;
}

template<typename T, int N>
inline
SquareMatrix<N,T>
operator-(const SquareMatrix<N,T>& m, const T x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(m);
  result -= x;
  return result;
}

template<typename T, int N>
inline
SquareMatrix<N,T>
operator-(const T x, const SquareMatrix<N,T>& m) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(x);
  result -= m;
  return result;
}

template<typename T, int N>
inline
SquareMatrix<N,T>
operator-(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(x);
  result -= y;
  return result;
}

template<int N, typename T>
inline
SquareMatrix<N,T>
operator*(const SquareMatrix<N,T>& m, const T x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  SquareMatrix<N,T> result(m);
  result *= x;
  return result;
}

template<int N, typename T>
inline
SquareMatrix<N,T>
operator*(const T x, const SquareMatrix<N,T>& m) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return m * x;
}

template<int N, typename T>
SquareMatrix<N,T>
operator*(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y);

template<int N, typename T>
inline
SquareMatrix<N,T>
operator/(const SquareMatrix<N,T>& m, const T x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
#ifdef DEBUG_SquareMatrix
  assert(x != 0);
#endif
  SquareMatrix<N,T> result(m);
  result /= x;
  return result;
}

template<int N, typename T>
SquareMatrix<N,T>
operator/(const T x, const SquareMatrix<N,T>& m);

//
// Math operators.
//

template<int N, typename T>
inline
T
computeSum(const SquareMatrix<N,T>& x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return std::accumulate(x.begin(), x.end(), T(0));
}
    
template<int N, typename T>
inline
T
computeProduct(const SquareMatrix<N,T>& x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return std::accumulate(x.begin(), x.end(), T(1), std::multiplies<T>());
}

template<int N, typename T>
inline
T
computeMinimum(const SquareMatrix<N,T>& x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return *std::min_element(x.begin(), x.end());
}
    
template<int N, typename T>
inline
T
computeMaximum(const SquareMatrix<N,T>& x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return *std::max_element(x.begin(), x.end());
}

template<int N, typename T>
T
computeDeterminant(const SquareMatrix<N,T>& x);

template<int N, typename T>
T
computeTrace(const SquareMatrix<N,T>& x);

template<int N, typename T>
SquareMatrix<N,T>
computeTranspose(const SquareMatrix<N,T>& x);

template<int N, typename T>
SquareMatrix<N,T>
computeInverse(const SquareMatrix<N,T>& x);

template<int N, typename T>
void
computeInverse(const SquareMatrix<N,T>& x, SquareMatrix<N,T>* y);

template<int N, typename T>
void
computeScaledInverse(const SquareMatrix<N,T>& x, SquareMatrix<N,T>* si);

template<int N, typename T>
SquareMatrix<N,T>
computeScaledInverse(const SquareMatrix<N,T>& x);

template<int N, typename T>
T
computeFrobeniusNormSquared(const SquareMatrix<N,T>& x);

template<int N, typename T>
inline
T
computeFrobeniusNorm(const SquareMatrix<N,T>& x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return std::sqrt(computeFrobeniusNormSquared(x));
}

template<int N, typename T>
inline
T
computeInnerProduct(const SquareMatrix<N,T>& x, const SquareMatrix<N,T>& y) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  T result = 0;
  const typename SquareMatrix<N,T>::const_iterator finish = x.end();
  for (typename SquareMatrix<N,T>::const_iterator i = x.begin(), 
         j = y.begin(); i != finish; ++i, ++j) {
    result += *i * *j;
  }
  return result;
}


template<int N, typename T>
inline
void
computeProduct(const SquareMatrix<N,T>& m, const FixedArray<N,T>& v, 
               FixedArray<N,T>* x) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  *x = 0;
  typename SquareMatrix<N,T>::const_iterator mi = x->begin();
  typename FixedArray<N,T>::const_iterator vi;
  const typename FixedArray<N,T>::const_iterator vi_end = v.end();
  // Loop over the rows.
  for (int m = 0; m != N; ++m) {
    // Loop over the columns.
    for (vi = v.begin(); vi != vi_end; ++vi, ++mi) {
      (*x)[m] += *mi * *vi;
    }
  }
}


//
// Equality
//

template<int N, typename T1, typename T2>
inline
bool
operator==(const SquareMatrix<N,T1>& a, const SquareMatrix<N,T2>& b) {
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return std::equal(a.begin(), a.end(), b.begin());
}

template<int N, typename T1, typename T2>
inline
bool
operator!=(const SquareMatrix<N,T1>& a, const SquareMatrix<N,T2>& b) { 
  LOKI_STATIC_CHECK(N > 3, N_greater_than_3);
  return !(a == b); 
}

//
// I/O
//
    
template<int N, typename T>
std::ostream&
operator<<(std::ostream& out, const SquareMatrix<N,T>& x);
    
template<int N, typename T>
std::istream&
operator>>(std::istream& in, SquareMatrix<N,T>& x);




//----------------------------------------------------------------------------
// 1x1 matrices.
//----------------------------------------------------------------------------


template<typename T>
class SquareMatrix<1,T> :
  public TensorTypes<T> {
private:

  //
  // private types
  //

  typedef TensorTypes<T> base_type;

public:

  //
  // public types
  //

  typedef typename base_type::value_type value_type;
  typedef typename base_type::pointer pointer;
  typedef typename base_type::const_pointer const_pointer;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::index_type index_type;

protected:

  //
  // Data
  //

  value_type _elem;

public:

  //
  // Constructors
  //

  SquareMatrix()
  {}

  ~SquareMatrix() 
  {}

  SquareMatrix(const SquareMatrix& x);

  template<typename T2>
  SquareMatrix(const SquareMatrix<1,T2>& x);

  SquareMatrix(const value_type e);

  SquareMatrix(const_pointer x);

  //
  // Assignment operators
  //

  SquareMatrix& 
  operator=(const SquareMatrix& x);

  SquareMatrix& 
  operator=(const value_type x);

  template<typename T2> 
  SquareMatrix&
  operator=(const SquareMatrix<1,T2>& x);

  //
  // Accessors and Manipulators
  //

  iterator
  begin() { 
    return &_elem; 
  }

  iterator 
  end() { 
    return &_elem + 1;
  }

  const_iterator 
  begin() const { 
    return &_elem;
  }

  const_iterator 
  end() const { 
    return &_elem + 1;
  }

  pointer
  data() { 
    return &_elem;
  }

  const_pointer
  data() const { 
    return &_elem; 
  }

  size_type
  size() const { 
    return 1;
  }

  value_type
  operator()(const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(i == 0);
#endif
    return _elem;
  }

  reference
  operator()(const index_type i) { 
#ifdef DEBUG_SquareMatrix
    assert(i == 0);
#endif
    return _elem;
  }

  value_type 
  operator[](const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(i == 0);
#endif
    return _elem;
  }

  value_type& 
  operator[](const index_type i) {
#ifdef DEBUG_SquareMatrix
    assert(i == 0);
#endif
    return _elem;
  }

  value_type
  operator()(const index_type i, const index_type j) const {
#ifdef DEBUG_SquareMatrix
    assert(i == 0 && j == 0);
#endif
    return _elem;
  }

  reference
  operator()(const index_type i, const index_type j) {
#ifdef DEBUG_SquareMatrix
    assert(i == 0 && j == 0);
#endif
    return _elem;
  }

  void
  get(pointer x) const;

  void
  set(const_pointer x);

  void
  get(reference e) const;

  void
  set(const value_type e);

  void
  negate();

  void
  transpose();

  //
  // Assignment operators with scalar operand.
  //

  SquareMatrix& 
  operator+=(const value_type x);

  SquareMatrix& 
  operator-=(const value_type x);

  SquareMatrix& 
  operator*=(const value_type x);

  SquareMatrix& 
  operator/=(const value_type x);

  SquareMatrix& 
  operator%=(const value_type x);

  //
  // Assignment operators with SquareMatrix operand
  //

  template<typename T2>
  SquareMatrix& 
  operator+=(const SquareMatrix<1,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator-=(const SquareMatrix<1,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator*=(const SquareMatrix<1,T2>& x);

  //
  // Unary operators
  //

  SquareMatrix
  operator+() {
    return *this;
  }
    
  SquareMatrix
  operator-();

};

//
// Binary operators
//
    
template<typename T>
inline
SquareMatrix<1,T>
operator+(const SquareMatrix<1,T>& m, const T x) {
  return SquareMatrix<1,T>(m[0] + x);
}

template<typename T>
inline
SquareMatrix<1,T>
operator+(const T x, const SquareMatrix<1,T>& m) {
  return m + x;
}

template<typename T>
inline
SquareMatrix<1,T>
operator+(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) {
  return SquareMatrix<1,T>(x[0] + y[0]);
}

template<typename T>
inline
SquareMatrix<1,T>
operator-(const SquareMatrix<1,T>& m, const T x) {
  return SquareMatrix<1,T>(m[0] - x);
}

template<typename T>
inline
SquareMatrix<1,T>
operator-(const T x, const SquareMatrix<1,T>& m) {
  return SquareMatrix<1,T>(x - m[0]);
}

template<typename T>
inline
SquareMatrix<1,T>
operator-(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) {
  return SquareMatrix<1,T>(x[0] - y[0]);
}

template<typename T>
inline
SquareMatrix<1,T>
operator*(const SquareMatrix<1,T>& m, const T x) {
  return SquareMatrix<1,T>(m[0] * x);
}

template<typename T>
inline
SquareMatrix<1,T>
operator*(const T x, const SquareMatrix<1,T>& m) {
  return m * x;
}

template<typename T>
inline
SquareMatrix<1,T>
operator*(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) {
  return SquareMatrix<1,T>(x[0]*y[0]);
}

template<typename T>
inline
SquareMatrix<1,T>
operator/(const SquareMatrix<1,T>& m, const T x) {
#ifdef DEBUG_SquareMatrix
  assert(x != 0);
#endif
  return SquareMatrix<1,T>(m[0] / x);
}

template<typename T>
inline
SquareMatrix<1,T>
operator/(const T x, const SquareMatrix<1,T>& m) {
#ifdef DEBUG_SquareMatrix
  assert(m[0] != 0);
#endif
  return SquareMatrix<1,T>(x / m[0]);
}

//
// Math operators.
//

template<typename T>
inline
T
computeSum(const SquareMatrix<1,T>& x) {
  return (x[0]);
}
    
template<typename T>
inline
T
computeProduct(const SquareMatrix<1,T>& x) {
  return (x[0]);
}

template<typename T>
inline
T
computeMinimum(const SquareMatrix<1,T>& x) {
  return x[0];
}
    
template<typename T>
inline
T
computeMaximum(const SquareMatrix<1,T>& x) {
  return x[0];
}

template<typename T>
inline
T
computeDeterminant(const SquareMatrix<1,T>& x) {
  return x[0];
}

template<typename T>
inline
T
computeTrace(const SquareMatrix<1,T>& x) {
  return x[0];
}

template<typename T>
inline
SquareMatrix<1,T>
computeTranspose(const SquareMatrix<1,T>& x) {
  return x;
}

template<typename T>
inline
SquareMatrix<1,T>
computeInverse(const SquareMatrix<1,T>& x) {
  SquareMatrix<1,T> y;
  computeInverse(x, &y);
  return y;
}

template<typename T>
inline
void
computeInverse(const SquareMatrix<1,T>& x, SquareMatrix<1,T>* y) {
#ifdef DEBUG_SquareMatrix
  assert(x[0] != 0);
#endif
  y[0] = 1 / x[0];
}

template<typename T>
inline
void
computeScaledInverse(const SquareMatrix<1,T>& x, SquareMatrix<1,T>* si) {
  *si = x;
}

template<typename T>
inline
SquareMatrix<1,T>
computeScaledInverse(const SquareMatrix<1,T>& x) {
  return SquareMatrix<1,T>(1);
}

template<typename T>
inline
T
computeFrobeniusNorm(const SquareMatrix<1,T>& x) {
  return std::abs(x[0]);
}

template<typename T>
inline
T
computeFrobeniusNormSquared(const SquareMatrix<1,T>& x) {
  return (x[0]*x[0]);
}

template<typename T>
inline
T
computeInnerProduct(const SquareMatrix<1,T>& x, const SquareMatrix<1,T>& y) {
  return x * y;
}

template<typename T>
inline
void
computeProduct(const SquareMatrix<1,T>& m, const FixedArray<1,T>& v, 
                FixedArray<1,T>* x) {
  (*x)[0] = m[0] * v[0];
}

//
// Equality
//

template<typename T1, typename T2>
inline
bool
operator==(const SquareMatrix<1,T1>& a, const SquareMatrix<1,T2>& b) {
  return (a[0] == b[0]);
}

template<typename T1, typename T2>
inline
bool
operator!=(const SquareMatrix<1,T1>& a, const SquareMatrix<1,T2>& b) { 
  return !(a == b); 
}

//
// I/O
//
    
template<typename T>
std::ostream&
operator<<(std::ostream& out, const SquareMatrix<1,T>& x);
    
template<typename T>
std::istream&
operator>>(std::istream& in, SquareMatrix<1,T>& x);



//----------------------------------------------------------------------------
// 2x2 matrices.
//----------------------------------------------------------------------------


template<typename T>
class SquareMatrix<2,T> :
  public TensorTypes<T> {
private:

  //
  // private types
  //

  typedef TensorTypes<T> base_type;

public:

  //
  // public types
  //

  typedef typename base_type::value_type value_type;
  typedef typename base_type::pointer pointer;
  typedef typename base_type::const_pointer const_pointer;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::index_type index_type;

protected:

  //
  // Data
  //

  value_type _elem[4];

public:

  //
  // Constructors
  //

  SquareMatrix()
  {}

  ~SquareMatrix() 
  {}

  SquareMatrix(const SquareMatrix& x);

  template<typename T2>
  SquareMatrix(const SquareMatrix<2,T2>& x);

  SquareMatrix(const value_type e00, const value_type e01,
               const value_type e10, const value_type e11);

  SquareMatrix(const_pointer x);

  SquareMatrix(const value_type x);

  //
  // Assignment operators
  //

  SquareMatrix& 
  operator=(const SquareMatrix& x);

  SquareMatrix& 
  operator=(const value_type x);

  template<typename T2> 
  SquareMatrix&
  operator=(const SquareMatrix<2,T2>& x);

  //
  // Accessors and Manipulators
  //

  iterator
  begin() { 
    return _elem; 
  }

  iterator 
  end() { 
    return _elem + 4;
  }

  const_iterator 
  begin() const { 
    return _elem;
  }

  const_iterator 
  end() const { 
    return _elem + 4;
  }

  pointer
  data() { 
    return _elem;
  }

  const_pointer
  data() const { 
    return _elem; 
  }

  size_type
  size() const { 
    return 4;
  }

  value_type
  operator()(const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 4);
#endif
    return _elem[i];
  }

  reference
  operator()(const index_type i) { 
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 4);
#endif
    return _elem[i];
  }

  value_type 
  operator[](const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 4);
#endif
    return _elem[i];
  }

  value_type& 
  operator[](const index_type i) {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 4);
#endif
    return _elem[i];
  }

  value_type
  operator()(const index_type i, const index_type j) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 2 && 0 <= j && j < 2);
#endif
    return _elem[i * 2 + j];
  }

  reference
  operator()(const index_type i, const index_type j) {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 2 && 0 <= j && j < 2);
#endif
    return _elem[i * 2 + j];
  }

  void
  get(pointer x) const;

  void
  set(const_pointer x);

  void
  get(reference e00, reference e01,
      reference e10, reference e11) const;

  void
  set(const value_type e00, const value_type e01,
      const value_type e10, const value_type e11);

  void
  negate();

  void
  transpose();

  //
  // Assignment operators with scalar operand.
  //

  SquareMatrix& 
  operator+=(const value_type x);

  SquareMatrix& 
  operator-=(const value_type x);

  SquareMatrix& 
  operator*=(const value_type x);

  SquareMatrix& 
  operator/=(const value_type x);

  SquareMatrix& 
  operator%=(const value_type x);

  //
  // Assignment operators with SquareMatrix operand
  //

  template<typename T2>
  SquareMatrix& 
  operator+=(const SquareMatrix<2,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator-=(const SquareMatrix<2,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator*=(const SquareMatrix<2,T2>& x);

  //
  // Unary operators
  //

  SquareMatrix
  operator+() {
    return *this;
  }
    
  SquareMatrix
  operator-();

};

//
// Binary operators
//
    
template<typename T>
inline
SquareMatrix<2,T>
operator+(const SquareMatrix<2,T>& m, const T x) {
  return SquareMatrix<2,T>(m[0] + x, m[1] + x, 
                           m[2] + x, m[3] + x);
}

template<typename T>
inline
SquareMatrix<2,T>
operator+(const T x, const SquareMatrix<2,T>& m) {
  return m + x;
}

template<typename T>
inline
SquareMatrix<2,T>
operator+(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) {
  return SquareMatrix<2,T>(x[0] + y[0], x[1] + y[1], 
                           x[2] + y[2], x[3] + y[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
operator-(const SquareMatrix<2,T>& m, const T x) {
  return SquareMatrix<2,T>(m[0] - x, m[1] - x, 
                           m[2] - x, m[3] - x);
}

template<typename T>
inline
SquareMatrix<2,T>
operator-(const T x, const SquareMatrix<2,T>& m) {
  return SquareMatrix<2,T>(x - m[0], x - m[1], 
                           x - m[2], x - m[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
operator-(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) {
  return SquareMatrix<2,T>(x[0] - y[0], x[1] - y[1], 
                           x[2] - y[2], x[3] - y[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
operator*(const SquareMatrix<2,T>& m, const T x) {
  return SquareMatrix<2,T>(m[0] * x, m[1] * x, 
                           m[2] * x, m[3] * x);
}

template<typename T>
inline
SquareMatrix<2,T>
operator*(const T x, const SquareMatrix<2,T>& m) {
  return m * x;
}

template<typename T>
inline
SquareMatrix<2,T>
operator*(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) {
  return SquareMatrix<2,T>(x[0]*y[0]+x[1]*y[2],
                           x[0]*y[1]+x[1]*y[3],
                           x[2]*y[0]+x[3]*y[2],
                           x[2]*y[1]+x[3]*y[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
operator/(const SquareMatrix<2,T>& m, const T x) {
#ifdef DEBUG_SquareMatrix
  assert(x != 0);
#endif
  return SquareMatrix<2,T>(m[0] / x, m[1] / x, 
                           m[2] / x, m[3] / x);
}

template<typename T>
inline
SquareMatrix<2,T>
operator/(const T x, const SquareMatrix<2,T>& m) {
#ifdef DEBUG_SquareMatrix
  assert(m[0] != 0 && m[1] != 0 && 
         m[2] != 0 && m[3] != 0);
#endif
  return SquareMatrix<2,T>(x / m[0], x / m[1], 
                           x / m[2], x / m[3]);
}

//
// Math operators.
//

template<typename T>
inline
T
computeSum(const SquareMatrix<2,T>& x) {
  return (x[0] + x[1] + 
          x[2] + x[3]);
}
    
template<typename T>
inline
T
computeProduct(const SquareMatrix<2,T>& x) {
  return (x[0] * x[1] * 
          x[2] * x[3]);
}

template<typename T>
inline
T
computeMinimum(const SquareMatrix<2,T>& x) {
  return *std::min_element(x.begin(), x.end());
}
    
template<typename T>
inline
T
computeMaximum(const SquareMatrix<2,T>& x) {
  return *std::max_element(x.begin(), x.end());
}

template<typename T>
inline
T
computeDeterminant(const SquareMatrix<2,T>& x) {
  return (x[0] * x[3] - x[1] * x[2]);
}

template<typename T>
inline
T
computeTrace(const SquareMatrix<2,T>& x) {
  return (x[0] + x[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
computeTranspose(const SquareMatrix<2,T>& x) {
  return SquareMatrix<2,T>(x[0], x[2],
                           x[1], x[3]);
}

template<typename T>
inline
SquareMatrix<2,T>
computeInverse(const SquareMatrix<2,T>& x) {
  SquareMatrix<2,T> y;
  computeInverse(x, &y);
  return y;
}

template<typename T>
inline
SquareMatrix<2,T>
computeInverse(const SquareMatrix<2,T>& x, const T det) {
  SquareMatrix<2,T> y;
  computeInverse(x, det, &y);
  return y;
}

template<typename T>
inline
void
computeInverse(const SquareMatrix<2,T>& x, SquareMatrix<2,T>* y) {
  const T det = computeDeterminant(x);
  computeInverse(x, det, y);
}

template<typename T>
inline
void
computeInverse(const SquareMatrix<2,T>& x, const T det, SquareMatrix<2,T>* y) {
#ifdef DEBUG_SquareMatrix
  assert(det != 0);
#endif
  y->set(x[3] / det, - x[1] / det,
         - x[2] / det, x[0] / det);
}

template<typename T>
inline
void
computeScaledInverse(const SquareMatrix<2,T>& x, SquareMatrix<2,T>* si) {
  (*si)[0] = x[3];  (*si)[1] = -x[1];
  (*si)[2] = -x[2]; (*si)[3] = x[0];
}

template<typename T>
inline
SquareMatrix<2,T>
computeScaledInverse(const SquareMatrix<2,T>& x) {
  return SquareMatrix<2,T>(x[3], - x[1],
                           - x[2], x[0]);
}

template<typename T>
inline
T
computeFrobeniusNorm(const SquareMatrix<2,T>& x) {
  return std::sqrt(computeFrobeniusNormSquared(x));
}

template<typename T>
inline
T
computeFrobeniusNormSquared(const SquareMatrix<2,T>& x) {
  return (x[0]*x[0] + x[1]*x[1] + 
          x[2]*x[2] + x[3]*x[3]);
}

template<typename T>
inline
T
computeInnerProduct(const SquareMatrix<2,T>& x, const SquareMatrix<2,T>& y) {
  return x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3];
}

template<typename T>
inline
void
computeProduct(const SquareMatrix<2,T>& m, const FixedArray<2,T>& v, 
                FixedArray<2,T>* x) {
  (*x)[0] = m[0] * v[0] + m[1] * v[1];
  (*x)[1] = m[2] * v[0] + m[3] * v[1];
}

//
// Equality
//

template<typename T1, typename T2>
inline
bool
operator==(const SquareMatrix<2,T1>& a, const SquareMatrix<2,T2>& b) {
  return (a[0] == b[0] && a[1] == b[1] && 
          a[2] == b[2] && a[3] == b[3]);
}

template<typename T1, typename T2>
inline
bool
operator!=(const SquareMatrix<2,T1>& a, const SquareMatrix<2,T2>& b) { 
  return !(a == b); 
}

//
// I/O
//
    
template<typename T>
std::ostream&
operator<<(std::ostream& out, const SquareMatrix<2,T>& x);
    
template<typename T>
std::istream&
operator>>(std::istream& in, SquareMatrix<2,T>& x);



//----------------------------------------------------------------------------
// 3x3 matrices.
//----------------------------------------------------------------------------


template<typename T>
class SquareMatrix<3,T> :
  public TensorTypes<T> {
private:

  //
  // private types
  //

  typedef TensorTypes<T> base_type;

public:

  //
  // public types
  //

  typedef typename base_type::value_type value_type;
  typedef typename base_type::pointer pointer;
  typedef typename base_type::const_pointer const_pointer;
  typedef typename base_type::iterator iterator;
  typedef typename base_type::const_iterator const_iterator;
  typedef typename base_type::reference reference;
  typedef typename base_type::const_reference const_reference;
  typedef typename base_type::size_type size_type;
  typedef typename base_type::difference_type difference_type;
  typedef typename base_type::index_type index_type;

protected:

  //
  // Data
  //

  value_type _elem[9];

public:

  //
  // Constructors
  //

  SquareMatrix()
  {}

  ~SquareMatrix() 
  {}

  SquareMatrix(const SquareMatrix& x);

  template<typename T2>
  SquareMatrix(const SquareMatrix<3,T2>& x);

  SquareMatrix(const value_type e00, const value_type e01, const value_type e02, 
               const value_type e10, const value_type e11, const value_type e12, 
               const value_type e20, const value_type e21, const value_type e22);

  SquareMatrix(const_pointer x);

  SquareMatrix(const value_type x);

  //
  // Assignment operators
  //

  SquareMatrix& 
  operator=(const SquareMatrix& x);

  SquareMatrix& 
  operator=(const value_type x);

  template<typename T2> 
  SquareMatrix&
  operator=(const SquareMatrix<3,T2>& x);

  //
  // Accessors and Manipulators
  //

  iterator
  begin() { 
    return _elem; 
  }

  iterator 
  end() { 
    return _elem + 9;
  }

  const_iterator 
  begin() const { 
    return _elem;
  }

  const_iterator 
  end() const { 
    return _elem + 9;
  }

  pointer
  data() { 
    return _elem;
  }

  const_pointer
  data() const { 
    return _elem; 
  }

  size_type
  size() const { 
    return 9;
  }

  value_type
  operator()(const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 9);
#endif
    return _elem[i];
  }

  reference
  operator()(const index_type i) { 
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 9);
#endif
    return _elem[i];
  }

  value_type 
  operator[](const index_type i) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 9);
#endif
    return _elem[i];
  }

  value_type& 
  operator[](const index_type i) {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 9);
#endif
    return _elem[i];
  }

  value_type
  operator()(const index_type i, const index_type j) const {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 3 && 0 <= j && j < 3);
#endif
    return _elem[i * 3 + j];
  }

  reference
  operator()(const index_type i, const index_type j) {
#ifdef DEBUG_SquareMatrix
    assert(0 <= i && i < 3 && 0 <= j && j < 3);
#endif
    return _elem[i * 3 + j];
  }

  void
  get(pointer x) const;

  void
  set(const_pointer x);

  void
  get(reference e00, reference e01, reference e02, 
      reference e10, reference e11, reference e12, 
      reference e20, reference e21, reference e22) const;

  void
  set(const value_type e00, const value_type e01, const value_type e02, 
      const value_type e10, const value_type e11, const value_type e12, 
      const value_type e20, const value_type e21, const value_type e22);

  void
  negate();

  void
  transpose();

  //
  // Assignment operators with scalar operand.
  //

  SquareMatrix& 
  operator+=(const value_type x);

  SquareMatrix& 
  operator-=(const value_type x);

  SquareMatrix& 
  operator*=(const value_type x);

  SquareMatrix& 
  operator/=(const value_type x);

  SquareMatrix& 
  operator%=(const value_type x);

  //
  // Assignment operators with SquareMatrix operand
  //

  template<typename T2>
  SquareMatrix& 
  operator+=(const SquareMatrix<3,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator-=(const SquareMatrix<3,T2>& x);

  template<typename T2>
  SquareMatrix& 
  operator*=(const SquareMatrix<3,T2>& x);

  //
  // Unary operators
  //

  SquareMatrix
  operator+() {
    return *this;
  }
    
  SquareMatrix
  operator-();

};

//
// Binary operators
//
    
template<typename T>
inline
SquareMatrix<3,T>
operator+(const SquareMatrix<3,T>& m, const T x) {
  return SquareMatrix<3,T>(m[0] + x, m[1] + x, m[2] + x, 
                           m[3] + x, m[4] + x, m[5] + x, 
                           m[6] + x, m[7] + x, m[8] + x);
}

template<typename T>
inline
SquareMatrix<3,T>
operator+(const T x, const SquareMatrix<3,T>& m) {
  return m + x;
}

template<typename T>
inline
SquareMatrix<3,T>
operator+(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) {
  return SquareMatrix<3,T>(x[0] + y[0], x[1] + y[1], x[2] + y[2],
                           x[3] + y[3], x[4] + y[4], x[5] + y[5],
                           x[6] + y[6], x[7] + y[7], x[8] + y[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
operator-(const SquareMatrix<3,T>& m, const T x) {
  return SquareMatrix<3,T>(m[0] - x, m[1] - x, m[2] - x, 
                           m[3] - x, m[4] - x, m[5] - x, 
                           m[6] - x, m[7] - x, m[8] - x);
}

template<typename T>
inline
SquareMatrix<3,T>
operator-(const T x, const SquareMatrix<3,T>& m) {
  return SquareMatrix<3,T>(x - m[0], x - m[1], x - m[2], 
                           x - m[3], x - m[4], x - m[5], 
                           x - m[6], x - m[7], x - m[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
operator-(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) {
  return SquareMatrix<3,T>(x[0] - y[0], x[1] - y[1], x[2] - y[2],
                           x[3] - y[3], x[4] - y[4], x[5] - y[5],
                           x[6] - y[6], x[7] - y[7], x[8] - y[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
operator*(const SquareMatrix<3,T>& m, const T x) {
  return SquareMatrix<3,T>(m[0] * x, m[1] * x, m[2] * x, 
                           m[3] * x, m[4] * x, m[5] * x, 
                           m[6] * x, m[7] * x, m[8] * x);
}

template<typename T>
inline
SquareMatrix<3,T>
operator*(const T x, const SquareMatrix<3,T>& m) {
  return m * x;
}

template<typename T>
inline
SquareMatrix<3,T>
operator*(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) {
  return SquareMatrix<3,T>
    (x[0]*y[0]+x[1]*y[3]+x[2]*y[6],
     x[0]*y[1]+x[1]*y[4]+x[2]*y[7],
     x[0]*y[2]+x[1]*y[5]+x[2]*y[8],
     x[3]*y[0]+x[4]*y[3]+x[5]*y[6],
     x[3]*y[1]+x[4]*y[4]+x[5]*y[7],
     x[3]*y[2]+x[4]*y[5]+x[5]*y[8],
     x[6]*y[0]+x[7]*y[3]+x[8]*y[6],
     x[6]*y[1]+x[7]*y[4]+x[8]*y[7],
     x[6]*y[2]+x[7]*y[5]+x[8]*y[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
operator/(const SquareMatrix<3,T>& m, const T x) {
#ifdef DEBUG_SquareMatrix
  assert(x != 0);
#endif
  return SquareMatrix<3,T>(m[0] / x, m[1] / x, m[2] / x, 
                           m[3] / x, m[4] / x, m[5] / x, 
                           m[6] / x, m[7] / x, m[8] / x);
}

template<typename T>
inline
SquareMatrix<3,T>
operator/(const T x, const SquareMatrix<3,T>& m) {
#ifdef DEBUG_SquareMatrix
  assert(m[0] != 0 && m[1] != 0 && m[2] != 0 &&
         m[3] != 0 && m[4] != 0 && m[5] != 0 &&
         m[6] != 0 && m[7] != 0 && m[8] != 0);
#endif
  return SquareMatrix<3,T>(x / m[0], x / m[1], x / m[2], 
                           x / m[3], x / m[4], x / m[5], 
                           x / m[6], x / m[7], x / m[8]);
}

//
// Math operators.
//

template<typename T>
inline
T
computeSum(const SquareMatrix<3,T>& x) {
  return (x[0] + x[1] + x[2] + 
          x[3] + x[4] + x[5] + 
          x[6] + x[7] + x[8]);
}
    
template<typename T>
inline
T
computeProduct(const SquareMatrix<3,T>& x) {
  return (x[0] * x[1] * x[2] * 
          x[3] * x[4] * x[5] * 
          x[6] * x[7] * x[8]);
}

template<typename T>
inline
T
computeMinimum(const SquareMatrix<3,T>& x) {
  return *std::min_element(x.begin(), x.end());
}
    
template<typename T>
inline
T
computeMaximum(const SquareMatrix<3,T>& x) {
  return *std::max_element(x.begin(), x.end());
}

template<typename T>
inline
T
computeDeterminant(const SquareMatrix<3,T>& x) {
  return (x[0]*x[4]*x[8] + x[1]*x[5]*x[6] + x[2]*x[3]*x[7] -
          x[2]*x[4]*x[6] - x[1]*x[3]*x[8] - x[0]*x[5]*x[7]);
}

template<typename T>
inline
T
computeTrace(const SquareMatrix<3,T>& x) {
  return (x[0] + x[4] + x[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
computeTranspose(const SquareMatrix<3,T>& x) {
  return SquareMatrix<3,T>(x[0], x[3], x[6],
                           x[1], x[4], x[7],
                           x[2], x[5], x[8]);
}

template<typename T>
inline
SquareMatrix<3,T>
computeInverse(const SquareMatrix<3,T>& x) {
  SquareMatrix<3,T> y;
  computeInverse(x, &y);
  return y;
}

template<typename T>
inline
SquareMatrix<3,T>
computeInverse(const SquareMatrix<3,T>& x, const T det) {
  SquareMatrix<3,T> y;
  computeInverse(x, det, &y);
  return y;
}

template<typename T>
inline
void
computeInverse(const SquareMatrix<3,T>& x, SquareMatrix<3,T>* y) {
  const T det = computeDeterminant(x);
  computeInverse(x, det, y);
}

template<typename T>
inline
void
computeInverse(const SquareMatrix<3,T>& x, const T det, SquareMatrix<3,T>* y) {
#ifdef DEBUG_SquareMatrix
  assert(det != 0);
#endif
  y->set((x[4] * x[8] - x[5] * x[7]) / det,
        (x[2] * x[7] - x[1] * x[8]) / det,
        (x[1] * x[5] - x[2] * x[4]) / det,
        (x[5] * x[6] - x[3] * x[8]) / det,
        (x[0] * x[8] - x[2] * x[6]) / det,
        (x[2] * x[3] - x[0] * x[5]) / det,
        (x[3] * x[7] - x[4] * x[6]) / det,
        (x[1] * x[6] - x[0] * x[7]) / det,
        (x[0] * x[4] - x[1] * x[3]) / det);
}

template<typename T>
inline
void
computeScaledInverse(const SquareMatrix<3,T>& x, SquareMatrix<3,T>* si) {
  (*si)[0] = x[4] * x[8] - x[5] * x[7];
  (*si)[1] = x[2] * x[7] - x[1] * x[8];
  (*si)[2] = x[1] * x[5] - x[2] * x[4];
  (*si)[3] = x[5] * x[6] - x[3] * x[8];
  (*si)[4] = x[0] * x[8] - x[2] * x[6];
  (*si)[5] = x[2] * x[3] - x[0] * x[5];
  (*si)[6] = x[3] * x[7] - x[4] * x[6];
  (*si)[7] = x[1] * x[6] - x[0] * x[7];
  (*si)[8] = x[0] * x[4] - x[1] * x[3];
}

template<typename T>
inline
SquareMatrix<3,T>
computeScaledInverse(const SquareMatrix<3,T>& x) {
  return SquareMatrix<3,T>(x[4] * x[8] - x[5] * x[7],
                           x[2] * x[7] - x[1] * x[8],
                           x[1] * x[5] - x[2] * x[4],
                           x[5] * x[6] - x[3] * x[8],
                           x[0] * x[8] - x[2] * x[6],
                           x[2] * x[3] - x[0] * x[5],
                           x[3] * x[7] - x[4] * x[6],
                           x[1] * x[6] - x[0] * x[7],
                           x[0] * x[4] - x[1] * x[3]);
}

template<typename T>
inline
T
computeFrobeniusNorm(const SquareMatrix<3,T>& x) {
  return std::sqrt(computeFrobeniusNormSquared(x));
}

template<typename T>
inline
T
computeFrobeniusNormSquared(const SquareMatrix<3,T>& x) {
  return (x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + 
          x[3]*x[3] + x[4]*x[4] + x[5]*x[5] + 
          x[6]*x[6] + x[7]*x[7] + x[8]*x[8]);
}

template<typename T>
inline
T
computeInnerProduct(const SquareMatrix<3,T>& x, const SquareMatrix<3,T>& y) {
  return x[0] * y[0] + x[1] * y[1] + x[2] * y[2] + x[3] * y[3] + x[4] * y[4]
    + x[5] * y[5] + x[6] * y[6] + x[7] * y[7] + x[8] * y[8];
}

template<typename T>
inline
void
computeProduct(const SquareMatrix<3,T>& m, const FixedArray<3,T>& v, 
                FixedArray<3,T>* x) {
  (*x)[0] = m[0] * v[0] + m[1] * v[1] + m[2] * v[2];
  (*x)[1] = m[3] * v[0] + m[4] * v[1] + m[5] * v[2];
  (*x)[2] = m[6] * v[0] + m[7] * v[1] + m[8] * v[2];
}

//
// Equality
//

template<typename T1, typename T2>
inline
bool
operator==(const SquareMatrix<3,T1>& a, const SquareMatrix<3,T2>& b) {
  return (a[0] == b[0] && a[1] == b[1] && a[2] == b[2] &&
          a[3] == b[3] && a[4] == b[4] && a[5] == b[5] &&
          a[6] == b[6] && a[7] == b[7] && a[8] == b[8]);
}

template<typename T1, typename T2>
inline
bool
operator!=(const SquareMatrix<3,T1>& a, const SquareMatrix<3,T2>& b) { 
  return !(a == b); 
}

//
// I/O
//
    
template<typename T>
std::ostream&
operator<<(std::ostream& out, const SquareMatrix<3,T>& x);
    
template<typename T>
std::istream&
operator>>(std::istream& in, SquareMatrix<3,T>& x);

END_NAMESPACE_ADS

#define __SquareMatrix_ipp__
#include "SquareMatrix.ipp"
#undef __SquareMatrix_ipp__

#endif

---