MatrixBase.hpp

Go to the documentation of this file.

#pragma ident "$Id: MatrixBase.hpp 2006 2009-07-02 14:44:05Z afarris $"



#ifndef GPSTK_MATRIX_BASE_HPP
#define GPSTK_MATRIX_BASE_HPP

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk 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 2.1 of the License, or
//  any later version.
//
//  The GPSTk 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 GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

#include "Vector.hpp"

namespace gpstk
{

   NEW_EXCEPTION_CLASS(MatrixException, Exception);
   NEW_EXCEPTION_CLASS(SingularMatrixException, MatrixException);

   template <class T, class BaseClass>
   class ConstMatrixBase
   {
   public:
      explicit 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
         throw(MatrixException)
         { 
#ifdef RANGECHECK
            if ((c >= cols()) || (r >= rows()))
            {
               MatrixException e("Invalid ConstMatrixBase index for colCopy");
               GPSTK_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
         throw(MatrixException)
         { 
#ifdef RANGECHECK
            if ((c >= cols()) || (r >= rows()))
            {
               MatrixException e("Invalid ConstMatrixBase index for rowCopy");
               GPSTK_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;
         }

   protected:
      inline T constMatrixRef(size_t i, size_t j) const
         throw(MatrixException)
         {
            const BaseClass& b = static_cast<const BaseClass&>(*this);
#ifdef RANGECHECK
            if ((i >= b.rows()) || (j > b.cols()))
            {
               MatrixException e("Invalid ConstMatrixBase index for ref");
               GPSTK_THROW(e);
            }
#endif
            return b(i,j);
         }
   };

   template <class T, class BaseClass>
   class RefMatrixBase : public ConstMatrixBase<T, BaseClass>
   {
   public:
      explicit 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;
         }

// this macro expansion technique has 2 major faults: 1. there can be no RANGECHECK,
// 2. since there is no range checking, one can have a mis-match in dimensions in the
// l-value vs the r-value, which makes no sense.
/*
#define MatBaseArrayAssignMacro(func) \
   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) func x(i,j); \
   return me;
*/              
/*
#define MatBaseArrayAssignMacroVecSource(func) \
   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) func x[i*me.cols()+j]; \
   return me; 
*/
/*
#define MatBaseAtomicAssignMacro(func) \
   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) func x; \
   return me;
*/
/*
//#define MatBaseNewAssignOperator(funcName, op) \
*/


/*
//   template <class E> BaseClass& funcName(const ConstMatrixBase<T, E>& x) 
//      { MatBaseArrayAssignMacro(op); } \
*/
/*
//   template <class E> BaseClass& funcName(const ConstVectorBase<T, E>& x) 
//      { MatBaseArrayAssignMacroVecSource(op); } \
*/
/*
//   BaseClass& funcName(const std::valarray<T>& x) \
//      { MatBaseArrayAssignMacroVecSource(op); } \
*/
/*
//   BaseClass& funcName(const T* x) \
//      { MatBaseArrayAssignMacroVecSource(op); } \
*/
/*
//   BaseClass& funcName(T x) \
//      { MatBaseAtomicAssignMacro(op); }
*/              
      //MatBaseNewAssignOperator(assignFrom, =);
      //MatBaseNewAssignOperator(operator+=, +=);
      //MatBaseNewAssignOperator(operator-=, -=);
//--------------------------------------------------------------------------
//      MatBaseNewAssignOperator(assignFrom, =);
   template <class E> BaseClass& assignFrom(const ConstMatrixBase<T, E>& x)
      throw(MatrixException)
      {
         //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)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix assignFrom(Vector)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix assignFrom(valarray)");
            GPSTK_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)
      throw(MatrixException)
      {
         //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)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator+=(Vector)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(+=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator+=(valarray)");
            GPSTK_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)
      throw(MatrixException)
      {
         //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)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator-=(Vector)");
            GPSTK_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)
      throw(MatrixException)
      {
         //MatBaseArrayAssignMacroVecSource(-=);
         BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
         if(x.size() != me.rows() * me.cols()) {
            MatrixException e("Invalid dimensions for Matrix operator-=(valarray)");
            GPSTK_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;
         }

         // unary minus: multiplies each element in this matrix by -1.
      //BaseClass& operator-()
      //   {
      //      const T x=T(-1);
      //      MatBaseAtomicAssignMacro(*=);
      //   }
      // unary minus must not return an l-value

      const BaseClass operator-() const
      {
         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) 
         throw(MatrixException)
         {
            BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
            if ( (row1 >= me.rows()) || (row2 >= me.rows()) )
            {
               MatrixException e("Invalid rows for swapRows");
               GPSTK_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) 
         throw(MatrixException)
         {
            BaseClass& me = static_cast<BaseClass&>(*this);
#ifdef RANGECHECK
            if ( (col1 >= me.cols()) || (col2 >= me.cols()) )
            {
               MatrixException e("Invalid columns for swapCols");
               GPSTK_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
         throw(MatrixException)
         {
//#ifdef RANGECHECK
            if (rowSize() > 0)
            {
               if ( (rowStart() >= sourceRowSize) || 
                    ((rowStart() + (rowSize()-1) * rowStride()) >= sourceRowSize))
               {
                  MatrixException e("Invalid row range for slice");
                  GPSTK_THROW(e);
               }
            }
            if (colSize() > 0)
            {
               if ( (colStart() >= sourceColSize) ||
                    ((colStart() + (colSize()-1) * colStride()) >= sourceColSize))
               {
                  MatrixException e("Invalid col range for slice");
                  GPSTK_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

---