MatrixBase.hpp

Go to the documentation of this file.

//==============================================================================
//
//  This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
//
//  The GNSSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 3.0 of the License, or
//  any later version.
//
//  The GNSSTk 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 Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GNSSTk; if not, write to the Free Software Foundation,
//  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin.
//  Copyright 2004-2022, The Board of Regents of The University of Texas System
//
//==============================================================================
 
//==============================================================================
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin, under contract to an agency or agencies
//  within the U.S. Department of Defense. The U.S. Government retains all
//  rights to use, duplicate, distribute, disclose, or release this software.
//
//  Pursuant to DoD Directive 523024
//
//  DISTRIBUTION STATEMENT A: This software has been approved for public
//                            release, distribution is unlimited.
//
//==============================================================================
 
#ifndef GNSSTK_MATRIX_BASE_HPP
#define GNSSTK_MATRIX_BASE_HPP
 
#include "Vector.hpp"
 
namespace gnsstk
{
 
 
   NEW_EXCEPTION_CLASS(MatrixException, Exception);
 
   NEW_EXCEPTION_CLASS(SingularMatrixException, MatrixException);
 
   template <class T, class BaseClass>
   class ConstMatrixBase
   {
   public:
      explicit ConstMatrixBase() {}
 
      virtual ~ConstMatrixBase()
      {}
 
      size_t size() const
      { return static_cast<const BaseClass*>(this)->size(); }
      size_t cols() const
      { return static_cast<const BaseClass*>(this)->cols(); }
      size_t rows() const
      { return static_cast<const BaseClass*>(this)->rows(); }
      T operator() (size_t i, size_t j) const
      { return constMatrixRef(i, j); }
 
      inline bool isSquare() const
      { return ((rows() == cols()) && (rows() != 0)); }
      inline bool isUT() const
      {
         if (!isSquare())
            return false;
         size_t i, j;
         for (i = 1; i < rows(); i++)
            for (j = 0; j < i; j++)
               if ((*this)(i,j) != T(0))
                  return false;
         return true;
      }
      inline bool isLT() const
      {
         if (!isSquare())
            return false;
         size_t i, j;
         for (i = 0; i < rows(); i++)
            for (j = i+1; j < cols(); j++)
               if ((*this)(i,j) != T(0))
                  return false;
         return true;
      }
 
      inline bool isDiagonal() const
      {
         if (!isSquare())
            return false;
         size_t i, j;
         for (i = 0; i < rows(); i++)
            for (j = 0; j < cols(); j++)
               if (i != j)
                  if ((*this)(i,j) != T(0))
                     return false;
         return true;
      }
      inline bool isSymmetric() const
      {
         if (!isSquare())
            return false;
         size_t i,j;
         for (i = 0; i < rows(); i++)
            for (j = i + 1; j < cols(); j++)
               if ((*this)(i,j) != (*this)(j,i))
                  return false;
         return true;
      }
 
      Vector<T> colCopy(size_t c, size_t r = 0) const
      {
#ifdef RANGECHECK
         if ((c >= cols()) || (r >= rows()))
         {
            MatrixException e("Invalid ConstMatrixBase index for colCopy");
            GNSSTK_THROW(e);
         }
#endif
         Vector<T> toReturn(rows() - r);
         size_t i;
         for (i = r; i < rows(); i++)
            toReturn(i - r) = (*this)(i, c);
         return toReturn;
      }
 
      Vector<T> rowCopy(size_t r, size_t c = 0) const
      {
#ifdef RANGECHECK
         if ((c >= cols()) || (r >= rows()))
         {
            MatrixException e("Invalid ConstMatrixBase index for rowCopy");
            GNSSTK_THROW(e);
         }
#endif
         Vector<T> toReturn(cols() - c);
         size_t i;
         for (i = c; i < cols(); i++)
            toReturn(i - c) = (*this)(r, i);
         return toReturn;
      }
 
      Vector<T> diagCopy(void) const
      {
         size_t i, n(cols());
         if(rows() < n) n = rows();
         Vector<T> toReturn(n);
         for(i=0; i<n; i++)
            toReturn(i) = (*this)(i,i);
         return toReturn;
      }
 
   protected:
      inline T constMatrixRef(size_t i, size_t j) const
      {
         const BaseClass& b = static_cast<const BaseClass&>(*this);
#ifdef RANGECHECK
         if ((i >= b.rows()) || (j > b.cols()))
         {
            MatrixException e("Invalid ConstMatrixBase index for ref");
            GNSSTK_THROW(e);
         }
#endif
         return b(i,j);
      }
   };
 
   template <class T, class BaseClass>
   class RefMatrixBase : public ConstMatrixBase<T, BaseClass>
   {
   public:
      explicit RefMatrixBase() {}
 
      virtual ~RefMatrixBase()
      {}
 
      T& operator() (size_t i, size_t j)
      { return static_cast<BaseClass*>(this)->operator()(i,j); }
 
      size_t size() const
      { return static_cast<const BaseClass*>(this)->size(); }
      size_t cols() const
      { return static_cast<const BaseClass*>(this)->cols(); }
      size_t rows() const
      { return static_cast<const BaseClass*>(this)->rows(); }
      BaseClass& zeroize()
      {
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i, j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               if (ABS(me(i,j)) < RefVectorBaseHelper::zeroTolerance)
                  me(i,j) = T(0);
         return me;
      }
      BaseClass& zeroizeRow(size_t r)
      {
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t j;
         for (j=0; j < me.cols(); j++)
             if (ABS(me(r, j)) < RefVectorBaseHelper::zeroTolerance)
                 me(r, j) = T(0);
         return me;
      }
      BaseClass& zeroizeCol(size_t c)
      {
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i;
         for (i=0; i < me.rows(); i++)
            if (ABS(me(i,c)) < RefVectorBaseHelper::zeroTolerance)
               me(i,c) = T(0);
         return me;
      }
 
 
 
         //MatBaseNewAssignOperator(assignFrom, =);
         //MatBaseNewAssignOperator(operator+=, +=);
         //MatBaseNewAssignOperator(operator-=, -=);
 
         //      MatBaseNewAssignOperator(assignFrom, =);
 
      template <class E> BaseClass& assignFrom(const ConstMatrixBase<T, E>& x)
      {
            //MatBaseArrayAssignMacro(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.rows() != me.rows() || x.cols() != me.cols()) {
            MatrixException e("Invalid dimensions for Matrix assignFrom(Matrix)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) = x(i,j);
         return me;
      }
      template <class E> BaseClass& assignFrom(const ConstVectorBase<T, E>& x)
      {
            //MatBaseArrayAssignMacroVecSource(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix assignFrom(Vector)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) = x[i*me.cols()+j];
         return me;
      }
      BaseClass& assignFrom(const std::valarray<T>& x)
      {
            //MatBaseArrayAssignMacroVecSource(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix assignFrom(valarray)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) = x[i*me.cols()+j];
         return me;
      }
      BaseClass& assignFrom(const T* x)
      {
            //MatBaseArrayAssignMacroVecSource(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) = x[i*me.cols()+j];          // no way to RANGECHECK on x[..]!
         return me;
      }
      BaseClass& assignFrom(T x)
      {
            //MatBaseAtomicAssignMacro(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) = x;
         return me;
      }
 
         //      MatBaseNewAssignOperator(operator+=, +=);
      template <class E> BaseClass& operator+=(const ConstMatrixBase<T, E>& x)
      {
            //MatBaseArrayAssignMacro(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.rows() != me.rows() || x.cols() != me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator+=(Matrix)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) += x(i,j);
         return me;
      }
      template <class E> BaseClass& operator+=(const ConstVectorBase<T, E>& x)
      {
            //MatBaseArrayAssignMacroVecSource(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator+=(Vector)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) += x[i*me.cols()+j];
         return me;
      }
      BaseClass& operator+=(const std::valarray<T>& x)
      {
            //MatBaseArrayAssignMacroVecSource(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator+=(valarray)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) += x[i*me.cols()+j];
         return me;
      }
      BaseClass& operator+=(const T* x)
      {
            //MatBaseArrayAssignMacroVecSource(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) += x[i*me.cols()+j];          // no way to RANGECHECK on x[..]!
         return me;
      }
      BaseClass& operator+=(T x)
      {
            //MatBaseAtomicAssignMacro(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) += x;
         return me;
      }
 
 
         //#define MatBaseNewAssignOperator(operator-=, -=)
      template <class E> BaseClass& operator-=(const ConstMatrixBase<T, E>& x)
      {
            //MatBaseArrayAssignMacro(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.rows() != me.rows() || x.cols() != me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator-=(Matrix)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) -= x(i,j);
         return me;
      }
      template <class E> BaseClass& operator-=(const ConstVectorBase<T, E>& x)
      {
            //MatBaseArrayAssignMacroVecSource(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator-=(Vector)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) -= x[i*me.cols()+j];
         return me;
      }
      BaseClass& operator-=(const std::valarray<T>& x)
      {
            //MatBaseArrayAssignMacroVecSource(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator-=(valarray)");
            GNSSTK_THROW(e);
         }
#endif
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) -= x[i*me.cols()+j];
         return me;
      }
      BaseClass& operator-=(const T* x)
      {
            //MatBaseArrayAssignMacroVecSource(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) -= x[i*me.cols()+j];          // no way to RANGECHECK on x[..]!
         return me;
      }
      BaseClass& operator-=(T x)
      {
            //MatBaseAtomicAssignMacro(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) -= x;
         return me;
      }
 
 
      BaseClass& operator*=(const T x)
      {
            //MatBaseAtomicAssignMacro(*=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) *= x;
         return me;
      }
 
      BaseClass& operator/=(const T x)
      {
            //MatBaseAtomicAssignMacro(/=);
         BaseClass& me = static_cast<BaseClass&>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) /= x;
         return me;
      }
 
      const BaseClass operator-()
      {
         const T x=T(-1);
         BaseClass me = static_cast<BaseClass>(*this);
         size_t i,j;
         for (i=0; i < me.rows(); i++)
            for (j=0; j < me.cols(); j++)
               me(i,j) *= x;
         return me;
      }
 
      BaseClass& swapRows(size_t row1, size_t row2)
      {
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if ( (row1 >= me.rows()) || (row2 >= me.rows()) )
         {
            MatrixException e("Invalid rows for swapRows");
            GNSSTK_THROW(e);
         }
#endif
         size_t i;
         T temp;
         for (i = 0; i < me.cols(); i++)
         {
            temp = me(row1, i);
            me(row1,i) = me(row2,i);
            me(row2,i) = temp;
         }
         return me;
      }
 
      BaseClass& swapCols(size_t col1, size_t col2)
      {
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if ( (col1 >= me.cols()) || (col2 >= me.cols()) )
         {
            MatrixException e("Invalid columns for swapCols");
            GNSSTK_THROW(e);
         }
#endif
         size_t i;
         T temp;
         for (i = 0; i < me.rows(); i++)
         {
            temp = me(i, col1);
            me(i, col1) = me(i, col2);
            me(i, col2) = temp;
         }
         return me;
      }
   };
 
   template <class T, class BaseClass>
   class MatrixSliceBase
   {
      size_t rowSize() const
      { return static_cast<const BaseClass*>(this)->rowSize(); }
      size_t rowStart() const
      { return static_cast<const BaseClass*>(this)->rowStart(); }
      size_t rowStride() const
      { return static_cast<const BaseClass*>(this)->rowStride(); }
      size_t colSize() const
      { return static_cast<const BaseClass*>(this)->colSize(); }
      size_t colStart() const
      { return static_cast<const BaseClass*>(this)->colStart(); }
      size_t colStride() const
      { return static_cast<const BaseClass*>(this)->colStride(); }
   protected:
      inline void matSliceCheck(size_t sourceRowSize,
                                size_t sourceColSize) const
      {
            //#ifdef RANGECHECK
         if (rowSize() > 0)
         {
            if ( (rowStart() >= sourceRowSize) ||
                 ((rowStart() + (rowSize()-1) * rowStride()) >= sourceRowSize))
            {
               MatrixException e("Invalid row range for slice");
               GNSSTK_THROW(e);
            }
         }
         if (colSize() > 0)
         {
            if ( (colStart() >= sourceColSize) ||
                 ((colStart() + (colSize()-1) * colStride()) >= sourceColSize))
            {
               MatrixException e("Invalid col range for slice");
               GNSSTK_THROW(e);
            }
         }
            //#endif
      }
   };
 
   template <class T, class BaseClass>
   class ConstMatrixSliceBase : public MatrixSliceBase<T, BaseClass>,
                                public ConstMatrixBase<T, BaseClass>
   {
   public:
      explicit ConstMatrixSliceBase() {}
   };
 
   template <class T, class BaseClass>
   class RefMatrixSliceBase : public MatrixSliceBase<T, BaseClass>,
                              public RefMatrixBase<T, BaseClass>
   {
   public:
      explicit RefMatrixSliceBase() {}
   };
 
 
}  // namespace
 
#include "MatrixBaseOperators.hpp"
 
#endif