CS_CRS.hpp

#include "Triplet.hpp"
#include "CRS.hpp"
#include "CS_ELEMENT.hpp"
#include "CS_ARRAY.hpp"

#include <cstdlib>
#include <vector>

//#define _DEBUG

#ifndef _H_CS_CRS
#define _H_CS_CRS

#ifdef _DEBUG
#include<iostream>
#endif

template< typename T >
class CS_CRS {

  private:

        typedef unsigned long int ULI;

  protected:

        CS_ELEMENT< ULI >* nnz;

        CS_ELEMENT< ULI >* nor;

        CS_ELEMENT< ULI >* noc;

        CS_ARRAY< ULI >* row_start;

        CS_ARRAY< T >* nzs;

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

        bool find( ULI col_index, ULI search_start, 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:

        CS_CRS( const ULI number_of_nonzeros, const ULI number_of_rows, const ULI number_of_cols, const T zero ):
                nnz( number_of_nonzeros ), nor( number_of_rows ) {
                nnz = new CS_ELEMENT< ULI >( number_of_nonzeros );
                nor = new CS_ELEMENT< ULI >( number_of_rows );
                noc = new CS_ELEMENT< ULI >( number_of_cols );
                zero_element = new CS_ELEMENT< T >( zero );
                row_start = new CS_ARRAY< ULI >( nor + 1, "Row start vector (CS_CRS.h)" );
                nzs = new CS_ARRAY< T>( nnz, "Non-zero storage vector (CS_CRS.h" );
                col_ind = new CS_ARRAY< ULI >( nnz, "Columm index vector (CS_CRS.h)" );
        }

        CS_CRS( CRS< T >& toCopy ) {
                zero_element = new CS_ELEMENT< T >( toCopy.zero_element );      
                nnz = new CS_ELEMENT< ULI >( toCopy.nnz );
                nor = new CS_ELEMENT< ULI >( toCopy.nor );
                noc = new CS_ELEMENT< ULI >( toCopy.noc );
                row_start = new CS_ARRAY< ULI >( nor + 1, "Row start vector (CS_CRS.h)" );
                nzs = new CS_ARRAY< T>( nnz, "Non-zero storage vector (CS_CRS.h" );
                col_ind = new CS_ARRAY< ULI >( nnz, "Columm index vector (CS_CRS.h)" );
                for( ULI i=0; i<(*nnz).const_access(); i++ ) {
                        nzs[ i ] = toCopy.nzs->const_access( i );
                        col_ind[ i ] = toCopy.col_ind->const_access( i );
                }
                for( ULI i=0; i<(*nor).const_access(); i++ )
                        row_start[ i ] = toCopy.row_start->const_access( i );
        }

        CS_CRS( CS_CRS< T >& toCopy ) {
                zero_element = toCopy.zero_element;
                nnz = toCopy.nnz;
                nor = toCopy.nor;
                noc = toCopy.noc;
                row_start = new CS_ARRAY< ULI >( nor + 1, "Row start vector (CS_CRS.h)" );
                nzs = new CS_ARRAY< T>( nnz, "Non-zero storage vector (CS_CRS.h" );
                col_ind = new CS_ARRAY< ULI >( nnz, "Columm index vector (CS_CRS.h)" );
                for( ULI i=0; i<(*nnz).const_access(); i++ ) {
                        nzs[ i ] = toCopy.nzs->const_access( i );
                        col_ind[ i ] = toCopy.col_ind->const_access( i );
                }
                for( ULI i=0; i<(*nor).const_access(); i++ )
                        row_start[ i ] = toCopy.row_start->const_access( i );
        }

        CS_CRS( std::vector< Triplet< T > > input, unsigned long int currow, unsigned long int curcol, T& zero ) {
                zero_element = new CS_ELEMENT< T >( zero );
                nnz = new CS_ELEMENT< ULI >( input.size() );
                nor = new CS_ELEMENT< ULI >( static_cast< ULI >( currow ) );
                noc = new CS_ELEMENT< ULI >( static_cast< ULI >( curcol ) );
        
                //build better datastructure
                std::vector< std::vector< Triplet< T >* > > ds( nor->const_access() );
                
                //move input there
                typename std::vector< Triplet< T > >::iterator in_it;
                in_it = input.begin();
                for( ; in_it != input.end(); in_it++ ) {
                        Triplet< T >* cur = &(*in_it);
                        const ULI currow = cur->i();
                        const T value = cur->value;
                        if( value == zero_element->const_access() ) { nnz->access() -= 1; continue; }
                        ds.at( currow ).push_back( cur );
                }

                //allocate arrays
                row_start = new CS_ARRAY< ULI >( nor->const_access() + 1, "Row start vector (CS_CRS.h)" );
                nzs = new CS_ARRAY< T>( nnz->const_access(), "Non-zero storage vector (CS_CRS.h" );
                col_ind = new CS_ARRAY< ULI >( nnz->const_access(), "Columm index vector (CS_CRS.h)" );

                ULI index = 0;
                for( ULI currow = 0; currow < nor->const_access(); currow++ ) {
                        row_start->access( 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);
                                nzs->access( index ) = cur.value;
                                col_ind->access( index ) = cur.j();
                                index++;
                        }
                }

                row_start->access( nor->const_access() ) = nnz->const_access();
        }

        CS_CRS( std::vector< Triplet< T > > input, T& zero ) {
                zero_element = new CS_ELEMENT< T >( zero );
                //find nnz
                nnz = new CS_ELEMENT< ULI >( input.size() );
        
                //find number of rows
                nor = new CS_ELEMENT< ULI >( static_cast< ULI >( 0 ) );
                noc = new CS_ELEMENT< ULI >( static_cast< ULI >( 0 ) );
                typename std::vector< Triplet< T > >::iterator in_it;
                in_it = input.begin();
                for( ; in_it!=input.end(); in_it++ ) {
                        const ULI currow = ( *in_it ).i();
                        const ULI curcol = ( *in_it ).j();
                        if ( currow > nor->const_access() ) nor->access() = currow;
                        if ( curcol > noc->const_access() ) noc->access() = curcol;
                }
#ifdef _DEBUG
                std::cout << "Max row index found: " << nor->const_access() << std::endl;
#endif
                nor->access() += 1;
                noc->access() += 1;

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

                //allocate arrays
                row_start = new CS_ARRAY< ULI >( nor->const_access() + 1, "Row start vector (CS_CRS.h)" );
                nzs = new CS_ARRAY< T>( nnz->const_access(), "Non-zero storage vector (CS_CRS.h" );
                col_ind = new CS_ARRAY< ULI >( nnz->const_access(), "Columm index vector (CS_CRS.h)" );

                //make CRS
                ULI index = 0;
                for( ULI currow = 0; currow < nor->const_access(); currow++ ) {
#ifdef _DEBUG
                        std::cout << "row_start[ " << currow << " ]=" << index << std::endl;
#endif
                        row_start->access( currow ) = index;
                        if( ds.at( currow ).size() == 0 ) continue;
                        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);
                                nzs->access( index ) = cur.value;
                                col_ind->access( index ) = cur.j();
                                index++;
                        }
                }
                row_start->access( nor->const_access() ) = nnz->const_access();
        }

        T& random_access( ULI i, ULI j ) {
                const ULI found_index;
                if ( find( j, row_start[ i ], row_start[ i+1 ], found_index ) )
                        return (*nzs)[ found_index ];
                else
                        return zero_element;
        }

        T* zax( T* x_orig ) {

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

//              ULI index = 0;
//              ULI row = 0;
//              ULI row_i = row_start->const_access( 1 );
#ifdef _DEBUG
                std::cout << "Debug: nor=" << nor.unrecorded_access() << std::endl;
#endif
                CS_ARRAY< T >* ret = new CS_ARRAY< T >( nor->unrecorded_access(), "z-vector in CS_CRS::zax" );
                for( ULI i=0; i < nor->unrecorded_access(); ++i )
                        ret->access( i ) = zero_element->unrecorded_access();
//              for( ; index < nnz->unrecorded_access(); index++ ) {
//                      while( index == row_i ) {
//                              ++row;
//                              row_i = row_start->const_access( row + 1 );
//                      }
// #ifdef _DEBUG
//                      std::cout << "index: " << index << " nzs[ index ]: " << nzs->unrecorded_access( index ) << " col_ind[ index ]: " << col_ind->unrecorded_access( index ) << " x[ col_ind[ index ] ]: " << x->unrecorded_access( col_ind->unrecorded_access( index ) ) << std::endl;
// #endif
// 
// #ifdef _DEBUG
//                      double toAdd = nzs->const_access( index );
//                      const unsigned long int temp = col_ind->const_access(index );
//                      toAdd *= x->const_access( temp );
// #else
//                      const double toAdd = nzs->const_access( index ) * x->const_access( col_ind->const_access( index ) );
// #endif
//                      ret->access( row ) += toAdd;
//              }
                for( ULI i=0; i < nor->unrecorded_access(); i++ ) {
                        for( ULI j=row_start->const_access( i ); j<row_start->const_access( i+1 ); j++ )
                                ret->access( i ) += nzs->const_access( j ) * x->const_access( col_ind->const_access( j ) );
                }
                T* real_ret = new T[ nor->unrecorded_access() ];
                ret->unrecordedCopyTo( real_ret );
                delete ret;
                delete x;
                return real_ret;
        }

        ~CS_CRS() {
                delete row_start;
                delete zero_element;
                delete nzs;
                delete col_ind;
                delete nnz;
                delete nor;
                delete noc;
        }

};

#endif

---