TMatrix.h

Go to the documentation of this file.

// SPDX-License-Identifier: LGPL-2.1-or-later OR GPL-2.0-or-later OR GPL-3.0-or-later OR LicenseRef-ACF-Commercial
#pragma once
 
 
// ACF includes
#include <istd/TArray.h>
#include <istd/CIndex2d.h>
#include <iser/CArchiveTag.h>
#include <imath/TVector.h>
 
 
namespace iser
{
       class IArchive;
}
 
 
namespace imath
{
 
 
template <int Width, int Height, typename Element = double>
class TMatrix
{
public:
       typedef istd::CIndex2d IndexType;
       typedef istd::CIndex2d SizesType;
       typedef Element ElementType;
       typedef imath::TVector<Height, Element> ColumnVector;
       typedef imath::TVector<Width, Element> RowVector;
 
       enum MatrixInitMode
       {
              MIM_ZERO,
              MIM_ONES,
              MIM_IDENTITY,
              MIM_UPPER_TRIANGLE,
              MIM_DOWNER_TRIANGLE
       };
 
       TMatrix();
 
       explicit TMatrix(MatrixInitMode mode);
 
       TMatrix(const TMatrix& matrix);
 
       void Reset();
 
       void Clear();
 
       bool IsDimensionsCountFixed() const;
 
       int GetDimensionsCount() const;
 
       bool SetDimensionsCount(int count);
 
       const SizesType& GetSizes() const;
 
       bool SetSizes(const SizesType& sizes);
 
       int GetSize(int dimension) const;
 
       bool SetSize(int dimension, int size);
 
       const ElementType& GetAt(const IndexType& index) const;
 
       const ElementType& GetAt(int x, int y) const;
 
       ElementType& GetAtRef(const IndexType& index);
 
       ElementType& GetAtRef(int x, int y);
 
       void SetAt(const IndexType& index, const ElementType& value);
 
       void SetAt(int x, int y, const ElementType& value);
 
       void InitToIdentity();
 
       /*
              Get maximal element.
       */
       double GetMaxElement() const;
 
       /*
              Get minimal element.
       */
       double GetMinElement() const;
 
       void GetNegated(TMatrix<Width, Height, Element>& result);
 
       void GetAdded(const TMatrix<Width, Height, Element>& matrix, TMatrix<Width, Height, Element>& result) const;
 
       void GetSubstracted(const TMatrix<Width, Height, Element>& matrix, TMatrix<Width, Height, Element>& result) const;
 
       template <int SecondWidth>
       void GetMultiplied(const TMatrix<SecondWidth, Width, Element>& matrix, TMatrix<SecondWidth, Height, Element>& result) const
       {
              for (int resultY = 0; resultY < Height; ++resultY){
                     for (int resultX = 0; resultX < SecondWidth; ++resultX){
                            double sum = 0;
                            for (int i = 0; i < Width; ++i){
                                   sum += m_elements[i][resultY] * matrix.m_elements[resultX][i];
                            }
 
                            result.m_elements[resultX][resultY] = sum;
                     }
              }
       }
 
       template <int SecondWidth>
       TMatrix<SecondWidth, Height, Element> GetMultiplied(const TMatrix<SecondWidth, Width, Element>& matrix) const
       {
              TMatrix<SecondWidth, Height, Element> retVal;
 
              GetMultiplied(matrix, retVal);
 
              return retVal;
       }
 
       void GetMultiplied(const TVector<Width, Element>& vector, TVector<Height, Element>& result) const;
 
       TVector<Height, Element> GetMultiplied(const TVector<Width, Element>& vector) const;
 
       void GetScaled(double value, TMatrix<Width, Height, Element>& result) const;
 
       void GetTransposed(TMatrix<Height, Width, Element>& result) const;
 
       TMatrix<Height, Width, Element> GetTransposed() const;
 
       void Transpose();
 
       double GetTrace() const;
 
       /*
              Get square of euclidean norm.
       */
       double GetFrobeniusNorm2() const;
 
       /*
              Get euclidean norm.
       */
       double GetFrobeniusNorm() const;
 
       bool GetTriangleDecomposed(
                            TMatrix<Width, Height, Element>& result,
                            TMatrix<Height, Height, Element>* matrixQPtr = NULL,
                            int maxColumns = -1,
                            double minHhNorm = I_BIG_EPSILON) const;
 
       bool GetDecompositionQDQ(TMatrix<Height, Height, Element>& matrixQ, TVector<Height, Element>& diagonalD, double tolerance = I_BIG_EPSILON, int maxIterations = 100) const;
 
       void GetColumnVector(int columnIndex, TVector<Height, Element>& result) const;
 
       void SetColumnVector(int columnIndex, const TVector<Height, Element>& columnVector);
 
       void GetRowVector(int rowIndex, TVector<Width, Element>& result);
 
       bool Serialize(iser::IArchive& archive);
 
       // operators
       TMatrix<Width, Height, Element> operator+(const TMatrix<Width, Height, Element>& matrix) const;
       TMatrix<Width, Height, Element> operator-(const TMatrix<Width, Height, Element>& matrix) const;
       TMatrix<Width, Height, Element> operator-();
 
       template <int SecondWidth>
       TMatrix<SecondWidth, Height, Element> operator*(const TMatrix<SecondWidth, Width, Element>& matrix) const
       {
              TMatrix<SecondWidth, Height, Element> retVal;
 
              GetMultiplied(matrix, retVal);
 
              return retVal;
       }
 
       TMatrix<Width, Height, Element> operator*(double value) const;
 
       TMatrix<Width, Height, Element>& operator+=(const TMatrix<Width, Height, Element>& matrix);
       TMatrix<Width, Height, Element>& operator-=(const TMatrix<Width, Height, Element>& matrix);
       TMatrix<Width, Height, Element>& operator*=(double value);
 
       bool operator==(const TMatrix<Width, Height, Element>& matrix) const;
       bool operator!=(const TMatrix<Width, Height, Element>& matrix) const;
 
       const ElementType& operator[](const IndexType& index) const;
       ElementType& operator[](const IndexType& index);
 
    TMatrix<Width, Height, Element>& operator=(const TMatrix<Width, Height, Element>& matrix) = default;
 
private:
       typedef Element Elements[Width][Height];
 
       Elements m_elements;
};
 
 
// inline methods
 
template <int Width, int Height, typename Element>
inline TMatrix<Height, Width, Element> TMatrix<Width, Height, Element>::GetTransposed() const
{
       TMatrix<Height, Width, Element> retVal;
 
       GetTransposed(retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
inline void TMatrix<Width, Height, Element>::Transpose()
{
       TMatrix<Height, Width, Element> retVal;
 
       GetTransposed(retVal);
 
       *this = retVal;
}
 
 
template <int Width, int Height, typename Element>
TVector<Height, Element> TMatrix<Width, Height, Element>::GetMultiplied(const TVector<Width, Element>& vector) const
{
       TVector<Height, Element> retVal;
 
       GetMultiplied(vector, retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
inline TMatrix<Width, Height, Element> TMatrix<Width, Height, Element>::operator+(const TMatrix<Width, Height, Element>& matrix) const
{
       TMatrix<Width, Height, Element> retVal;
 
       GetAdded(matrix, retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
inline TMatrix<Width, Height, Element> TMatrix<Width, Height, Element>::operator-(const TMatrix<Width, Height, Element>& matrix) const
{
       TMatrix<Width, Height, Element> retVal;
 
       GetSubstracted(matrix, retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
inline TMatrix<Width, Height, Element> TMatrix<Width, Height, Element>::operator-()
{
       TMatrix<Width, Height, Element> retVal;
 
       GetNegated(retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
inline TMatrix<Width, Height, Element> TMatrix<Width, Height, Element>::operator*(double value) const
{
       TMatrix retVal;
 
       GetScaled(value, retVal);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>& TMatrix<Width, Height, Element>::operator+=(const TMatrix<Width, Height, Element>& matrix)
{
       GetAdded(matrix, *this);
 
       return *this;
}
 
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>& TMatrix<Width, Height, Element>::operator-=(const TMatrix<Width, Height, Element>& matrix)
{
       GetSubstracted(matrix, *this);
 
       return *this;
}
 
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>& TMatrix<Width, Height, Element>::operator*=(double value)
{
       GetScaled(value, *this);
 
       return *this;
}
 
 
template <int Width, int Height, typename Element>
inline TMatrix<Width, Height, Element> operator*(double value, const imath::TMatrix<Width, Height, Element>& matrix)
{
       return matrix * value;
}
 
 
// public methods
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>::TMatrix()
{
}
 
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>::TMatrix(MatrixInitMode mode)
{
       switch (mode){
       case MIM_ZERO:
              for (int x = 0; x < Width; ++x){
                     for (int y = 0; y < Height; ++y){
                            m_elements[x][y] = 0;
                     }
              }
              break;
 
       case MIM_ONES:
              for (int x = 0; x < Width; ++x){
                     for (int y = 0; y < Height; ++y){
                            m_elements[x][y] = 1;
                     }
              }
              break;
 
       case MIM_IDENTITY:
              for (int x = 0; x < Width; ++x){
                     for (int y = 0; y < Height; ++y){
                            if (x == y){
                                   m_elements[x][y] = 1;
                            }
                            else{
                                   m_elements[x][y] = 0;
                            }
                     }
              }
              break;
 
       case MIM_UPPER_TRIANGLE:
              for (int x = 0; x < Width; ++x){
                     for (int y = 0; y < Height; ++y){
                            if (x <= y){
                                   m_elements[x][y] = 1;
                            }
                            else{
                                   m_elements[x][y] = 0;
                            }
                     }
              }
              break;
 
       case MIM_DOWNER_TRIANGLE:
              for (int x = 0; x < Width; ++x){
                     for (int y = 0; y < Height; ++y){
                            if (x >= y){
                                   m_elements[x][y] = 1;
                            }
                            else{
                                   m_elements[x][y] = 0;
                            }
                     }
              }
              break;
       }
}
 
 
template <int Width, int Height, typename Element>
TMatrix<Width, Height, Element>::TMatrix(const TMatrix& matrix)
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     m_elements[x][y] = matrix.m_elements[x][y];
              }
       }
}
 
 
// init operations
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::Reset()
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     m_elements[x][y] = 0;
              }
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::Clear()
{
       Reset();
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::IsDimensionsCountFixed() const
{
       return true;
}
 
 
template <int Width, int Height, typename Element>
int TMatrix<Width, Height, Element>::GetDimensionsCount() const
{
       return 2;
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::SetDimensionsCount(int count)
{
       return (count == 2);
}
 
 
template <int Width, int Height, typename Element>
const typename TMatrix<Width, Height, Element>::SizesType& TMatrix<Width, Height, Element>::GetSizes() const
{
       static SizesType retVal(Width, Height);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::SetSizes(const SizesType& sizes)
{
       return (sizes[0] == Width) && (sizes[1] == Height);
}
 
 
template <int Width, int Height, typename Element>
int TMatrix<Width, Height, Element>::GetSize(int dimension) const
{
       switch (dimension){
       case 0:
              return Width;
 
       case 1:
              return Height;
 
       default:
              return 0;
       }
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::SetSize(int dimension, int size)
{
       switch (dimension){
       case 0:
              return (size == Width);
 
       case 1:
              return (size == Height);
 
       default:
              return false;
       }
}
 
 
template <int Width, int Height, typename Element>
const typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::GetAt(const IndexType& index) const
{
       return m_elements[index[0]][index[1]];
}
 
 
template <int Width, int Height, typename Element>
const typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::GetAt(int x, int y) const
{
       return m_elements[x][y];
}
 
 
template <int Width, int Height, typename Element>
typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::GetAtRef(const IndexType& index)
{
       return m_elements[index[0]][index[1]];
}
 
 
template <int Width, int Height, typename Element>
typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::GetAtRef(int x, int y)
{
       return m_elements[x][y];
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::SetAt(const IndexType& index, const ElementType& value)
{
       m_elements[index[0]][index[1]] = value;
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::SetAt(int x, int y, const ElementType& value)
{
       m_elements[x][y] = value;
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::InitToIdentity()
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     if (x != y){
                            m_elements[x][y] = 0;
                     }
                     else{
                            m_elements[x][y] = 1;
                     }
              }
       }
}
 
 
template <int Width, int Height, typename Element>
double TMatrix<Width, Height, Element>::GetMaxElement() const
{
       double retVal = m_elements[0][0];
 
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     double value = m_elements[x][y];
 
                     if (value > retVal){
                            retVal = value;
                     }
              }
       }
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
double TMatrix<Width, Height, Element>::GetMinElement() const
{
       double retVal = m_elements[0][0];
 
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     double value = m_elements[x][y];
 
                     if (value < retVal){
                            retVal = value;
                     }
              }
       }
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetNegated(TMatrix<Width, Height, Element>& result)
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     result.m_elements[x][y] = -m_elements[x][y];
              }
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetAdded(const TMatrix<Width, Height, Element>& matrix, TMatrix<Width, Height, Element>& result) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     result.m_elements[x][y] = m_elements[x][y] + matrix.m_elements[x][y];
              }
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetSubstracted(const TMatrix<Width, Height, Element>& matrix, TMatrix<Width, Height, Element>& result) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     result.m_elements[x][y] = m_elements[x][y] - matrix.m_elements[x][y];
              }
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetMultiplied(const TVector<Width, Element>& vector, TVector<Height, Element>& result) const
{
       for (int resultY = 0; resultY < Height; ++resultY){
              double sum = 0;
              for (int i = 0; i < Width; ++i){
                     sum += m_elements[i][resultY] * vector[i];
              }
 
              result[resultY] = sum;
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetScaled(double value, TMatrix<Width, Height, Element>& result) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     result.m_elements[x][y] = m_elements[x][y] * value;
              }
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetTransposed(TMatrix<Height, Width, Element>& result) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     result.SetAt(y, x, GetAt(x, y));
              }
       }
}
 
 
template <int Width, int Height, typename Element>
double TMatrix<Width, Height, Element>::GetTrace() const
{
       int commonSize = qMin(Width, Height);
 
       double retVal = 0;
       for (int i = 0; i < commonSize; ++i){
              retVal += m_elements[i][i];
       }
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
double TMatrix<Width, Height, Element>::GetFrobeniusNorm2() const
{
       double retVal = 0.0;
 
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     double value = m_elements[x][y];
                     retVal += double(value * value);
              }
       }
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
double TMatrix<Width, Height, Element>::GetFrobeniusNorm() const
{
       return qSqrt(GetFrobeniusNorm2());
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::GetDecompositionQDQ(TMatrix<Height, Height, Element>& matrixQ, TVector<Height, Element>& diagonalD, double tolerance, int maxIterations) const
{
       if (Width != Height){
              return false;
       }
 
       TMatrix<Width, Height, Element> matrixR = *this;
       matrixQ.Clear();
 
       for (int i = 0; i < maxIterations; ++i){
              TMatrix<Width, Height, Element> tempMatrixR;
              TMatrix<Height, Height, Element> tempMatrixQ;
              if (!matrixR.GetTriangleDecomposed(tempMatrixR, &tempMatrixQ)){
                     return false;
              }
              matrixR = tempMatrixR.GetMultiplied(tempMatrixQ);
 
              matrixQ = tempMatrixQ.GetMultiplied(matrixQ);
 
              double residue = 0;
              istd::CIndex2d index;
              for (index[0] = 0; index[0] < Width; ++index[0]){
                     for (index[1] = 0; index[1] < Width; ++index[1]){
                            if (index[0] != index[1]){
                                   double element = matrixR.GetAt(index);
 
                                   residue += element * element;
                            }
                     }
              }
 
              if (residue <= tolerance){
                     for (int i = 0; i < Height; ++i){
                            diagonalD[i] = matrixR.GetAt(i, i);
                     }
 
                     return true;
              }
       }
 
       return false;
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetColumnVector(int columnIndex, TVector<Height, Element>& result) const
{
       Q_ASSERT(columnIndex < Width);
       Q_ASSERT(columnIndex >= 0);
 
       for (int y = 0; y < Height; ++y){
              result[y] = m_elements[columnIndex][y];
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::SetColumnVector(int columnIndex, const TVector<Height, Element>& columnVector)
{
       Q_ASSERT(columnIndex < Width);
       Q_ASSERT(columnIndex >= 0);
 
       for (int y = 0; y < Height; ++y){
              m_elements[columnIndex][y] = columnVector[y];
       }
}
 
 
template <int Width, int Height, typename Element>
void TMatrix<Width, Height, Element>::GetRowVector(int rowIndex, TVector<Width, Element>& result)
{
       Q_ASSERT(rowIndex < Height);
 
       for (int x = 0; x < Width; ++x){
              result[x] = m_elements[x][rowIndex];
       }
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::GetTriangleDecomposed(
                     TMatrix<Width, Height, Element>& result,
                     TMatrix<Height, Height, Element>* matrixQPtr,
                     int maxColumns,
                     double minHhNorm) const
{
       if (&result != this){
              result = *this;
       }
 
       int columnsCount = qMin(Height - 1, Width);
       if ((maxColumns >= 0) && (maxColumns < columnsCount)){
              columnsCount = maxColumns;
       }
 
       for (int hhIndex = 0; hhIndex < columnsCount; ++hhIndex){
              double hhNorm2 = 0;
              for (istd::CIndex2d normIndex(hhIndex, hhIndex); normIndex[1] < Height; ++normIndex[1]){
                     double element = result.GetAt(normIndex);
 
                     hhNorm2 += element * element;
              }
 
              if (hhNorm2 < minHhNorm){
                     return false;
              }
 
              double element0 = result[istd::CIndex2d(hhIndex, hhIndex)];
 
              double hhLength = (element0 >= 0)? qSqrt(hhNorm2): -qSqrt(hhNorm2);   // destination diagonal value of current processed column, sign is choosen optimal for maximal householder vector length
 
              double hhVector0 = element0 + hhLength;   // element 0 of householder vector, rest of ist will be taken directly from matrix elements
              double hhVectorNorm2 = hhNorm2 + hhVector0 * hhVector0 - element0 * element0;       // sqare of norm of householder vector, first element replaced from original column norm sqare
 
              // correct next columns
              istd::CIndex2d matrixIndex;
              for (matrixIndex[0] = hhIndex + 1; matrixIndex[0] < Width; ++matrixIndex[0]){
                     matrixIndex[1] = hhIndex;
                     double dotProduct = result[matrixIndex] * hhVector0;
                     for (matrixIndex[1]++; matrixIndex[1] < Height; ++matrixIndex[1]){
                            dotProduct += result.GetAt(matrixIndex) * result.GetAt(istd::CIndex2d(hhIndex, matrixIndex[1]));
                     }
 
                     double reflexionFactor = 2 * dotProduct / hhVectorNorm2;
                     matrixIndex[1] = hhIndex;
                     result[matrixIndex] -= hhVector0 * reflexionFactor;
                     for (matrixIndex[1]++; matrixIndex[1] < Height; ++matrixIndex[1]){
                            result.GetAtRef(matrixIndex) -= result[istd::CIndex2d(hhIndex, matrixIndex[1])] * reflexionFactor;
                     }
              }
 
              // correct optional matrix
              if (matrixQPtr != NULL){
                     istd::CIndex2d matrix2Index;
                     for (matrix2Index[0] = 0; matrix2Index[0] < Height; ++matrix2Index[0]){
                            matrix2Index[1] = hhIndex;
                            double dotProduct = matrixQPtr->GetAt(matrix2Index) * hhVector0;
                            for (++matrix2Index[1]; matrix2Index[1] < Height; ++matrix2Index[1]){
                                   dotProduct += matrixQPtr->GetAt(matrix2Index) * result[istd::CIndex2d(hhIndex, matrix2Index[1])];
                            }
 
                            double reflexionFactor = 2 * dotProduct / hhVectorNorm2;
                            matrix2Index[1] = hhIndex;
                            matrixQPtr->GetAtRef(matrix2Index) -= hhVector0 * reflexionFactor;
                            for (matrix2Index[1]++; matrix2Index[1] < Height; ++matrix2Index[1]){
                                   matrixQPtr->GetAtRef(matrix2Index) -= result.GetAt(istd::CIndex2d(hhIndex, matrix2Index[1])) * reflexionFactor;
                            }
                     }
              }
 
              // correct current column
              result.SetAt(istd::CIndex2d(hhIndex, hhIndex), -hhLength);
              for (int i = hhIndex + 1; i < Height; ++i){
                     result.SetAt(istd::CIndex2d(hhIndex, i), 0);
              }
       }
 
       return true;
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::Serialize(iser::IArchive& archive)
{
       bool retVal = true;
 
       static iser::CArchiveTag sizeTag("Size", "Size of matrix", iser::CArchiveTag::TT_GROUP);
       static iser::CArchiveTag cellsTag("Cells", "Values of matrix cells", iser::CArchiveTag::TT_GROUP);
 
       retVal = retVal && archive.BeginTag(sizeTag);
       if (archive.IsStoring()){
              int width = Width;
              int height = Height;
 
              retVal = retVal && archive.Process(width);
              retVal = retVal && archive.Process(height);
       }
       else{
              int width = 0;
              int height = 0;
 
              retVal = retVal && archive.Process(width);
              retVal = retVal && archive.Process(height);
 
              if (!retVal || (width != Width) || (height != Height)){
                     return false;
              }
       }
 
       retVal = retVal && archive.EndTag(sizeTag);
 
       retVal = retVal && archive.BeginTag(cellsTag);
       for (int y = 0; y < Height; ++y){
              for (int x = 0; x < Width; ++x){
                     retVal = retVal && archive.Process(m_elements[x][y]);
              }
       }
       retVal = retVal && archive.EndTag(cellsTag);
 
       return retVal;
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::operator==(const TMatrix<Width, Height, Element>& matrix) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     if (m_elements[x][y] != matrix.m_elements[x][y]){
                            return false;
                     }
              }
       }
 
       return true;
}
 
 
template <int Width, int Height, typename Element>
bool TMatrix<Width, Height, Element>::operator!=(const TMatrix<Width, Height, Element>& matrix) const
{
       for (int x = 0; x < Width; ++x){
              for (int y = 0; y < Height; ++y){
                     if (m_elements[x][y] != matrix.m_elements[x][y]){
                            return true;
                     }
              }
       }
 
       return false;
}
 
 
template <int Width, int Height, typename Element>
const typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::operator[](const IndexType& index) const
{
       return m_elements[index[0]][index[1]];
}
 
 
template <int Width, int Height, typename Element>
typename TMatrix<Width, Height, Element>::ElementType& TMatrix<Width, Height, Element>::operator[](const IndexType& index)
{
       return m_elements[index[0]][index[1]];
}
 
 
} // namespace imath