CS_ZZ_ICRS.hpp

//#define _DEBUG

#ifndef _H_CS_ZZ_ICRS
#define _H_CS_ZZ_ICRS

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

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

//#_UNSIGNED

#ifdef _DEBUG
#include<iostream>
#endif

template< typename T >
class CS_ZZ_ICRS {

  private:

#ifdef _UNSIGNED

        typedef unsigned long int ULI;
#else

        typedef long int ULI;
#endif

  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 bool compareTriplets( const Triplet< T >& one, const Triplet< T >& two ) {
                if( one.i() < two.i() )
                        return true;
                if ( one.i() > two.i() )
                        return false;
                return one.j() < two.j();
        }

  public:

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

        CS_ZZ_ICRS( CS_ZZ_ICRS< T >& toCopy ) {
                zero_element = toCopy.zero_element;
                nnz = toCopy.nnz;
                noc = toCopy.noc;
                nor = toCopy.nor;
                nzs = new CS_ARRAY< T >( nnz, "Nonzeros array( CS_ZZ_ICRS, copied )" );
                c_ind = new CS_ARRAY< T >( nnz, "Array storing the column index differences (CS_ZZ_ICRS, copied )" );
                r_ind = new CS_ARRAY< T >( nnz, "Array storing the row index differences (CS_ZZ_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_ZZ_ICRS( std::vector< Triplet< T > >& input, const ULI m, const ULI n, const T zero ) {
                zero_element = zero;

                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 );

                //analyse the number of row jumps
                unsigned long int 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.const_access(), "Array of nonzeros (CS_ZZ_ICRS)" );
                c_ind = new CS_ARRAY< ULI >( nnz.const_access() + 1, "Array of column index differences (CS_ZZ_ICRS)" );
                r_ind = new CS_ARRAY< ULI >( r_ind_temp.size() + 1, "Array of row index differences (CS_ZZ_ICRS)" );

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

                //copy row-jump vector
                for( unsigned long int i=0; i<r_ind_temp.size(); i++ )
                        r_ind->access( i ) = r_ind_temp[ i ];

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

                for( ULI i=0; i<nnz.const_access(); ++i ) {
                        const Triplet< T > cur = input[ i ];
                        const unsigned long int currow = cur.i();
                        const unsigned long int curcol = cur.j();
                        if( currow == prev_row ) {
                                c_ind->access( i ) = curcol - prev_col;
#ifdef _DEBUG
                                std::cout << "1. c_ind[ " << i << " ] = " << c_ind->access( i ) << " = " << curcol << " - " << prev_col << ":\t\t " << currow << "," << curcol << std::endl;
#endif
                                nzs->access( i ) = cur.value;
                                prev_col = curcol;
                        } else {
                                //reverse = !reverse;
                                prev_row = currow;
                                ULI stop = i;
                                ULI j = i;
                                
                                while( j < nnz.const_access() && input[ j ].i() == currow )
                                        j++; //go to next row
                                j--; // j now points at last element at new row

                                //process the first new element
                                nzs->access( i ) = input[ j ].value;
                                c_ind->access( i ) = noc.const_access() + ( input[ j ].j() - prev_col );
                                prev_col = input[ j ].j();

#ifdef _DEBUG
                                std::cout << "2. c_ind[ " << i << " ] = " << c_ind->const_access( i ) << " = " << input[ j ].j() << " - " << input[ j-1 ].j() << ":\t\t " << input[ j ].i() << "," << input[ j ].j() << std::endl;
#endif
                                i++; //go to next indices
                                j--;

                                for( ; j>=stop; j-- ) { //reversively add the rest to data structure
                                        nzs->access( i ) = input[ j ].value;
                                        c_ind->access( i ) = prev_col - input[ j ].j();
                                        prev_col = input[ j ].j();
#ifdef _DEBUG
                                        std::cout << "3. c_ind[ " << i << " ] = " << c_ind->const_access( i ) << " = " << input[ j ].j() << " - " << input[ j-1 ].j() << ":\t\t " << input[ j ].i() << "," << input[ j ].j() << std::endl;
#endif
                                        i++; //increment target index
                                }

                                if( i >= nnz.const_access() ) break;

                                //at this point, the reverse row is completely added
                                //and i should be at the next row index.
#ifdef _UNSIGNED
                                c_ind->access( i ) = noc.const_access() + input[ i ].j();
#else
                                c_ind->access( i ) = noc.const_access() + ( prev_col - input[ i ].j() );
#endif

#ifdef _DEBUG
                                std::cout << "4. c_ind[ " << i << " ] = " << c_ind->unrecorded_access( i ) << " = " << noc.unrecorded_access() << " + " << input[ i ].j() << " - " << input[ stop ].j() << ":\t\t " << input[ i ].i() << "," << input[ i ].j() << std::endl;
#endif
                                nzs->access( i ) = input[ i ].value;
                                prev_col = input[ i ].j();
                                prev_row = input[ i ].i();
                        }
                }
        }

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

        T* zax( T* x_orig ) {

                CS_ARRAY< T >* x = new CS_ARRAY< T >( noc.unrecorded_access(), "x-vector in CS_ZZ_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_ZZ_ICRS::zax" );

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

                //free variable in terms of cache because of pointer arithmatic
                ULI row_ind_pointer = 0;//ULI* row_ind_pointer = r_ind;
                ULI index = 0;
                ULI row = 0;
                ULI column = c_ind->const_access( 0 ); //ULI* column = c_ind;
                //ULI* c_ind_p = &( c_ind[ 1 ] );

#ifdef _DEBUG
                std::cout << "nor = " << noc.unrecorded_access() << std::endl; std::cout << "noc = " << noc.unrecorded_access() << std::endl;
#endif

                while( index < nnz.unrecorded_access() ) {
                        while( column < noc.unrecorded_access() ) {
#ifdef _DEBUG
                                std::cout << "ret[ " << row << " ] += nzs[ " << index << " ] * x[ " << column << " ], \t(" << c_ind->unrecorded_access( index + 1 ) << ")" << std::endl;
#endif
                                ret->access( row ) += nzs->const_access( index ) * x->const_access( column );
                                ++index;
                                //*column++ simulation:
                                column += c_ind->const_access( index );
                        }

                        if( index >= nnz.unrecorded_access() ) break;

                        //row += *row_ind_pointer++;
                        //*row++ simulation:
                        row += r_ind->const_access( row_ind_pointer++ );
#ifdef _DEBUG
                        std::cout << column.unrecorded_access() << " > " << noc.unrecorded_access() << std::endl;
#endif
                        column -= noc.unrecorded_access();
#ifdef _UNSIGNED
                        while( column < noc.unrecorded_access() ) {
#else
                        while( column >= 0 ) {
#endif
#ifdef _DEBUG
                                std::cout << "ret[ " << row << " ] += nzs[ " << index << " ] * x[ " << column << " ], \t(-" << c_ind->unrecorded_access( index ) << ")" << std::endl;
#endif
                                ret->access( row ) += nzs->const_access( index ) * x->const_access( column ); //*column++
                                ++index;
#ifdef _UNSIGNED
                                if( c_ind->const_access( index ) <= column ) //this takes waaaay too much time. seems unsolvable in unsigned, switch to signed ints
#endif
                                        column -= c_ind->const_access( index );//*c_ind_p++;
#ifdef _UNSIGNED
                                else {
                                        column = c_ind->const_access( index ) - noc.unrecorded_access();//*c_ind_p++ - noc;
                                        break;
                                }
                        }
#else
                        }

                        column += noc.unrecorded_access();
#endif
                        row += r_ind->const_access( row_ind_pointer++ );//*row_ind_pointer++;
                }

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

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

};

#endif

---