CS_ICRS.hpp

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

#include "CS_ELEMENT.hpp"
#include "CS_ARRAY.hpp"

//#define _DEBUG

#ifndef _H_CS_ICRS
#define _H_CS_ICRS

#ifdef _DEBUG
#include<iostream>
#endif

template< typename T >
class CS_ICRS {

  private:

        typedef unsigned long int ULI;

  protected:

        CS_ELEMENT< ULI > nnz;

        CS_ELEMENT< ULI > nor;

        CS_ELEMENT< ULI > noc;

        CS_ARRAY< T >* nzs;

        CS_ARRAY< ULI >* c_ind;

        CS_ARRAY< ULI >* r_ind;

        CS_ELEMENT< T > zero_element;

        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.i() < two.i() )
                        return -1;
                if ( one.i() > two.i() )
                        return 1;
                if ( one.j() < two.j() )
                        return -1;
                if ( one.j() > two.j() )
                        return 1;
                return 0;
        }

  public:

        CS_ICRS( const ULI number_of_nonzeros, const ULI number_of_rows, const ULI number_of_cols, const T zero ):
                nnz( number_of_nonzeros ), noc( number_of_cols ), nor( number_of_rows ), zero_element( zero ) {
                nzs = new CS_ARRAY< T>( nnz.const_access(), "Nonzeros array (ICRS)" );
                c_ind = new CS_ARRAY< ULI >( nnz.const_access(), "Array storing the column index differences (ICRS)" );
                r_ind = new CS_ARRAY< ULI >( nnz.const_access(), "Array storing the row index differences (ICRS)" );
        }

        CS_ICRS( CS_ICRS< T >& toCopy ) {
                zero_element = toCopy.zero_element;
                nnz = toCopy.nnz;
                noc = toCopy.noc;
                nor = toCopy.nor;
                nzs = new CS_ARRAY< T >( nnz.const_access(), "Nonzeros array( CS_ICRS, copied )" );
                c_ind = new CS_ARRAY< ULI >( nnz.const_access(), "Array storing the column indices (CS_ICRS, copied)" );
                r_ind = new CS_ARRAY< ULI >( nnz.const_access(), "Array storing the row index indices (CS_ICRS, copied)" );
                for( CS_ELEMENT< ULI > i=0; i<nnz; i = i + 1 ) {
                        nzs->access( i.const_access() ) = toCopy.nzs->const_access( i.const_access() );
                        c_ind->access( i.const_access() ) = toCopy.c_ind->const_access( i.const_access() );
                        r_ind->access( i.const_access() ) = toCopy.r_ind->const_access( i.const_access() );
                }
        }

        CS_ICRS( std::vector< Triplet< T > >& input, const ULI m, const ULI n, T& zero ) {
                zero_element = zero;
                //find nnz
                nnz = input.size();

                nor = m;
                noc = n;
        
                typename std::vector< Triplet< T > >::iterator in_it;

                //WARNING: noc*nor typically overflows on 32 bits!
                //         immediately start recording differences
                //         instead of trying to directly calculate
                //         the index as i*noc+j.
        
                //sort( input.begin(), input.end(), compareTriplets );
                qsort( &input[ 0 ], input.size(), sizeof( Triplet< T > ), &compareTriplets );

                //filter out zeros and count the number of row jumps
                ULI prev_row = 0;
                std::vector< ULI > r_ind_temp;
                typename std::vector< Triplet< T > >::iterator it = input.begin();
                for( ; it!=input.end(); it++ ) {
                        if( (*it).i() > prev_row ) { r_ind_temp.push_back( (*it).i() - prev_row ); prev_row = (*it).i(); }
                }
                nnz = input.size();

                //allocate arrays
                nzs = new CS_ARRAY< T >( nnz.unrecorded_access(), "Array of nonzeros (CS_ICRS)" );
                c_ind = new CS_ARRAY< ULI >( nnz.unrecorded_access() + 1, "Array of column index indices (CS_ICRS)" );
                r_ind = new CS_ARRAY< ULI >( r_ind_temp.size() + 1, "Array of row index indices (CS_ICRS)" );

                //set last entry
                c_ind->unrecorded_access( nnz.unrecorded_access() ) = noc.unrecorded_access();
                //r_ind does not have to be set; although the last element is read, it is actually never used.

                //copy row-jump vector
                for( ULI i=0; i<r_ind_temp.size(); i++ )
                        r_ind->unrecorded_access( i ) = r_ind_temp[ i ];

                //make ICRS
                ULI prev_col = 0;
                prev_row = 0;

                for( ULI i=0; i<nnz.unrecorded_access(); ++i ) {
                        const Triplet< T > cur = input[ i ];
                        const ULI currow = cur.i();
                        const ULI curcol = cur.j();
                        if( currow == prev_row ) {
                                c_ind->unrecorded_access( i ) = curcol - prev_col;
                        } else {
                                c_ind->unrecorded_access( i ) = noc.unrecorded_access() + ( curcol - prev_col );
                                prev_row = currow;
                        }
                        nzs->unrecorded_access( i ) = cur.value;
                        prev_col = curcol;

#ifdef _DEBUG
                        std::cout << currow << "," << curcol << "(" << cur.value << ") maps to " << c_ind[ i ] << std::endl;
#endif
                }
        }

        /*CS_ICRS( std::vector< Triplet< T > >& input, T& zero ) {
                zero_element = zero;
                //find nnz
                nnz = input.size();
        
                //find number of cols
                noc = nor = static_cast< ULI >( 0 );
                typename std::vector< Triplet< T > >::iterator in_it;
                in_it = input.begin();
                for( ; in_it!=input.end(); in_it++ ) {
                        const CS_ELEMENT< ULI > currow = ( *in_it ).i();
                        const CS_ELEMENT< ULI > curcol = ( *in_it ).j();
                        if ( currow > nor ) nor = currow;
                        if ( curcol > noc ) noc = curcol;
                }

                ++nor; ++noc;

#ifdef _DEBUG
                std::cout << "Max col index found: " << noc << std::endl;
                std::cout << "Max row index found: " << nor << std::endl;
#endif
        
                //WARNING: noc*nor typically overflows on 32 bits!
                //         immediately start recording differences
                //         instead of trying to directly calculate
                //         the index as i*noc+j.
        
                sort( input.begin(), input.end(), compareTriplets );

                //allocate arrays
                nzs = new CS_ARRAY< T >( nnz.const_access(), "Array of nonzeros (CS_ICRS)" );
                ind = new CS_ARRAY< ULI >( nnz.const_access(), "Array of index differences (CS_ICRS)" );
                ULI prev_col = 0;
                ULI prev_row = 0;

                //make CRS
                for( CS_ELEMENT< ULI > i=0; i<nnz; ++i ) {
                        const Triplet< T > cur = input[ i.const_access() ];
                        if( cur.value == zero_element.const_access() ) { nnz.access() -= 1; continue; }

                        const ULI curdiff = ( cur.i() - prev_row ) * ( noc.const_access() - 1 ) + cur.j() - prev_col;
                        ind->access( i.const_access() ) = curdiff;
                        nzs->access( i.const_access() ) = cur.value;
                        prev_row = cur.i();
                        prev_col = cur.j();

                        //the below overflows like crazy
                        //nzs[ i ] = cur.value;
                        //ind[ i ] = cur.i() * noc + cur.j();
#ifdef _DEBUG
                        std::cout << cur.i() << "," << cur.j() << " maps to " << ind[ i.const_access() ] << " / " << ( nor*noc ) << std::endl;
#endif
                }
        }*/

        /*T& random_access( const ULI& i, const ULI& j ) { TODO
        }*/

        T* zax( T* x_orig ) {

                CS_ARRAY< T >* x = new CS_ARRAY< T >( noc.unrecorded_access(), "x-vector in CS_ICRS::zax" );
                x->unrecordedCopyFrom( x_orig );

                T* true_ret = new T[ nor.unrecorded_access() ];
                CS_ARRAY< T >* ret = new CS_ARRAY< T >( nor.unrecorded_access(), "z-vector in CS_ICRS::zax" );

                for( ULI i=0; i < nor.unrecorded_access(); ++i )
                        ret->access( i ) = zero_element.unrecorded_access();

                ULI row_ind_pointer = 0; //ULI* row_ind_pointer = &(r_ind[0]) ( = r_ind, when using raw array types)
                ULI index = 0;
                ULI row = 0;
                ULI column = c_ind->const_access( 0 ); //ULI* column = c_ind;

                while( index < nnz.unrecorded_access() ) {
                        while( column < noc.unrecorded_access() ) {
                                ret->access( row ) += nzs->const_access( index ) * x->const_access( column ); //*column++
                                index += 1;
                                column += c_ind->const_access( index );//colum += *c_ind_p++
                        }
                        row += r_ind->const_access( row_ind_pointer ); // row += *row_ind_pointer++; i.e., no cache effects due to indirect access
                        row_ind_pointer++;
                        column -= noc.unrecorded_access();
                }

                ret->unrecordedCopyTo( true_ret );
                delete ret;
                delete x;
                return true_ret;
        }

        T* zaxtrace( T* x_orig ) {

                CS_ARRAY< T >* x = new CS_ARRAY< T >( noc.unrecorded_access(), "x-vector in CS_ICRS::zax" );
                x->unrecordedCopyFrom( x_orig );

                T* true_ret = new T[ nor.unrecorded_access() ];
                CS_ARRAY< T >* ret = new CS_ARRAY< T >( nor.unrecorded_access(), "z-vector in CS_ICRS::zax" );

                for( ULI i=0; i < nor.unrecorded_access(); ++i )
                        ret->access( i ) = zero_element.unrecorded_access();

                ULI row_ind_pointer = 0; //ULI* row_ind_pointer = &(r_ind[0]) ( = r_ind, when using raw array types)
                ULI index = 0;
                ULI row = 0;
                ULI column = c_ind->const_access( 0 ); //ULI* column = c_ind;

                while( index < nnz.unrecorded_access() ) {
                        while( column < noc.unrecorded_access() ) {
                                std::cout << "TRACE:"<<row<<","<<column<<std::endl;
                                ret->access( row ) += nzs->const_access( index ) * x->const_access( column ); //*column++
                                index += 1;
                                column += c_ind->const_access( index );//colum += *c_ind_p++
                        }
                        row += r_ind->const_access( row_ind_pointer ); // row += *row_ind_pointer++; i.e., no cache effects due to indirect access
                        row_ind_pointer++;
                        column -= noc.unrecorded_access();
                }

                ret->unrecordedCopyTo( true_ret );
                delete ret;
                delete x;
                return true_ret;
        }

        ~CS_ICRS() {
                delete nzs;
                delete c_ind;
                delete r_ind;
        }

};

#endif

---