/Users/jab/Documents/Teaching/02585/GEL2_and_demos/GEL/src/LinAlg/Matrix.h

Go to the documentation of this file.

#if !defined(MATRIX_H_HAA_AGUST_2001)
#define MATRIX_H_HAA_AGUST_2001

#include <cassert>
#include <iostream>
#include <cmath>
#include "Vector.h"
#include "CGLA/ArithMatFloat.h"

namespace LinAlg
{

        template <class T>
        class CMatrixType 
        {
                friend class CVectorType<T>;

        private:
                int nRows;
                int nCols;
                int nElems;
                T* Data;
                T** Ilife;
        
                void SetScalar(const T& Scalar) 
                {
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        *Itt=Scalar;
                                        ++Itt;
                                }
                } 
        
                void CleanUp(void)
                {
                        if(Data!=NULL) 
                                delete [] Data; 
                        if(Ilife!=NULL) 
                                delete [] Ilife;
                };
        
                void Allocate(void)
                { 
                        Data= new T[nElems]; 
                        Ilife=new T*[nRows];    
                        for(int cRow=0;cRow<nRows;cRow++) 
                                Ilife[cRow]=&Data[cRow*nCols];
                };
        
        public:

                CMatrixType<T>() : nRows(0),nCols(0),nElems(0),Data(NULL),Ilife(NULL) {};

                CMatrixType<T>(const int Row,const int Col) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate();}

                CMatrixType<T>(const int Row,const int Col,const T& Scalar) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate(); SetScalar(Scalar);};

                CMatrixType<T>(const CMatrixType<T>& Rhs) : nRows(Rhs.nRows),nCols(Rhs.nCols),nElems(Rhs.nRows*Rhs.nCols) 
                { 
                        Allocate();
                        memcpy(Data,Rhs.Data,nElems*sizeof(T));
                }

                template <class VVT, class HVT, class MT, unsigned int ROWS>
                CMatrixType<T>(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs): 
                        nRows(Rhs.get_v_dim()),nCols(Rhs.get_h_dim()),
                        nElems(nRows*nCols) 
                { 
                        Allocate();
                        for(int i=0;i<nElems;++i)
                                Data[i] = Rhs.get()[i];
                }

                CMatrixType<T>(const CVectorType<T>& Rhs) : nRows(Rhs.Length()),nCols(1),nElems(Rhs.Length()) 
                { 
                        Allocate();
                        memcpy(Data,Rhs.Data,nElems*sizeof(T));
                }
                virtual ~CMatrixType<T>(){CleanUp();};

                CMatrixType<T>& operator=(const T& Rhs){ SetScalar(Rhs);return *this;}

    CMatrixType<T>& operator=(const CMatrixType<T>& Rhs)
                { 
                        Resize(Rhs.nRows,Rhs.nCols); 
                        memcpy(Data,Rhs.Data,nElems*sizeof(T));
                        return *this;
                }

                template <class VVT, class HVT, class MT, unsigned int ROWS>
                CMatrixType<T>& operator=(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs)
                { 
                        Resize(Rhs.get_v_dim(),Rhs.get_h_dim()); 
                        for(int i=0;i<nElems;++i)
                                Data[i] = Rhs.get()[i];
                        return *this;
                }

                CMatrixType<T>& operator=(const CVectorType<T>&Rhs)
                {
                        Resize(Rhs.Lenght(),1);
                        memcpy(Data,Rhs.Dat,nElems*sizeof(T));
                        return *this;
                }



                const int Rows(void) const {return nRows;};
                const int Cols(void) const {return nCols;};
                void Resize(const int Row,const int Col)
                {
                        if(Row!=nRows || Col!=nCols)
                                {
                                        CleanUp();
                                        nRows=Row;
                                        nCols=Col;
                                        nElems=Row*Col;
                                        Data= new T[nElems];
                                        Ilife=new T*[nRows];
                                        for(int cRow=0;cRow<nRows;cRow++)
                                                {
                                                        Ilife[cRow]=&Data[cRow*nCols];
                                                }
                                }
                }
        
                T* operator[](const int Row) 
                {
                        assert(Row>=0);
                        assert(Row<nRows);
                
                        return Ilife[Row];
                }
        
                const T* operator[](const int Row) const 
                {
                        assert(Row>=0);
                        assert(Row<nRows);
                
                        return Ilife[Row];
                }

                const T& get(const int Row,const int Col) const
                {
                        assert(Row>=0 && Row<nRows);
                        assert(Col>=0 && Col<nCols);
                        return Ilife[Row][Col];
                }

                void set(const int Row,const int Col,const T val) 
                {
                        assert(Row>=0 && Row<nRows);
                        assert(Col>=0 && Col<nCols);
                        Ilife[Row][Col]=val;
                }
        
                CMatrixType<T> operator+(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
                CMatrixType<T> operator-(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};
                CMatrixType<T> operator*(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret*=Rhs;};
                CMatrixType<T> operator/(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret/=Rhs;};
                CMatrixType<T>& operator+=(const T& Rhs)
                {

                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)+=Rhs;
                                        ++Itt;
                                }
                
                        return *this; 
                }
        
                CMatrixType<T>& operator-=(const T& Rhs)
                {
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)-=Rhs;
                                        ++Itt;
                                }
                
                        return *this; 
                }
        
                CMatrixType<T>& operator*=(const T& Rhs)
                {
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)*=Rhs;
                                        ++Itt;
                                }
                
                        return *this; 
                }
        
                CMatrixType<T>& operator/=(const T& Rhs)
                {
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)/=Rhs;
                                        ++Itt;
                                }
                
                        return *this; 
                }

                CMatrixType<T> operator+(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
                CMatrixType<T> operator-(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};

                CMatrixType<T>& operator+=(const CMatrixType<T>& Rhs)
                {
                        assert(nRows==Rhs.nRows);
                        assert(nCols==Rhs.nCols);

                        T* IttRhs=Rhs.Data;
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)+=(*IttRhs);
                                        ++Itt;
                                        ++IttRhs;
                                }

                        return *this; 
                }
        
                CMatrixType<T>& operator-=(const CMatrixType<T>& Rhs)
                {
                        assert(nRows==Rhs.nRows);
                        assert(nCols==Rhs.nCols);
                        T* IttRhs=Rhs.Data;
                        T* Itt=Data;
                        T* Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)-=(*IttRhs);
                                        ++Itt;
                                        ++IttRhs;
                                }
                        return *this; 
                }

                CMatrixType<T> operator* (const CMatrixType<T>& Rhs) const
                {
                        assert(nCols==Rhs.nRows);
                        CMatrixType<T> Ret(nRows,Rhs.nCols);

                        for(int cRow=0;cRow<nRows;cRow++)
                                {
                                        for(int cCol=0;cCol<Rhs.nCols;cCol++)
                                                {
                                                        T* Ret_ij=&(Ret[cRow][cCol]);
                                                        *Ret_ij=0;
                                                        for(int cI=0;cI<nCols;cI++)
                                                                {
                                                                        (*Ret_ij)+=(Ilife[cRow])[cI]*Rhs[cI][cCol];
                                                                }
                                                }
                                }
                        return Ret;
                }

                CVectorType<T> operator* (const CVectorType<T>& Rhs) const 
                {
                        assert(nCols==Rhs.nElems);
                        CVectorType<T> Ret(nRows,0);

                        T* ThisItt=Data;
                        T* Stop=&(Data[nElems]);
                        T* RhsItt;
                        T* RetItt=Ret.Data;
                        while(ThisItt!=Stop)
                                {
                                        RhsItt=Rhs.Data;
                                        for(int cCol=nCols;cCol>0;--cCol)
                                                {
                                                        (*RetItt)+=(*RhsItt)*(*ThisItt);
                                                        ++ThisItt;
                                                        ++RhsItt;
                                                }
                                        ++RetItt;
                                }
                        return Ret;
                }
        

                CMatrixType<T>& Transpose(void)
                {
                        T* NewData =new T[nElems];
                        T*RowP;
                        int cRow;
                        for(cRow=0;cRow<nRows;cRow++)
                                {
                                        RowP=Ilife[cRow];
                                        for(int cCol=0;cCol<nCols;cCol++)
                                                {
                                                        NewData[cCol*nRows+cRow]=RowP[cCol];
                                                }
                                }
                
                        CleanUp();
                
                        Data=NewData;
                        int temp=nRows;
                        nRows=nCols;
                        nCols=temp;
                        Ilife=new T*[nRows];    
                        for(cRow=0;cRow<nRows;cRow++) 
                                {
                                        Ilife[cRow]=&Data[cRow*nCols];
                                }
                
                        return *this;
                }

                void Transposed(CMatrixType<T>& AT) const
                {
                        AT.Resize(nCols,nRows);
                        for(int cRow=nRows-1;cRow>=0;--cRow)
                                {
                                        for(int cCol=nCols-1;cCol>=0;--cCol)
                                                {
                                                        AT[cCol][cRow]=Ilife[cRow][cCol];
                                                }
                                }
                }

                CMatrixType<T> Transposed(void) const
                {
                        CMatrixType<T> Ret(nCols,nRows);
                        Transposed(Ret);
                        return Ret;
                }



                void ElemMult(const CMatrixType<T>& Rhs)
                {
                        T*Itt=Data;
                        T*Stop=&(Data[nElems]);
                        T*RhsItt=Rhs.Data;
                        while(Itt!=Stop)
                                {
                                        (*Itt)*=(*RhsItt);
                                        ++Itt;
                                        ++RhsItt;
                                }
                }

                void ElemDiv(const CMatrixType<T>& Rhs)
                {
                        T*Itt=Data;
                        T*Stop=&(Data[nElems]);
                        T*RhsItt=Rhs.Data;
                        while(Itt!=Stop)
                                {
                                        (*Itt)/=(*RhsItt);
                                        ++Itt;
                                        ++RhsItt;
                                }
                }

                void ElemSqr(void)
                {
                        T*Itt=Data;
                        T*Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)*=(*Itt);
                                        ++Itt;
                                }
                }

                void ElemSqrt(void)
                {
                        T*Itt=Data;
                        T*Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        (*Itt)=sqrt(*Itt);
                                        ++Itt;
                                }
                }


                const T NormMax(void) const
                {
                        T Ret=0;
                        T Abs;
                        const T*Itt=Data;
                        const T*Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        Abs=(*Itt)>0?(*Itt):-(*Itt);
                                        if(Ret<Abs)
                                                Ret=Abs;

                                        ++Itt;
                                }
                        return Ret;
                }

                const T NormFrobenius(void) const
                {
                        T Ret=0;
                        const T*Itt=Data;
                        const T*Stop=&(Data[nElems]);
                        while(Itt!=Stop)
                                {
                                        Ret+=(*Itt)*(*Itt);
                                        ++Itt;
                                }
                        return sqrt(Ret);
                }



                void GetRow(CVectorType<T>& Row,const int cRow) const
                {
                        Row.Resize(nCols);
                        memcpy(Row.Data,Ilife[cRow],nCols*sizeof(T));
                }

                void SetRow(CVectorType<T>& Row,const int cRow)
                {
                        assert(Row.Length()==nCols);
                        memcpy(Ilife[cRow],Row.Data,nCols*sizeof(T));
                }

                void RowTransfere(CMatrixType<T>& A, const int from,const int to) const
                {
                        assert(A.Cols()==nCols);
                        memcpy(A.Ilife[to],Ilife[from],nCols*sizeof(T));
                }

                void GetCol(CVectorType<T>&Col, const int nCol) const
                {
                        Col.Resize(nRows);
                        const T* Itt=&Data[nCol];
                        T* VecItt=&Col[0];
                        for(int cRow=0;cRow<nRows;cRow++)
                                {
                                        *VecItt=*Itt;
                                        ++VecItt;
                                        Itt+=nCols;
                                }
                }

                void SetCol(CVectorType<T>&Col, const int nCol)
                {
                        assert(nRows==Col.Length());
                        assert(nCol<nCols);

                        T* Itt=&Data[nCol];
                        T* VecItt=&Col[0];
                        for(int cRow=0;cRow<nRows;cRow++)
                                {
                                        *Itt=*VecItt;
                                        ++VecItt;
                                        Itt+=nCols;
                                }
                }

                void GetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol, const int RowDim, const int ColDim) const
                {
                        assert(StartRow+RowDim<=nRows);
                        assert(StartCol+ColDim<=nCols);
                        A.Resize(RowDim,ColDim);
                        const T* Itt=&(Ilife[StartRow][StartCol]);

                        for(int cRow=0;cRow<A.Rows();cRow++)
                                {
                                        memcpy(A[cRow],Itt,A.Cols()*sizeof(T));
                                        Itt+=nCols;
                                }
                }

                void SetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol)
                {
                        assert(A.Rows()+StartRow<=nRows);
                        assert(A.Cols()+StartCol<=nCols);
                        T* Itt=&(Ilife[StartRow][StartCol]);

                        for(int cRow=0;cRow<A.Rows();cRow++)
                                {
                                        memcpy(Itt,A[cRow],A.Cols()*sizeof(T));
                                        Itt+=nCols;
                                }
                }




        };

        template<class T>
        inline CMatrixType<T> operator+(const T& Lhs,const CMatrixType<T>&Rhs ) 
        {
                return Rhs+Lhs;
        }

        template<class T>
        inline CMatrixType<T> operator*(const T& Lhs,const CMatrixType<T>&Rhs ) 
        {
                return Rhs*Lhs;
        }

        template <class T>
        std::ostream& operator<<(std::ostream &s, const CMatrixType<T> &A)
        {
    int nRows=A.Rows();
    int nCols=A.Cols();
        
    for (int cRow=0; cRow<nRows; cRow++)
                        {
        for (int cCol=0; cCol<nCols; cCol++)
                                        {
            s << A[cRow][cCol] << " ";
                                        }
        s << "\n";
                        }
    return s;
        }

        template <class T>
        std::istream& operator>>(std::istream &s, CMatrixType<T> &A)
        {
    int nRows;
                int nCols;

    s >> nRows >> nCols;

    if ( nRows!=A.Rows() || nCols!=A.Cols() )
                        A.Resize(nRows,nCols);

                T* RowP;


    for (int cRow=0; cRow<nRows;cRow++)
                        {
                                RowP=A[cRow];
        for (int cCol=0;cCol<nCols;cCol++)
                                        {
            s >>  RowP[cCol];
                                        }
                        }

    return s;
        }

        typedef CMatrixType<double> CMatrix;
}
#endif // !defined(MATRIX_H_HAA_AGUST_2001)