STK_Matrix.h

Go to the documentation of this file.

/*--------------------------------------------------------------------*/
/*     Copyright (C) 2004-2007  Serge Iovleff

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

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this program; if not, write to the
    Free Software Foundation, Inc.,
    59 Temple Place,
    Suite 330,
    Boston, MA 02111-1307
    USA

    Contact : Serge.Iovleff@stkpp.org
*/

/*
 * Project:  stkpp::Arrays
 * Purpose:  Define the Matrix classes.
 * Author:   Serge Iovleff, serge.iovleff@stkpp.org
 *
 **/

#ifndef STK_MATRIX_H
#define STK_MATRIX_H

#include "STK_Vector.h"
#include "STK_Point.h"

#include "STK_Array2D.h"
#include "../../Sdk/include/STK_ContainerTraits.h"

namespace STK
{

typedef Array2D<Real> Matrix;

template<>
class Array2D<Real> : public ITArray2D< Real, Array2D<Real> >
{
  typedef ITArray2D<Real, Array2D<Real> > _ITArray2DType;

  public:
    Array2D( Range const& I = Range(), Range const& J = Range())
           : _ITArray2DType(I, J)
    {
      // initialize vertically the container
      this->initializeCols(J);
    }

    Array2D( Range const& I, Range const& J, Real const& v)
           : _ITArray2DType(I, J)
    {
      // initialize vertically the container
      this->initializeCols(J);
      // initialize with v
      for (Integer j=J.first(); j<=J.last(); j++)
      {
        Real* p(this->data(j));
        for (Integer i=I.first(); i<=I.last(); i++) p[i]= v;
      }
    }

    Array2D( const Array2D &T, bool ref=false)
           : _ITArray2DType(T, ref)
    { 
      if (!ref)
      {
        // initialize the Cols and Rows
        this->initializeCols(T.rangeHo());
        for (Integer j=T.firstCol(); j<=T.lastCol(); j++)
        {
          // ptr on the rows
          Real* p = this->data(j);
          Real* pt= T.data(j);
          for (Integer i=T.firstRow(); i<=T.lastRow(); i++) p[i]=pt[i];
        }
      }
    }

    Array2D( const Array2D& T, Range const& I, Range const& J)
           : _ITArray2DType(T, I, J)
    { ;}

    Array2D( Real** q, Range const& I, Range const& J)
           : _ITArray2DType(q, I, J)
    { ;}

    virtual ~Array2D()
    { ;}

    inline Real& elt(Integer i, Integer j)
    { return this->data(j)[i];}

    inline Real const& elt(Integer i, Integer j) const
    { return this->data(j)[i];}

    inline Array2D sub(Range const& I, Range const& J) const
    { return Array2D(*this, I, J);} 

    inline Vector col( Integer j) const
    { return Vector(this->data(j), this->rangeVe(), j);}

    inline Vector col(Range const& I, Integer j) const
    { return Vector(this->data(j), I, j);}

    inline Point row( Integer const& i) const
    { return Point(*this, this->rangeHo(), i);}

    inline Point row(Integer i, Range const& J) const
    { return Point(*this, J, i);}

    inline Range compRangeVe( Integer const& icol) const
    { return this->rangeVe();}

    void swapRows( Integer const& pos1, Integer pos2)
    {
#ifdef STK_BOUNDS_CHECK
      // check conditions
      if (this->firstRow() > pos1)
      { throw out_of_range("Array2D::swapRows(pos1, pos2) "
                                "this->firstRow() > pos1");
      }
      if (this->lastRow() < pos1)
      { throw out_of_range("Array2D::swapRows(pos1, pos2) "
                                "this->lastRow() < pos1");
      }
      if (this->firstRow() > pos2)
      { throw out_of_range("Array2D::swapRows(pos1, pos2) "
                                "this->firstRow() > pos2");
      }
      if (this->lastRow() < pos2)
      { throw out_of_range("Array2D::swapRows(pos1, pos2) "
                                "this->lastRow() < pos2");
      }
#endif
      // swap
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      {
        Real aux(this->asLeaf().elt(pos1, j));
        this->asLeaf().elt(pos1, j) = this->asLeaf().elt(pos2, j);
        this->asLeaf().elt(pos2, j) = aux;
      }
    }

    inline Array2D& operator=(const Array2D &T)
    {
//      // clear any allocated memory
//      this->clear();
//      // set AllocatorBase part
//      this->setRef(true);
//      this->setPtrData(T.ptrData(), T.rangeData());
//      // ITContainer2D part
//      this->setRange(T.rangeVe_, T.rangeHo_);
//      // Array2D part
//      rangeCols_ = T.rangeCols();
//      this->setCapacityHo( T.capacityHo());
      // Resize if necessary.
      if ( (this->sizeVe() != T.sizeVe())
         ||(this->sizeHo() != T.sizeHo())
         )
        this->resize(T.rangeVe(), T.rangeHo());
      // Copy without overlapping
      if (T.firstRow()>=this->firstRow())
      {
        if (T.firstCol()>=this->firstCol())
        {
          for ( Integer jt=T.firstCol(), j=this->firstCol()
              ; jt<=T.lastCol()
              ; j++, jt++
              )
          {
            Real *p =this->data(j), *pt =T.data(jt);
            for ( Integer it=T.firstRow(), i=this->firstRow()
                ; it<=T.lastRow()
                ; i++, it++
                )
              p[i] = pt[it];
          }
          return *this;
        }
        // T.firstCol()<this->firstCol()
        for ( Integer jt=T.lastCol(), j=this->lastCol()
            ; jt>=T.firstCol()
            ; j--, jt--
            )
        {
          Real *p =this->data(j), *pt =T.data(jt);
          for ( Integer it=T.firstRow(), i=this->firstRow()
              ; it<=T.lastRow()
              ; i++, it++
              )
            p[i] = pt[it];
        }
        return *this;
      }
      // T.firstRow()<this->firstRow()
      if (T.firstCol()>=this->firstCol())
      {
        for ( Integer jt=T.firstCol(), j=this->firstCol()
            ; jt<=T.lastCol()
            ; j++, jt++
            )
        {
          Real *p =this->data(j), *pt =T.data(jt);
          for ( Integer it=T.lastRow(), i=this->lastRow()
              ; it>=T.firstRow()
              ; i--, it--
              )
            p[i] = pt[it] ;
        }
        return *this;
      }
      // T.firstCol()<this->firstCol()
      for ( Integer jt=T.lastCol(), j=this->lastCol()
          ; jt>=T.firstCol()
          ; j--, jt--
          )
      {
        Real *p =this->data(j), *pt =T.data(jt);
        for ( Integer it=T.lastRow(), i=this->lastRow()
            ; it>=T.firstRow()
            ; i--, it--
            )
          p[i] = pt[it];
      }
      return *this;
    }

    inline Array2D& operator=(Real const& v)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      { Real *p =this->data(j);
        for (Integer i=this->firstRow(); i<=this->lastRow(); i++)
          p[i] = v;
      }
      return *this;
    }

    inline Array2D& operator+=(Real const& v)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      { Real *p =this->data(j);
        for (Integer i=this->firstRow(); i<=this->lastRow(); i++)
          p[i] += v;
      }
      return *this;
    }

    inline Array2D& operator-=(Real const& v)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      { Real *p =this->data(j);
        for (Integer i=this->firstRow(); i<=this->lastRow(); i++)
          p[i] += v;
      }
      return *this;
    }

    inline Array2D& operator*=(Real const& v)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      { Real *p =this->data(j);
        for (Integer i=this->firstRow(); i<=this->lastRow(); i++)
          p[i] += v;
      }
      return *this;
    }

    inline Array2D& operator/=(Real const& v)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); j++)
      { Real *p =this->data(j);
        for (Integer i=this->firstRow(); i<=this->lastRow(); i++)
          p[i] /= v;
      }
      return *this;
    }

    template< class Exp>
    inline Array2D& operator=(const Exp& rhs)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); ++j)
        for (Integer i=this->firstRow(); i<=this->lastRow(); ++i)
        {
          this->elt(i,j) = rhs(i,j);
        }
      return *this;
    }

    template< class Exp>
    inline Array2D& operator+=(const Exp& rhs)
    {
      for (Integer j=this->firstCol(); j<=this->lastCol(); ++j)
        for (Integer i=this->firstRow(); i<=this->lastRow(); ++i)
        {
          this->elt(i,j) += rhs(i,j);
        }
      return *this;
    }

    template< class Exp>
    inline Array2D& operator-=(const Exp& rhs)
    {
        for (Integer j=this->firstCol(); j<=this->lastCol(); ++j)
          for (Integer i=this->firstRow(); i<=this->lastRow(); ++i)
          {
            this->elt(i,j) -= rhs(i,j);
          }
      return *this;
    }

    template< class Exp>
    inline Array2D& operator/=(const Exp& rhs)
    {
        for (Integer j=this->firstCol(); j<=this->lastCol(); ++j)
          for (Integer i=this->firstRow(); i<=this->lastRow(); ++i)
          {
            this->elt(i,j) /= rhs(i,j);
          }
      return *this;
    }

    template< class Exp>
    inline Array2D& operator*=(const Exp& rhs)
    {
        for (Integer j=this->firstCol(); j<=this->lastCol(); ++j)
          for (Integer i=this->firstRow(); i<=this->lastRow(); ++i)
          {
            this->elt(i,j) *= rhs(i,j);
          }
      return *this;
    }
};

} // namespace STK

#endif
// STK_MATRIX_H