FFLASMatrix.cpp

Go to the documentation of this file.

/*
template <class Field >
template <class TMatrix>
 FFpackMatrix<Field >::FFpackMatrix(const TMatrix &mat,
                                                const int characteristic
                                                //const Field* _F       
                                                                        ):
                                                                                characteristic_m(characteristic),
                                                                                //F( _F ), 
                                                                                F(new Field((long)characteristic)),
                                                                                rows_m( mat.getRowNum() ), 
                                                                                cols_m( mat.getColNum() ),
                                                                                name(  mat.getName() )
{
        checkDimensions();

        data_m          = NULL;
        blocksPerRow    = 0;


        data_m = new typename Field::Element[rows_m*cols_m];
        assert(data_m!=NULL);
        
        int epsExp=0;
        int matrixpos=0;
        
        for (int row=0;row<rows_m;row++)
        {
                for (int col=0;col<cols_m;col++)
                {
                        //data_m[cols_m*row+col]=(mat.getVal(row,col)).getValue(epsExp);
                        data_m[matrixpos] = (mat.getVal(row,col)).getValue(epsExp);
                        matrixpos++;
                }
        }
}*/

template <class Field >
template <class TMatrix>
 FFpackMatrix<Field >::FFpackMatrix(const TMatrix &mat,
                                                int     characteristic,
                                                 bool transpose ):
                                                                                characteristic_m(characteristic),
                                                                                F( new Field((long)characteristic) ),
                                                                                rows_m( mat.getRowNum() ), 
                                                                                cols_m( mat.getColNum() ), 
                                                                                name( mat.getName() )
{
        
        checkDimensions();

        data_m          = NULL;
        blocksPerRow    = 0;


        data_m = new typename Field::Element[rows_m*cols_m];
        assert(data_m!=NULL);

        int epsExp=0;

        int matrixpos=0;
        if (!transpose)
                for (int row=0;row<rows_m;row++)
                {
                        for (int col=0;col<cols_m;col++)
                        {
                                
                                //data_m[cols_m*row+col]=(mat.getVal(row,col)).getValue(epsExp);
                                data_m[matrixpos]=(mat.getVal(row,col)).getValue(epsExp);
                                matrixpos++;
                        }
                }
        else
        {
                rows_m = mat.getColNum();
                cols_m = mat.getRowNum();

                for (int row=0;row<rows_m;row++)
                {
                        for (int col=0;col<cols_m;col++)
                        {
                                //data_m[cols_m*row+col]=(mat.getVal(col,row)).getValue(epsExp);
                                data_m[matrixpos]=(mat.getVal(col,row)).getValue(epsExp);
                                matrixpos++;
                        }
                }
        }
}

template <class Field >
        bool             FFpackMatrix< Field >::isZero()        const
{

        for (int pos=0;pos <rows_m*cols_m;pos++)
        {                       
                if (data_m[pos]!=0)
                        return false;
        }
                
        return true;
}


template <class Field >
 FFpackMatrix< Field >::FFpackMatrix(int        _rows, 
                                                 int    _cols,
                                                 int    characteristic, 
                                                 string _name      ):
                                                                                characteristic_m(characteristic),
                                                                                F(  new Field( (long)characteristic )  ),
                                                                                rows_m(_rows), 
                                                                                cols_m(_cols), 
                                                                                name(_name)
        {
                checkDimensions();

                data_m          = NULL;
                blocksPerRow    = 0;

                data_m = new typename Field::Element[_rows*_cols+1];

                data_m[_rows*_cols]=0;

                assert( data_m!=NULL );

                for (int pos=0;pos <rows_m*cols_m;pos++)
                {                       
                        data_m[pos]=0;
                }

                //memset(data_m ,0x00,rows_m*cols_m*sizeof(typename Field::Element));
        }





template <class Field >
FFpackMatrix< Field >::
FFpackMatrix(typename Field::Element    *_data, 
                                        int             _rows,
                                        int             _cols,
                                        int characteristic, 
                                        string  _name     ):
                                                                        characteristic_m(characteristic),
                                                                        F( new Field( (long)characteristic) ),
                                                                        rows_m( _rows ), 
                                                                        cols_m( _cols ), 
                                                                        data_m( _data ),
                                                                        name( _name )
{
        checkDimensions();

}



template <class Field >
 FFpackMatrix< Field >* FFpackMatrix< Field >::getRightKernel(int &RightKernelDim)
{
        typename Field::Element *NS=NULL;
        size_t NullSpaceDim;

        //size_t lda=rows_m; - falsch!
        size_t lda=cols_m;
        size_t ldn;

        FFPACK::NullSpaceBasis(*F,FFLAS::FflasRight, rows_m, cols_m, data_m, lda, NS, ldn, NullSpaceDim);

        RightKernelDim = NullSpaceDim;

        if (NullSpaceDim==0)
                return NULL;
        return new FFpackMatrix< Field >(NS, cols_m, NullSpaceDim, characteristic_m, string("Kernel("+name+")") );
        
}


/*  leading dimension - Erklärung:

        bei Spaltenorientierter Speicherung von Daten ist LDA die Anzahl der Elemente einer Spalte.( = #Zeilen! )

        Entsprechend ist bei Zeilenorientieruter Speicherung von Daten LDA die Anzahl der Elemente einer Zeile. ( = #Spalten !)

        */
template <class Field >
FFpackMatrix< Field >* FFpackMatrix< Field >::getLeftInverse()
{
        assert(cols_m==cols_m);


        int NullSpaceDim=0;

        size_t matSize = rows_m*cols_m; // A dimension
        size_t lda = cols_m; // A leading dimension, if rows_m==cols_m it does not matter, what ld(A) really is.
        size_t ldx = cols_m; // X leading dimension, if rows_m==cols_m it does not matter, what ld(X) really is.

        
        if (rows_m ==0   )
        {
                std::cerr << " matrix is empty!" << std::endl;
                return NULL;
        }
        #ifdef DEBUG
                std::cerr << "getLeftInverse: matSize " << matSize << std::endl;
                std::cerr << "getLeftInverse:data_m " << *data_m << std::endl;
                std::cerr << "getLeftInverse:lda " << lda<< std::endl;
                std::cerr << "getLeftInverse:cols_m" << cols_m << std::endl;
                std::cerr << "getLeftInverse:rows_m" << rows_m << std::endl;
        #endif
        typename Field::Element * INVDATA= NULL;

        INVDATA = new typename Field::Element [ matSize ];

        for (size_t i=0; i< matSize; i++ ) INVDATA[i]=0;

        #ifdef DEBUG
                std::cerr << "Invert2;" << std::endl;
        #endif
        // Invert2( const Field& F, const size_t M,[in] typename Field::Element * A, const size_t lda,
        //      [out] typename Field::Element * X, const size_t ldx,    [out] int& nullity)
        FFPACK::Invert2(*F, lda, data_m, lda ,  INVDATA, ldx, NullSpaceDim);
        //FFPACK::Invert(*F, lda, data_m, lda, NullSpaceDim);
        #ifdef DEBUG
                printValue(std::cerr);
                
                std::cerr << "Invert2 finished" << std::endl;
        #endif

        if (NullSpaceDim!=0)
        {
                std::cerr << "FFpackMatrix< Field >* FFpackMatrix< Field > getLeftInverse :";
                std::cerr << " matrix is singular!" << std::endl;
                std::cerr << "NullSpaceDim" << NullSpaceDim << std::endl;
                return NULL;
        }
        else
        {
                assert(INVDATA!=0);
        }

        return new FFpackMatrix< Field >(INVDATA, cols_m, rows_m, characteristic_m, string("LeftInverse("+name+")") );
        //return this;
}

template <class Field >
typename Field::Element*        FFpackMatrix< Field >::solveLGS(const typename Field::Element* rightHandSide, int & rankRef)
{
        
        typename Field::Element* solution  = new typename Field::Element[ cols_m ];

        /*
        DOKUMENTATION nicht aktuell!!
        * @brief Rectangular system solver
        const Field& F          * @param Field The computation domain
        const FFLAS_SIDE Side   * @param Side Determine wheter the resolution is left or right looking
        const size_t M,                 * @param M row dimension of A
        const size_t N,         * @param N col dimension of A
        const size_t NRHS,      * @param NRHS number of columns (if Side = FflasLeft) or row (if Side = FflasRight) of the matrices X and B
        Field::Element *A,      * @param A input matrix
        const size_t lda,               * @param lda leading dimension of A
        Field::Element *X,
        const size_t ldx,
        const Field::Element *B, * @param B Right/Left hand side matrix. Initially contains B, finally contains the solution X.
        const size_t ldb,               * @param ldb leading dimension of B
        int * info                      * @info Success of the computation: 0 if successfull, >0 if system is inconsistent
                                        * @return the rank of the system
        
        ld: leading dimension:
                Suppose that you have a matrix A of size 100x100 which is stored in an array 100x100. 
                In this case LDA is the same as N. Now suppose that you want to work only on the submatrix A(91:100 , 1:100); 
                in this case the number of rows is 10 but LDA=100. Assuming the fortran column-major ordering (which is the case in LAPACK), 
                the LDA is used to define the distance in memory between elements of two consecutive columns which have the same row index. 
                If you call B = A(91:100 , 1:100) then B(1,1) and B(1,2) are 100 memory locations far from each other. 
        */

        /*column-major ordering
        bei Spalten-Orientierter Speicherung von Daten ist LDA die Anzahl der Elemente einer Spalte.(= Anzahl Zeilen!)

        Entsprechend ist bei Zeilen-orientierter Speicherung von Daten LDA die Anzahl der Elemente einer Zeile. (= Anzahl Spalten!)

        */

        size_t NRHS = 1; // ist unklar , was NHRS  ist , angeblich die Anzahl der Spalten von B. 
        size_t ldx = 1;
        size_t ldb = 1;
        size_t lda = cols_m;

        int info = 1;
        
        rankRef = FFPACK::fgesv(        *F, 
                                                FFPACK::FflasLeft, 
                                                rows_m, 
                                                cols_m, 
                                                NRHS, 
                                                data_m, 
                                                lda, 
                                                solution, 
                                                ldx, 
                                                rightHandSide, 
                                                ldb, 
                                                &info           );

        if (info != 0)
        {
                //      cerr << "LGSSolve: no solutions" << endl;
                return NULL;
        }
                
        return solution;
}


template <class Field >
  unsigned int FFpackMatrix< Field >::getRank()
{
        return (unsigned int)FFPACK::Rank (*F, rows_m, cols_m, data_m, rows_m);
}



template <class Field >
 typename Field::Element FFpackMatrix< Field >::getVal(int row, int col) const
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif 

        return data_m[cols_m*row+col];
}


template <class Field >
 void   FFpackMatrix<Field >::setVal(int row, int col,typename Field::Element elem)
{
        #ifdef SAFE
                checkBounds(row,col);
        #endif 

        data_m[cols_m*row+col]=elem;
}

template <class Field >
 void  FFpackMatrix<Field >::checkDimensions() const
{
        if (rows_m<=0   || cols_m<=0 )
        {
                std::cerr << "rows_m " << rows_m << std::endl;
                std::cerr << "cols_m " << cols_m << std::endl;
        }
        assert( rows_m>=0 );
        assert( cols_m>=0 );
}

template <class Field >
 void  FFpackMatrix<Field >::checkBounds(int row,int col) const
{
        if (row<0  || row>=rows_m || col<0  || col>=cols_m)
        {
                std::cerr << "row " << row << " rows_m " << rows_m << std::endl;
                std::cerr << "col " << col << " cols_m " << cols_m << std::endl;
                assert(false);
        }
}

template <class Field>
bool FFpackMatrix<Field>::operator==(const FFpackMatrix& 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 Field >
 FFpackMatrix<Field >::~FFpackMatrix()
        {
                if (data_m!=NULL)
                delete[] data_m;

                if (F!=NULL)  delete F;
                
        };


template <class Field >
void  FFpackMatrix<Field >::printValue(std::ostream & os) const
{
        
        os << std::endl << "{" ;
        for ( int i=0; i<rows_m; i++)
        {
                if (i>0)
                {
                        os << ", " << std::endl;
                }
                os << "{";
                
                for ( int j=0; j<cols_m; j++)
                {
                        if (j>0)
                        {
                                os << ", ";
                        }
                        os << ::std::setw( 2 ) << getVal(i,j);
                        
                }
                os << "}" ;
        }               
        os << std::endl << "} ";
}