Hilbert.hpp

/*
 * Copyright (c) 2007-2014, A. N. Yzelman,   Utrecht University 2007-2011;
 *                                                    KU Leuven 2011-2014.
 *                          R. H. Bisseling, Utrecht University 2007-2014.
 * 
 * This file is part of the Sparse Library.
 * 
 * This library was developed under supervision of Prof. dr. Rob H. Bisseling at
 * Utrecht University, from 2007 until 2011. From 2011-2014, development continued 
 * at KU Leuven, where Prof. dr. Dirk Roose contributed significantly to the ideas 
 * behind the newer parts of the library code.
 * 
 *     The Sparse Library is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by the
 *     Free Software Foundation, either version 3 of the License, or (at your
 *     option) any later version.
 * 
 *     The Sparse Library 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 General Public License
 *     for more details.
 * 
 *     You should have received a copy of the GNU General Public License along
 *     with the Sparse Library. If not, see <http://www.gnu.org/licenses/>.
 */


/*
 * File created by:
 *     A. N. Yzelman, Dept. of Mathematics, Utrecht University, 2010.
 */


#include <vector>
#include "SparseMatrix.hpp"
#include "HilbertTriplet.hpp"
#include "HilbertTripletCompare.hpp"
#include "Triplet.hpp"

#ifndef _H_HILBERT
#define _H_HILBERT

#define _HILBERT_COMPRESSED_BICRS

#ifdef _DEBUG
        #include <iostream>
#endif

#ifdef _HILBERT_COMPRESSED_BICRS
        #include "CBICRS.hpp"
#else
        #include "BICRS.hpp"
#endif

template< typename T >
class Hilbert: public SparseMatrix< T, ULI > {

   private:

   protected:

        ULI minexp;

        std::vector< HilbertTriplet< T > > ds;

        Matrix< T > *ads;

        static Matrix< T >* getDataStructure( std::vector< Triplet< T > > &tds, const ULI m, const ULI n, T zero ) {
#ifdef _HILBERT_COMPRESSED_BICRS
                return CBICRS_factory< T >::getCBICRS( tds, m, n, zero );
#else
                return new BICRS< T >( tds, m, n, zero );
#endif
        }

   public:

        virtual ~Hilbert() {
                if( ads != NULL )
                        delete ads;
        }

        Hilbert() {
                ads = NULL;
        }

        Hilbert( std::string file, T zero = 0 ) {
                ads = NULL;
                this->loadFromFile( file, zero );
        }
        
        Hilbert( std::vector< Triplet< T > >& input, ULI m, ULI n, T zero ) {
                ads = NULL;
                load( input, m, n, zero );
        }

        virtual void load( std::vector< Triplet< T > >& input, const ULI m, const ULI n, const T zero ) {
                this->zero_element = 0;
                this->nor = m;
                this->noc = n;
                this->nnz = input.size();
                for( ULI i=0; i<this->nnz; i++ ) {
                        ds.push_back( HilbertTriplet< T >( input[ i ].i(), input[ i ].j(), input[ i ].value ) );
                        ds[ i ].calculateHilbertCoordinate();
                }
                HilbertTripletCompare< T > compare;
                std::sort( ds.begin(), ds.end(), compare );
#ifdef _DEBUG
                for( ULI i=0; i<this->nnz; i++ )
                        std::cout << ds[i].getMostSignificantHilbertBits() << " " << ds[i].getLeastSignificantHilbertBits() << std::endl;
#endif
                //create temporary structure
                std::vector< Triplet< T > > tds;
                typename std::vector< HilbertTriplet< T > >::iterator it = ds.begin();
                for( ; it!=ds.end(); ++it )
                        tds.push_back( Triplet< T >( it->i(), it->j(), it->value ) );
                //create actual structure
                ads = getDataStructure( tds, m, n, zero );
        }

        virtual void getFirstIndexPair( ULI &row, ULI &col ) {
                row = ds[ 0 ].i();
                col = ds[ 0 ].j();
        }

        virtual void zxa( const T* x, T* z ) {
                ads->zxa( x, z );
        }

        virtual void zax( const T* x, T* z ) {
                ads->zax( x, z );
        }

        virtual size_t bytesUsed() {
                return ads->bytesUsed();
        }

        void saveBinary( const std::string fn ) {
                HilbertTriplet< T >::save( fn, ds, this->nor, this->noc );
        }

        void loadBinary( const std::string fn ) {
                std::cerr << "Warning: assuming binary file was saved by a HTS scheme, i.e., that it is pre-ordered using Hilbert coordinates." << std::endl;
                std::vector< Triplet< T > > tds = Triplet< T >::load( fn, this->nor, this->noc );
                ads = getDataStructure( tds, this->nor, this->noc, 0 );
        }

};

#endif