CS_ARRAY.hpp

#include<cstdlib>
#include<string>
#include<iostream>
#include<sstream>

#include "CACHE.hpp"
#include "CS_MATRIX.hpp"

#ifndef _H_CS_ARRAY
#define _H_CS_ARRAY

template < typename T > class CS_ARRAY {
   private:
        CS_ARRAY() {}

   protected:
        T* _array;

        unsigned long int _size;

        unsigned long int FAC;

   public:

        //static stuff first

        static std::string fDim( const CS_ARRAY& obj ) {
                std::ostringstream oss;
                oss << "( " << obj.getSize() << " x 1 )";
                return oss.str();
        }

        //end static stuff


        std::string _descr;

        CS_ARRAY( unsigned long int size, std::string descr ) {
                if ( descr == "" ) {
                        std::cerr << "Parameter description is mandatory during CS_ARRAY initialisation!" << std::endl;
                        exit( 1 );
                }

                if ( size == 0 ) {
                        std::cerr << "Size cannot be 0 during allocation of " << descr << std::endl;
                        exit( 1 );
                }

                _array = new T[size];
                _size = size;
                _descr = descr;
                FAC = 0;
        }

        ~CS_ARRAY() {
                delete [] _array;
        }

        unsigned long int size() { return _size; }

        inline
        T& operator[] ( unsigned long int i ) { return access( i ); }

        const T& const_access( unsigned long int i ) {
                if( i < _size ) {
                        CACHE::getInstance()->access( &( _array[ i ] ), sizeof( T ), FAC );
                        return _array[ i ];
                } else {
                        std::cerr << _descr << ": Index out of bounds! (" << i << "/" << this->_size << ")" << std::endl;
                        exit( !EXIT_SUCCESS );
                }
        }

        T& unrecorded_access( const unsigned long int i ) { return _array[ i ]; }

        const T& unrecorded_const_access( const unsigned long int i ) const { return _array[ i ]; }

        T& access( unsigned long int i ) {
                if( i < _size ) {
                        CACHE::getInstance()->access( &( _array[ i ] ), sizeof( T ), FAC );
                        return _array[ i ];
                } else {
                        std::cerr << _descr << ": Index out of bounds! (" << i << "/" << this->_size << ")" << std::endl;
                        exit( !EXIT_SUCCESS );
                }
        }

        unsigned long int getSize() const { return _size; }

        void unrecordedCopyTo( T* target ) {
                for( unsigned long int i=0; i<_size; i++ )
                        target[ i ] = _array[ i ];
        }

        void unrecordedCopyFrom( T* target ) {
                for( unsigned long int i=0; i<_size; i++ )
                        _array[ i ] = target[ i ];
        }

        void yaxpy( const CS_MATRIX< T >& A, const CS_ARRAY< T >& x ) {
                if( A.cols() != x.getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << x._descr << " (x) for y=Ax+y!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " times " << CS_ARRAY< T >::fDim( x ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                if( A.rows() != getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << _descr << " (y) for y=Ax+y!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " times " << CS_ARRAY< T >::fDim( *this ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                for( unsigned int i=0; i<A.rows(); i++ )
                        for( unsigned int j=0; j<A.cols(); j++ ) {
                                //preserve order of data access
                                const T _a = A( i, j );
                                const T _x = x[ j ];
                                const T _y = access( i );
                                access( i ) = _a * _x + _y;
                        }
        }

        void zaxpy( const CS_MATRIX< T >& A, const CS_ARRAY< T >& x, const CS_ARRAY< T >& y ) {
                if( A.cols() != x.getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << x._descr << " (x) for z=Ax+y!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " X " << CS_ARRAY< T >::fDim( x ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                if( A.rows() != y.getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << y._descr << " (y) for z=Ax+y!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " X " << CS_ARRAY< T >::fDim( y ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                if( A.rows() != getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << _descr << " (z) for z=Ax+y!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " X " << CS_ARRAY< T >::fDim( this ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                if( y.getSize() != getSize() ) {
                        std::cerr << y._descr << " (y) does not agree with " << _descr << " (z) for z=Ax+y!" << std::endl;
                        std::cerr << CS_ARRAY< T >::fDim( y ) << " X " << CS_ARRAY< T >::fDim( *this ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }
 
                for( unsigned int i=0; i<A.rows(); i++ ) {
                        //preserve order of access
                        const T ao = A( i, 0 );
                        const T xo = x[ 0 ];
                        const T yo = y [ i ];
                        access( i ) = ao*xo+yo; // z[i] = A[i][0]*x[0]+y[i]
                        for( unsigned int j=1; j<A.cols(); j++ ) {
                                const T _a = A( i, j );
                                const T _x = x[ j ];
                                const T _z = const_access( i );
                                access( i ) = _a*_x + _z; // z[i] += A[i][j]*x[j]
                        }
                }
        }

        void zax( const CS_MATRIX< T >& A, const CS_ARRAY< T >& x ) {
                if( A.cols() != x.getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << x._descr << " (x) for z=Ax!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " X " << CS_ARRAY< T >::fDim( x ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                if( A.rows() != getSize() ) {
                        std::cerr << A._descr << " (A) does not agree with " << _descr << " (z) for z=Ax!" << std::endl;
                        std::cerr << CS_MATRIX< T >::fDim( A ) << " X " << CS_ARRAY< T >::fDim( *this ) << std::endl;
                        exit( !EXIT_SUCCESS );
                }

                for( unsigned int i=0; i<A.rows(); i++ ) {
                        //preserve order of access
                        const T _a = A.const_access( i, 0 );
                        const T _x = x.const_access( 0 );
                        access( i ) = _a*_x;
                        for( unsigned int j=1; j<A.cols(); j++ ) {
                                const T _a = A.const_access( i, j );
                                const T _x = x.const_access( j );
                                const T _z = access( i );
                                access( i ) = _a * _x + _z;
                        }
                }
        }

};

#endif

---