TMatrix.cpp

Go to the documentation of this file.

#include <iomanip>

template <class TRing>
inline void TMatrix<TRing>::init( )
{
        assert (rows_m>=0);
        assert (cols_m>=0);

        data            = NULL;

        size_t size = rows_m;
        size*=cols_m;

        if ( (size) > 0 )
        {
                data            = new TNum[ size ];
                assert( data!=NULL );
        }
        
        pcolumn = NULL;
        
        if (rows_m>0)
                pcolumn = new TNum*[rows_m];

        for (unsigned int i=0; i<this->rows_m; i++)
                pcolumn[i] = &data[ getDataPos(i,0)];
}


//---------------------------------------construction /destruction------------------------------------------------------
template <class TRing>
TMatrix<TRing>::TMatrix(unsigned int _zeilen, unsigned  int _spalten,  
                        const TRing* _ring,     std::string _name  ):   rows_m(_zeilen), 
                                                                        cols_m(_spalten), 
                                                                        ring(_ring), 
                                                                        name(_name)
{
        
        init();
        // ab hier ist es anders.
        fillZero();
}

template <class TRing>
TMatrix<TRing>::TMatrix(unsigned int _zeilen, unsigned  int _spalten,  
                        const TRing* _ring, TVector<TRing> & vec,std::string _name  ):  rows_m(_zeilen), 
                                                                                        cols_m(_spalten), 
                                                                                        ring(_ring), 
                                                                                        name(_name)
{
        init();

        size_t size = rows_m;
        size*=cols_m;
        // ab hier ist es anders.
        assert(size==vec.getSize() )    ;

        for ( size_t i=0; i<getSize(); i++ )
                        data[i] = vec.getVal(i);

}




template <class TRing>
template <class TemplateMatrix>
TMatrix<TRing>::TMatrix(const TemplateMatrix &mat, const TRing* _ring):
                                                                        rows_m( mat.getRowNum() ), 
                                                                        cols_m( mat.getColNum() ), 
                                                                        ring( _ring), 
                                                                        name( mat.getName() )
                                                                        
{
        init();
        
        
        // ab hier ist es anders.
        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<cols_m; col++)
                {
                        data[ getDataPos(row,col) ] = ring->Convert( mat.getVal(row,col) );
                }
        }
}

/*
template <class TRing> template <class _istream>
TNum Convert(int entry)
{
        if ( ring )

                ring->Convert();
        else
                return TNum(entry);
}*/


template <class TRing> template <class _istream>
void TMatrix<TRing>::createFromStream(_istream & _dataStream)
{
        #ifdef DEBUG
                std::cerr << " TMatrix::createFromStream(): " << std::endl << std::endl; 
        #endif

        std::string matrixStr = extractNextBracedData( _dataStream );

        std::stringstream matrixStream(matrixStr);

        //matrixStream << ;

        rows_m = countSubGroups( matrixStream );

        std::cerr << "rows_m"<<rows_m << std::endl;
        

        extractChar(matrixStream, '{');

        std::string matrixRow = extractNextBracedData( matrixStream );

        std::stringstream smatrixRow(matrixRow);
        

        cols_m = countElements(smatrixRow);

        std::cerr << "cols_m" << cols_m << std::endl;
        
        // ab hier ist es gleich

        init();

        // ab hier ist es anders

        int elemNr = 0;

        matrixStream.clear();
        matrixStream.seekg(0);
        
        extractChar(matrixStream, '{');
        
        std::cerr << "matrixStream" << matrixStream.str() << std::endl;

        for (unsigned int row=0;  row<rows_m ;  row++)
        {
                matrixRow   = extractNextBracedData(matrixStream);

                std::cerr << "matrixRow" << matrixRow << std::endl;     
                std::stringstream       smatrixRow2(matrixRow) ;


                extractChar(smatrixRow2, '{');

                for (unsigned int col=0;  col<cols_m ;  col++)
                {
                        int currEntry = 0;
                        smatrixRow2 >> currEntry;
                
                        data[ elemNr ] = ring->Convert(currEntry);
                        elemNr++;

                        if (col !=cols_m-1 ) 
                                extractChar( smatrixRow2, ',' );
                        else extractChar( smatrixRow2, '}' );
                }
                if (row !=rows_m-1 ) 
                                extractChar( matrixStream, ',' );
                else extractChar( matrixStream, '}' );
        }
        std::cerr << "matrix initialization OK" << std::endl;
}



template <class TRing> 
TMatrix<TRing>::TMatrix(string matrixStr,  const TRing* _ring, std::string _name)   :  ring(_ring),  name(_name)
{
        std::stringstream  matrixStream;

        matrixStream << matrixStr;

        createFromStream( matrixStream );
}

template <class TRing>
 template <class _istream>
TMatrix<TRing>::TMatrix(  const TRing* _ring, _istream& dataStream, std::string _name)  :       ring( _ring ), 
                                                                                                name( _name )
{
        createFromStream(dataStream);
}


template <class TRing>
  TMatrix<TRing>::TMatrix(const TMatrix<TRing>& mat):   rows_m ( mat.rows_m  ) ,
                                                        cols_m ( mat.cols_m ) ,
                                                        ring (mat.ring) 
                                                        
{

                // hier kann man es auch gleich gestalten wie oben!
                init();

                for ( size_t i=0; i<(size_t)rows_m*cols_m; i++ )
                        data[i] = mat.data[i];          
}



template <class TRing>
TMatrix<TRing>::~TMatrix()
{
        if ( data!=NULL )
                delete[] data;

        if ( pcolumn!=NULL )
                delete[] pcolumn;
};


//---------------------------------------safety------------------------------------------------------

template <class TRing> 
size_t  TMatrix<TRing>::getSize() const 
{       
        size_t size =rows_m;
        size *=cols_m;
        return size;
};

template <class TRing> 
TVector<TRing>*         TMatrix<TRing>::convertToVector  () const
{
        TVector<TRing>* matrixAsVector = new TVector<TRing>(getSize(),ring,name);
        
        for ( size_t i=0; i<getSize(); i++ )
                        matrixAsVector->setVal(i,data[i] );

        return matrixAsVector;
}


template <class TRing> 
void  TMatrix<TRing>::checkBounds(unsigned int row, unsigned int col) const
{
        if ( row>=rows_m  || col>=cols_m)
                        {
                                std::cerr << "TMatrix: bound error !"           << std::endl;
                                std::cerr << "row " << row << " rows_m " << rows_m      << std::endl;
                                std::cerr << "col " << col << " cols_m " << cols_m      << std::endl;   
                                exit(0);
                        }
}


//---------------------------------------init------------------------------------------------------


template <class TRing>
void TMatrix<TRing>::fillZero()
{               
        // http://de.wikibooks.org/wiki/C%2B%2B-Programmierung:_Standard_Template_Library

        //std::fill(data, data+rows_m*cols_m, TNum::Zero);

        
        for (size_t pos=0; pos< rows_m*cols_m; pos++)
        {
                data[ pos ] = TNum::Zero;
        }
        //memset(data ,0x00, rows_m*cols_m*sizeof(TNum) );
}

template <class TRing>
void TMatrix<TRing>::randomInit(long * _pRandomSeed )
{
        for (size_t pos=0; pos< rows_m*cols_m; pos++)
        {
                data[ pos ] = ring->ConvertScalar( random(_pRandomSeed, ring->getCharacteristic()-1)   );
        }
}




//------------------------------------data access ----------------------------------------------------------------------

template <class TRing>
inline void TMatrix<TRing>::setVal(const unsigned int row, const unsigned int col, const TNum z)
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif
        data[ getDataPos(row, col)]=z;
}


template <class TRing>
inline typename TRing::ElementType TMatrix<TRing>::getVal(const unsigned int row,const unsigned int col) const
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif
        return data[ getDataPos(row, col)];
}


template <class TRing>
inline const typename TRing::ElementType& TMatrix<TRing>::getConstValRef(const unsigned int row,const unsigned int col) const
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif
        return data[ getDataPos(row, col) ];
}


template <class TRing>
inline typename TRing::ElementType & TMatrix<TRing>::getValRef(const unsigned int row,const unsigned int col) 
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif
        return data[ getDataPos(row, col) ];
}

template <class TRing>
inline typename TRing::ElementType* TMatrix<TRing>::getValAddr(const unsigned int row,const unsigned int col) 
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif
        return &( data[  getDataPos(row, col) ] );
}

template <class TRing>
typename TRing::ElementType& TMatrix<TRing>::getRowRef(unsigned int row)
{
        assert(row>0 && row<rows_m );
        return &data[ getDataPos(row, 0)];
}






//------------------------------------operators----------------------------------------------------------------------


template <class TRing>
TMatrix<TRing>*  TMatrix<TRing>::operator*(const TMatrix<TRing>& mat) const
{

        assert( cols_m==mat.rows_m );
        
        TMatrix<TRing>* erg = new TMatrix<TRing>(rows_m, mat.cols_m, this->ring);


        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<mat.cols_m; col++)
                {
                        TNum val = TNum::Zero;
                        for (unsigned int i=0; i<mat.rows_m; i++)
                        {
                                //   getVal(row,i);
                                mat.getVal(i,col);
                                ring->addInPlace( val  ,ring->multiply( getVal(row,i), mat.getVal(i,col)) );
                        }
                        erg->setVal( row, col, val);
                }
        }
        return erg;
}


template <class TRing>
TMatrix<TRing>*  TMatrix<TRing>::operator|(const TMatrix& mat)  const
{
        assert(this->getRowNum()==mat.getRowNum() );

        TMatrix<TRing>* erg = new TMatrix<TRing>(rows_m, cols_m+ mat.cols_m , this->ring);

        for (unsigned int row=0;row<rows_m; row++ )
        {
                for (unsigned int col=0; col< cols_m; col++)
                {
                        erg->setVal(row,col,getVal(row,col));
                }
                for (unsigned int col=0; col< mat.cols_m; col++ )
                {
                        erg->setVal(row, this->cols_m + col, mat.getVal(row,col) );
                }
        }
        return erg;
}



template <class TRing>
TMatrix<TRing>*         TMatrix<TRing>::getSubMatrix(unsigned int firstrow, unsigned int firstcol,unsigned int lastrow,unsigned int lastcol)    const
{
        assert(  firstrow<rows_m );     assert(  lastrow<rows_m );
        assert(  firstcol<cols_m );     assert(  lastcol<cols_m );
        
        unsigned int newrows = lastrow-firstrow + 1;
        unsigned int newcols = lastcol-firstcol + 1;

        TMatrix<TRing>* erg = new TMatrix<TRing>(newrows, newcols, this->ring);

        for (unsigned int row=firstrow; row<=lastrow; row++ )
        {
                for (unsigned int col=firstcol; col<= lastcol; col++)
                {
                        erg->setVal( row-firstrow, col-firstcol, getVal(row,col) );
                }
        }
        return erg;
}


template <class TRing>
TVector<TRing>*  TMatrix<TRing>::operator*(TVector<TRing>& vec) const
{
        
        assert(vec.getSize()==(size_t)cols_m);
        
        TVector<TRing>* erg = new TVector<TRing>(rows_m, this->ring);

        for (unsigned int row=0;row<rows_m;row++)
        {
                TNum val=TNum::Zero;
                for (size_t i=0; i<vec.getSize(); i++)
                {
                        TNum val1 = getVal(row,i);
                        TNum val2 = vec.getVal(i);

                        ring->multiply(val1,val2);
                        ring->addInPlace(val  ,ring->multiply( getVal(row,i), vec.getVal(i)) );
                }
                erg->setVal( row, val);
        }
        return erg;
}


template <class TRing>
bool TMatrix<TRing>::operator==(const TMatrix& mat) const
{       
        if ( rows_m != mat.getRowNum() || cols_m != mat.getColNum() )
                return false;

        
        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<cols_m; col++)
                {
                        if ( getVal(row,col) != mat.getVal(row,col) )
                                return false;
                }
        }
        return true;
}

template <class TRing>
bool TMatrix<TRing>::isIdentity() const
{       
        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<cols_m; col++)
                {
                        if (col!=row)
                        {
                                if ( getVal(row,col) != typename TRing::ElementType(0) )
                                        return false;
                        }
                        else
                        {
                                if ( getVal(row,col) != typename TRing::ElementType(1) )
                                        return false;
                        }
                        
                }
        }
        return true;
}


template <class TRing>
TMatrix<TRing>&  TMatrix<TRing>::operator=(const TMatrix<TRing>& mat)
{

        if (this!=&mat)
        {
                assert( ring==mat.ring );
                delete [] data;

                cols_m  = mat.cols_m;
                rows_m  = mat.rows_m;

                data = new TNum[rows_m*cols_m];

                for (size_t i=0; i<rows_m*cols_m; i++)
                        data[i] = mat.data[i];

                pcolumn = NULL;
                pcolumn = new TNum*[rows_m];

                for (unsigned int i=0; i<this->rows_m; i++)
                        pcolumn[i] = &data[ getDataPos(i,0) ];

                assert( data!=NULL );
        }
        return *this;
}


template <class TRing>
bool   TMatrix<TRing>::isZero() const
{
        for (size_t i=0; i<rows_m*cols_m; i++)
        {
                if (   data[i].isNotZero()  ) 
                        return false;
        }
        return true;
}

template <class TRing>
bool   TMatrix<TRing>::isNotZero() const
{
        for (size_t i=0; i<rows_m*cols_m; i++)
        {
                if (   data[i].isNotZero()  ) 
                        return true;
        }
        return false;
}

//------------------------------------extract data ----------------------------------------------------------------------

template <class TRing>
BlockMatrix<TRing>*  TMatrix<TRing>::getBlocked(unsigned int rowBlockSize, unsigned int colBlockSize) const
{

        BlockMatrix<TRing>*     erg      = new BlockMatrix<TRing>(rowBlockSize, colBlockSize, rows_m ,cols_m, ring);

        unsigned int Block = 0;

        for (register unsigned int h=0; h < erg->getBlocksPerCol(); h++) // run over y-coordonate of the blocks
        {
                unsigned int upperMostBlockRow = erg->getRowBlockSize()*h;

                for (register unsigned int l=0; l< erg->getBlocksPerRow(); l++) // run over x-coordonate of the blocks
                { 
                        unsigned int leftMostBlockCol=erg->getColBlockSize()*l;

                        unsigned int BlockPos=0;

                        // run over y-coordonate in the block
                        for (unsigned int row = upperMostBlockRow; row< upperMostBlockRow + erg->getRowBlockSize(); row++)
                        {
                                // run over x-coordonate in the block
                                for (unsigned int col=leftMostBlockCol; col<leftMostBlockCol+erg->getColBlockSize(); col++)
                                {
                                        erg->setBlockVal(Block,BlockPos, getVal(row,col) );
                                        BlockPos++;
                                }
                        }
                        Block++;
                }
                
        }
        return erg;
        
}


template <class TRing>
TMatrix<TRing>&          TMatrix<TRing>::transposeInPlace()
{

        data    = NULL;

        size_t tmp=rows_m;
        tmp*=cols_m;

        TNum* newData = new TNum[ tmp ];

        assert(newData!=NULL);


        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<cols_m; col++)
                {
                        newData[ col*rows_m+ row ] = data[ getDataPos(row,col) ];
                }
        }
         tmp=rows_m;

        rows_m=cols_m;

        cols_m=tmp;

        delete[] data;
        data = newData;

        delete[] pcolumn;

        pcolumn = NULL;
        pcolumn = new TNum*[rows_m];

        for (unsigned int i=0; i<this->rows_m; i++)
                pcolumn[i] = &data[ getDataPos(i,0) ];

        return *this;
}


template <class TRing>
TMatrix<TRing>*  TMatrix<TRing>::getTransposed() const
{
        TMatrix<TRing>* erg = new TMatrix<TRing>(cols_m, rows_m, this->ring);

        for (unsigned int row=0; row<rows_m; row++)
        {
                for (unsigned int col=0; col<cols_m; col++)
                {
                        erg->setVal(col, row, getVal(row,col));
                }
        }
        return erg;
        
}


template <class TRing>
typename TRing::ElementType* TMatrix<TRing>::getCol(unsigned int _col) const
{
        assert(_col>=0 && _col<cols_m );
        TNum* columndata= new TNum[rows_m];
        for (unsigned int i=0;i<rows_m;i++)
                columndata[i]=data[  getDataPos(i, _col) ];
        return columndata;
}

template <class TRing>
typename TRing::ElementType* TMatrix<TRing>::getRow(unsigned int _row) const
{
        assert(_row>=0 && _row<rows_m );
        TNum* rowdata= new TNum[cols_m];
        for (unsigned int i=0;i<cols_m;i++)
                rowdata[i]=data[  getDataPos(_row, i) ];

        return rowdata;
}

template <class TRing>
TVector< TRing>* TMatrix<TRing>::getColAsVector(unsigned int _col) const
{
        assert(_col>=0 && _col<cols_m );
        TVector< TRing>* colVector= new TVector<TRing>(rows_m,ring, "column");
        for (unsigned int i=0;i<rows_m;i++)
                colVector->setVal(i, data[  getDataPos(i, _col) ] );
        return colVector;
}

template <class TRing>
TVector< TRing>* TMatrix<TRing>::getRowAsVector(unsigned int _row) const
{
        assert(_row>=0 && _row<rows_m );
        TVector< TRing>* rowVector= new TVector<TRing>(cols_m,ring, "row");
        for (unsigned int i=0;i<cols_m;i++)
                rowVector->setVal(i, data[  getDataPos(_row, i) ] );

        return rowVector;
}

template <class TRing>
typename TRing::ElementType* TMatrix<TRing>::getDiag() const
{
        assert(rows_m==cols_m);

        TNum* diag= new TNum[cols_m];
        for (unsigned int i=0; i<cols_m; i++)
                diag[i] = data[ getDataPos(i, i) ];

        return diag;
}

template <class TRing>
TVector<TRing>* TMatrix<TRing>::getTrace() const
{
        assert(rows_m==cols_m);

        TVector<TRing>* erg = new TVector<TRing>(cols_m, this->ring);

        for (unsigned int i=0; i<cols_m; i++)
                erg->setVal( i,data[ getDataPos(i, i) ] );

        return erg;
}

template <class TRing>
vector<bool>    TMatrix<TRing>::getFreeVariables()
{
        vector<bool> varFree(cols_m); 

        unsigned int row, col;

        for(row=0, col=0; col<cols_m ; col++)
        {       
                if (row>=rows_m)
                {       
                        varFree[col]=true;
                }
                else
                {
                        if (this->pcolumn[row][col].isZero() )
                        {
                                varFree[col]=true;
                        }
                        else 
                        {
                                row++;
                                varFree[col]=false;
                        }
                }
                
        }
        return varFree;
}

template <class TRing>
vector<bool>    TMatrix<TRing>::getLeadingVariables()
{
        
        vector<bool> leadingVariables(rows_m); 

        unsigned int row, col;

        int lastrow=0;
        for(row=0, col=0; col<cols_m ; col++)
        {       
                if (row>=rows_m)
                {       
                        //varFree[col]=true;
                }
                else
                {
                        if (this.pcolumn[row][col].isZero() )
                        {
                                //varFree[col]=true;
                        }
                        else 
                        {
                                row++;
                                //varFree[col]=false;
                                lastrow=row+1;
                                leadingVariables[row]=col;
                        }
                }
                
        }
        leadingVariables.resize(lastrow);
        return leadingVariables;
}
//---------------------------------------compute properties------------------------------------------------------
// http://www.math.lsa.umich.edu/~hochster/419/ker.im.html
// ahhh schlaue Leute: 
// http://www.math.tu-dresden.de/~ganter/inf2004/Matrixslides.pdf
// hier steht, wie man den Kern berechnet, wenn man eine QR-Zerlegung besitzt:
// http://en.wikipedia.org/wiki/Left_null_space
// DAS ist die beste Beschreibung: http://www.jstor.org/stable/2686433?seq=3
// Matrix auf 'reduced row echolon' -Form bringen und dann die führenden Einsen auf die Diagonale setzen,
// und dann die Einheitsmatrix abziehen!
// #
//# Subspaces and Echelon Forms
//# David C. Lay
//# The College Mathematics Journal, Vol. 24, No. 1 (Jan., 1993), pp. 57-62 
template <class TRing>
        TMatrix<TRing>* TMatrix<TRing>::getKernel() const
{
        
        TMatrix<TRing> matCopy(*this);

        //std::cerr << "cols_m - getRank(true)" <<  cols_m - getRank(false) << std::endl;
        unsigned int dimKernel = cols_m - matCopy.getRank(true);
        //unsigned int dimKernel = cols_m - getRank(false);
        if ( dimKernel==0 )
        {
                #ifdef DEBUG
                std::cerr << "dimKernel ==0" << std::endl;
                #endif
                return new TMatrix<TRing>( cols_m, dimKernel, ring );
        }
        //getRank(true);
        //unsigned int dimKernel = cols_m - getRank(true);
        
        
        TMatrix<TRing>* kernel = new TMatrix<TRing>(cols_m,dimKernel,ring);
        
        // 1. check if x_i is free parameter and store this info in "varFree"   
        vector<bool> varFree = matCopy.getFreeVariables(); 

        //init kernel Matrix
        kernel->fillZero();
        // 2. set current row to last row with non-zero-elements 

        int row = cols_m-dimKernel-1; 
        unsigned int free_varNumber=0; // next free parameter number

        for (int currentVariableNumber = cols_m-1;  currentVariableNumber>=0;   currentVariableNumber--)
        {
                //std::cerr  << "currentVariableNumber " << currentVariableNumber << std::endl;
                //if current variable is free parameter
                if (varFree[currentVariableNumber])
                {
                        // set Kernel_(variableNumber,free_varNumber)=1
                        kernel->setVal(currentVariableNumber,free_varNumber, 1 );
                        free_varNumber++; // (next free parameter number) ++
                }
                //if current variable is NOT free parameter
                else
                {
                        // then set current variable
                        //=M[currentrow][currentVariableNumber+1]*(x_(currentVariableNumber+1))+...+M[currentrow][cols_m-1]*(x_(cols_m-1))
                        //  x_(currentVariableNumber+1) to x_(cols_m-1) already determined and are in Kernel matrix
                        for (unsigned int j=currentVariableNumber+1; j<cols_m; j++)
                        {
                                if ( !matCopy.pcolumn[row][j].isZero() )
                                // wenn Eintrag j<>0, dann muss zu x_currentVariableNumber (Zeile currentVariableNumber im Kernel)
                                // {- div( coeff( Val(row,j), coeff(x_curr) ) * x_j} dazuaddoert werden.
                                // und x_j wurde schon berechnet und steht in der j_ten Zeile der Kernel-Matrix !
                                {
                                        for(unsigned int i=0;i<dimKernel;i++)
                                        {
                                                TNum    kernel_cur_i    = kernel->getVal(currentVariableNumber,i);

                                                TNum    coeff           = ring->multInv(matCopy.pcolumn[row][currentVariableNumber] );

                                                coeff = ring->multiply(coeff, matCopy.pcolumn[row][j]);

                                                coeff = ring->addInv(coeff);

                                                kernel_cur_i = ring->add( kernel_cur_i, 
                                                                          ring->multiply(coeff, kernel->getVal(j, i ) ) 
                                                                        );

                                                kernel->setVal(currentVariableNumber,i, kernel_cur_i);
                                        }
                                }

                        }
                        // change to previous row..
                        row--;
                }
        }
        return kernel;

}


template <class TRing>
unsigned int TMatrix<TRing>::getRank() const
{       
        TMatrix<TRing> matCopy (*this);
        
        return matCopy.getRank(false);
}



template <class TRing>
TMatrix<TRing>*         TMatrix<TRing>::getRowEchelonForm() const
{
        TMatrix<TRing>* matCopy = new TMatrix<TRing>(*this);

        matCopy->getRank(true);
        
        return matCopy;
}


template <class TRing>
TVector<TRing>*         TMatrix<TRing>::solveLGS(TVector<TRing> &rightHandSide)
{

        
        return NULL;
}

template <class TRing>
        unsigned int TMatrix<TRing>::getRank(bool inplace)
        {
                TNum*   dataCopy = NULL;

                if (!inplace)
                {
                        dataCopy = new TNum[rows_m*cols_m];

                        for (size_t i=0; i<rows_m*cols_m; i++)
                                dataCopy[i] =  data[i];
                }
                else
                {
                        dataCopy = data;
                }

                //TNum** pcolumn=new TNum*[rows_m];
                TNum**  pcolumn2 = new TNum*[rows_m];
                TNum*   pfaktor  = new TNum[rows_m];

                for (unsigned int i=0; i<rows_m; i++)
                {
                        pcolumn[i]  = NULL;
                        pcolumn2[i] = NULL;
                }
                TNum*           temprow;
                TNum            pivot;

                for (unsigned int i=0; i<this->rows_m; i++)
                        pcolumn[i] = &dataCopy[ getDataPos(i, 0) ];

                unsigned int row        = 0;
                unsigned int u;
                unsigned int w  = 0;
                unsigned int col        = 0;
                unsigned int temprank= 0;
                
                while (row<this->rows_m && col<this->cols_m)
                {
                        //search for a row with 1 in column 'col' (pivotrow)
                        //for(ring=row;pivotrow<0&&ring<this->rows_m;ring++)
                        for(unsigned int l=row; l<this->rows_m; l++)
                        {
                                if (pcolumn[l][col].isNotZero() )
                                {
                                        // was soll dass denn hier???
                                        pivot=pcolumn[l][col] ;
                                        //pivotrow=l;
                                        //found pivotrow ->,  switch rows_m 'row' and 'pivitrow'
                                        if (l>row)
                                        {
                                                temprow         = pcolumn[l];
                                                pcolumn[l]      = pcolumn[row];
                                                pcolumn[row]    = temprow;
                                        }
                                        temprank++;
                                        
                                        //find relevant rows_m and compute factors
                                        for ( u=row+1, w=0;  u<this->rows_m; u++)
                                        {
                                                if ( pcolumn[u][col].isNotZero() )
                                                {
                                                
                                                        //notice in pcolumn2  where  pcolumn[u][col]!=0
                                                
                                                        //pfaktor[w]   = pcolumn[u][col] / pivot 
                                                        pfaktor[w]   = ring->multiply( ring->multInv(pivot), pcolumn[u][col]  );
                                                        pcolumn2[w++] = pcolumn[u];
                                                }
                                        }

                                        //create zeros:
                                        for(unsigned int v=col;v<this->cols_m;v++)
                                        {
                                                for(u=0;u<w;u++)
                                                {
                                                        pcolumn2[u][v] = ring->add( TNum(pcolumn2[u][v]), 
                                                                                   ring->addInv( ring->multiply( pfaktor[u], 
                                                                                                                          pcolumn[row][v] 
                                                                                                                )
                                                                                                        )
                                                                                );
                                                }
                                        }
                                        //end create zeros
                                        
                                        row++;
                                        /*break;*/goto weiter;
                                }
                        }
                        weiter:
                        col++;
                        
                }
                if (!inplace)
                {
                        for (unsigned int i=0; i<this->rows_m; i++)
                                pcolumn[i] = &data[ getDataPos(i, 0) ];

                        if ( dataCopy!=NULL )
                                { delete[] dataCopy; dataCopy=NULL; }
                }
        
                if (pcolumn2!=NULL) { delete[] pcolumn2; pcolumn2=NULL; }
                if (pfaktor!=NULL) { delete[] pfaktor; pfaktor=NULL; }

                return temprank;
        }



        template <class TRing>
        TMatrix<TRing>* TMatrix<TRing>::getRowBasis() const
        {
                TMatrix<TRing>*   independentRows = getRowEchelonForm(); //Funktion gibt es auch in FFPAck, allerdings weiss ich nicht, wie man diese benutzt.

                independentRows->setName("independent Rows");


                unsigned int independentRowCount        = independentRows->getRank();
        
                TMatrix<TRing>*   rowBasis = new TMatrix<TRing>(independentRowCount, independentRows->getColNum(),independentRows->getRing() );
        
                // da die RowEchelon-form eventuell nicht korrekt sortiert ist, laufe ueber alle Zeilen der ZeilenstufenFormMatrix;
                int currentRowIdx = 0;
                for (unsigned int z=0; z < independentRows->getRowNum() ; z++)
                {
                        TVector<typename TMatrix<TRing>::RingType > * rowAsVector = independentRows->getRowAsVector(z);
                        if ( rowAsVector->isNotZero() )
                        {
                                for (unsigned int colidx=0; colidx < independentRows->getColNum() ;colidx++)
                                {
                                        rowBasis->setVal(currentRowIdx,colidx, rowAsVector->getVal(colidx) );
                                }
                                currentRowIdx++;
                        }
                        delete rowAsVector;
                }
        
                delete independentRows;
        
                return rowBasis;
        }





//---------------------------------------IO------------------------------------------------------
 
        template <class TRing>
        void TMatrix<TRing>::printValue(std::ostream & os) const
        {
                os << "{" << std::endl;
                for (unsigned int i=0; i<rows_m; i++)
                {
                        if (i>0)
                        {
                                os << ", " << std::endl;
                        }
                        os << "{";
                        for (unsigned int j=0; j<cols_m; j++)
                        {
                                if (j>0)
                                {
                                        os << ", ";
                                }
                                os <<  ::std::setw( 2 ) << data[i*cols_m+j];
                                
                        }
                        os << "}" ;
                }               
                os << std::endl << "} ";
        }

        template <class TRing>
        void TMatrix<TRing>::printValueInMacaulayStyle(std::ostream & os) const
        {
                if (rows_m> 0)
                {
                        os <<  " matrix " ;
                        printValue(os);
                }
                else  
                {
                        os <<  " transpose matrix " ;
                        TMatrix<TRing>* tm = getTransposed();
                        tm->printValue(os);
                        delete tm;
                }
        }

        template <class TRing>
        void TMatrix<TRing>::printInMacaulayStyle(std::ostream & os, bool bAssignment) const
        {
                os << std::endl;
                if (bAssignment) 
                        os <<  name << " = " ;
                else 
                        os <<  "--" << name << std::endl;

                printValueInMacaulayStyle(os);
 
                if (bAssignment) 
                        os << ";";
                os << std::endl;
        }


template <class TRing> 
std::ostream &  operator<<(std::ostream & out, const TMatrix<TRing>& z)
{  
        z.printInMacaulayStyle(out);
        return out;
}