CRS.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, 2008.
 */


#include "Triplet.hpp"
#include "SparseMatrix.hpp"
#include <vector>
#include <assert.h>

//#define _DEBUG

#ifndef _H_CRS
#define _H_CRS

#ifdef _DEBUG
#include<iostream>
#endif

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

  private:

  protected:

        ULI* row_start;

        T* ds;

        ULI* col_ind;

        static int compareTriplets( const void * left, const void * right ) {
                const Triplet< T > one = **( (Triplet< T > **)left );
                const Triplet< T > two = **( (Triplet< T > **)right );
                if ( one.j() < two.j() )
                        return -1;
                if ( one.j() > two.j() )
                        return 1;
                return 0;
        }     

        bool find( const ULI col_index, const ULI search_start, const ULI search_end, ULI &ret ) {
                for( ULI i=search_start; i<search_end; i++ )
                        if( col_ind[ i ] == col_index ) {
                                ret = i;
                                return true;
                        }
                return false;
        }

  public:

        CRS() {}

        CRS( std::string file, T zero = 0 ) {
                this->loadFromFile( file, zero );
        }
        
        CRS( const ULI number_of_nonzeros, const ULI number_of_rows, const ULI number_of_cols, T zero ) {
                this->nnz = number_of_nonzeros;
                this->nor = number_of_rows;
                this->noc = number_of_cols;
                this->zero_element = zero;
                row_start = new ULI[ this->nor + 1 ];
                ds = new T[ this->nnz ];
                col_ind = new ULI[ this->nnz ];
        }

        CRS( CRS< T >& toCopy ) {
                this->zero_element = toCopy.zero_element;
                this->nnz = toCopy.nnz;
                this->nor = toCopy.nor;
                row_start = new ULI[ this->nor + 1 ];
                ds = new T[ this->nnz ];
                col_ind = new ULI[ this->nnz ];
                for( ULI i=0; i<this->nnz; i++ ) {
                        ds[ i ] = toCopy.ds[ i ];
                        col_ind[ i ] = toCopy.col_ind[ i ];
                }
                for( ULI i=0; i<this->nor; i++ )
                        row_start[ i ] = toCopy.row_start[ i ];
        }
        
        CRS( std::vector< Triplet< T > > input, ULI m, ULI n, T zero ) {
                load( input, m, n, zero );
        }

        virtual void load( std::vector< Triplet< T > >& input, ULI m, ULI n, T zero ) {
                std::cout << "\tLoading in a vector of " << input.size() << " triplets into CRS..." << std::endl;

                this->zero_element = zero;
                //find nnz
                this->nnz = input.size();
        
                this->nor = m;
                this->noc = n;

                //build better datastructure
                std::vector< std::vector< Triplet< T >* > > ds( this->nor );
                
                //move input there
                typename std::vector< Triplet< T > >::iterator in_it = input.begin();
                for( ; in_it != input.end(); ++in_it ) {
                        //Triplet< T >* cur = &(*in_it);
                        const ULI currow = in_it->i();
                        const T value = in_it->value;
                        if( value == this->zero_element ) { 
                                this->nnz--;
                                continue;
                        }
                        ds.at( currow ).push_back( &(*in_it) );
                }

                //allocate arrays
                row_start = new ULI[ this->nor + 1 ];
                this->ds = new T[ this->nnz ];
                col_ind = new ULI[ this->nnz ];

                //make CRS
                ULI index = 0;
                for( ULI currow = 0; currow < this->nor; currow++ ) {
                        row_start[ currow ] = index;
                        if( ds.at( currow ).size() == 0 ) continue;
                        qsort( &( ds.at( currow )[ 0 ] ), ds.at( currow ).size(), sizeof( Triplet< T >* ), &compareTriplets );
                        typename std::vector< Triplet< T >* >::iterator row_it = ds.at( currow ).begin();
                        for( ; row_it!=ds.at( currow ).end(); row_it++ ) {
                                const Triplet< T > cur = *(*row_it);
                                this->ds[ index ] = cur.value;
                                col_ind[ index ] = cur.j();
                                index++;
                        }
                }
                row_start[ this->nor ] = this->nnz;
                std::cout << "\t" << index << " nonzeroes loaded into CRS structure.\n" << std::endl;
                assert( index == this->nnz );
        }

        T& random_access( ULI i, ULI j ) {
                ULI found_index;
                if ( find( j, row_start[ i ], row_start[ i+1 ], found_index ) ) {
#ifdef _DEBUG
                        std::cout << "Searched col_ind between " << row_start[ i ] << " and " << row_start[ i + 1 ] << ", found: " << std::endl;
                        std::cout << "Element (" << i << "," << j << ") found on index " << found_index << ", returning " << ds[ found_index ] << std::endl;
#endif
                        return ds[ found_index ];
                } else
                        return this->zero_element;
        }

        virtual void getFirstIndexPair( ULI &row, ULI &col ) {
                row = this->row_start[ 0 ];
                col = this->col_ind[ 0 ];
        }

        virtual void zxa( const T*__restrict__ x, T*__restrict__ z ) {
                T*__restrict__ ds_p = ds;
                ULI*__restrict__ col_ind_p = col_ind;
                ULI row, index;
                for( row = 0; row < this->nor; row++, x++ ) {
                        for( index = row_start[ row ]; index < row_start[ row + 1 ]; index++ )
                                z[ *col_ind_p++ ] += *ds_p++ * *x;
                }
        }

        virtual void zax( const T*__restrict__ x, T*__restrict__ z ) {
                T*__restrict__ ds_p = ds;
                ULI*__restrict__ col_ind_p = col_ind;
                ULI row, index;
                for( row = 0; row < this->nor; row++, z++ ) {
                        for( index = row_start[ row ]; index < row_start[ row + 1 ]; index++ )
                                *z += *ds_p++ * x[ *col_ind_p++ ];
                }
        }

        template< size_t k >
        void ZaX( const T*__restrict__ const *__restrict__ const X, T *__restrict__ const *__restrict__ const Z ) {
                T*__restrict__ ds_p = ds;
                ULI*__restrict__ col_ind_p = col_ind;
                ULI row, index;
                //local buffer of pointers
                const T *__restrict__ pDataX[ k ];
                      T *__restrict__ pDataZ[ k ];
                //fill local pointer buffer
                for( size_t s = 0; s < k; ++s ) {
                        pDataX[s] = X[s];
                        pDataZ[s] = Z[s];
                }
                for( row = 0; row < this->nor; ++row ) {
                        for( index = row_start[ row ]; index < row_start[ row + 1 ]; ++index ) {
                                for( size_t s = 0; s < k; ++s ) {
                                        assert( pDataZ[s] >= Z[s] );
                                        assert( pDataZ[s] < Z[s]+this->nor );
                                        assert( *col_ind_p < this->noc );
                                        assert( pDataX[s] == X[s] );
                                        *(pDataZ[s]) += *ds_p * pDataX[s][ *col_ind_p ];
                                }
                                ++col_ind_p;
                                ++ds_p;
                        }
                        for( size_t s = 0; s < k; ++s ) {
                                ++(pDataZ[s]);
                        }
                }
        }

        template< size_t k >
        void ZXa( const T *__restrict__ const *__restrict__ const X, T *__restrict__ const *__restrict__ const Z ) {
                T*__restrict__ ds_p = ds;
                ULI*__restrict__ col_ind_p = col_ind;
                ULI row, index;
                const T *__restrict__ pDataX[ k ];
                      T *__restrict__ pDataZ[ k ];
                for( size_t s = 0; s < k; ++s ) {
                        pDataX[s] = X[s];
                        pDataZ[s] = Z[s];
                }
                for( row = 0; row < this->nor; ++row ) {
                        for( index = row_start[ row ]; index < row_start[ row + 1 ]; ++index ) {
                                for( size_t s = 0; s < k; ++s ) {
                                        pDataZ[s][ *col_ind_p ] += *ds_p * *(pDataX[s]);
                                }
                                ++ds_p;
                                ++col_ind_p;
                        }
                        for( size_t s = 0; s < k; ++s ) {
                                ++(pDataX[s]);
                        }
                }
        }

        ULI* rowJump() { return row_start; }

        ULI* columnIndices() { return col_ind; }

        T* values() { return ds; } 

        virtual ~CRS() {
                delete [] row_start;
                delete [] ds;
                delete [] col_ind;
        }

        virtual size_t bytesUsed() {
                return sizeof( ULI ) * ( this->nor + this->nnz ) + sizeof( T ) * this->nnz;
        }

};
 
#undef _DEBUG
#endif