/***************************************************************************
                          AdEllpack_impl.h  -  description
                             -------------------
    begin                : Aug 27, 2018
    copyright            : (C) 2018 by Tomas Oberhuber
    email                : tomas.oberhuber@fjfi.cvut.cz
 ***************************************************************************/

/* See Copyright Notice in tnl/Copyright */

#include <TNL/Matrices/AdEllpack.h>
#include <TNL/Containers/Vector.h>
#include <TNL/Math.h>
#include <TNL/param-types.h>

#pragma once

namespace TNL {
namespace Matrices {

/*
 * Auxiliary list implementation
 */
warpList::warpList()
{
    this->head = new warpInfo;
    this->tail = new warpInfo;
    this->head->previous = NULL;
    this->head->next = this->tail;
    this->tail->previous = this->head;
    this->tail->next = NULL;

    this->numberOfWarps = 0;
}

bool warpList::addWarp( const int offset,
                        const int rowOffset,
                        const int localLoad,
                        const int* reduceMap )
{
    warpInfo* temp = new warpInfo();
    if( !temp )
        return false;
    temp->offset = offset;
    temp->rowOffset = rowOffset;
    temp->localLoad = localLoad;
    for( int i = 0; i < 32; i++ )
        temp->reduceMap[ i ] = reduceMap[ i ];
    temp->next = this->tail;
    temp->previous = this->tail->previous;
    temp->previous->next = temp;
    this->tail->previous = temp;

    this->numberOfWarps++;

    return true;
}

warpInfo* warpList::splitInHalf( warpInfo* warp )
{
    warpInfo* firstHalf = new warpInfo();
    warpInfo* secondHalf = new warpInfo();
    int localLoad = ( warp->localLoad / 2 ) + ( warp->localLoad % 2 == 0 ? 0 : 1 );

    int rowOffset = warp->rowOffset;

    // first half split
    firstHalf->localLoad = localLoad;
    firstHalf->rowOffset = warp->rowOffset;
    firstHalf->offset = warp->offset;

    firstHalf->reduceMap[ 0 ] = 1;
    firstHalf->reduceMap[ 1 ] = 0;
    for( int i = 1; i < 16; i++ )
    {
        if( warp->reduceMap[ i ] == 1 )
        {
            rowOffset++;
            firstHalf->reduceMap[ 2 * i ] = 1;
            firstHalf->reduceMap[ 2 * i + 1 ] = 0;
        }
        else
        {
            firstHalf->reduceMap[ 2 * i ] = 0;
            firstHalf->reduceMap[ 2 * i + 1 ] = 0;
        }
    }

    // second half split
    secondHalf->rowOffset = rowOffset;
    if( warp->reduceMap[ 16 ] == 1 )
        secondHalf->rowOffset = rowOffset + 1;
    secondHalf->offset = 32 * firstHalf->localLoad + firstHalf->offset;
    secondHalf->reduceMap[ 0 ] = 1;
    secondHalf->reduceMap[ 1 ] = 0;
    for( int i = 1; i < 16; i++ )
    {
        if( warp->reduceMap[ i + 16 ] == 1 )
        {
            secondHalf->reduceMap[ 2 * i ] = 1;
            secondHalf->reduceMap[ 2 * i + 1 ] = 0;
        }
        else
        {
            secondHalf->reduceMap[ 2 * i ] = 0;
            secondHalf->reduceMap[ 2 * i + 1 ] = 0;
        }
    }
    secondHalf->localLoad = localLoad;

    // and warps must be connected to the list
    firstHalf->next = secondHalf;
    secondHalf->previous = firstHalf;
    firstHalf->previous = warp->previous;
    warp->previous->next = firstHalf;
    secondHalf->next = warp->next;
    warp->next->previous = secondHalf;

    // if original load was odd, all next warp offsets must be changed!
    if( ( warp->localLoad % 2 ) != 0 )
    {
        warp = secondHalf->next;
        while( warp != this->tail )
        {
            warp->offset = warp->previous->offset + warp->previous->localLoad * 32;
            warp = warp->next;
        }
    }
    this->numberOfWarps++;

    // method returns the first of the two splitted warps
    return firstHalf;
}

warpList::~warpList()
{
    while( this->head->next != NULL )
    {
        warpInfo* temp = new warpInfo;
        temp = this->head->next;
        this->head->next = temp->next;
        delete temp;
    }
    delete this->head;
}


/*
 * Adaptive Ellpack implementation
 */

template< typename Real,
          typename Device,
          typename Index >
AdEllpack< Real, Device, Index >::AdEllpack()
:
warpSize( 32 )
{}

template< typename Real,
          typename Device,
          typename Index >
String AdEllpack< Real, Device, Index >::getTypeVirtual() const
{
    return this->getType();
}

template< typename Real,
          typename Device,
          typename Index >
String AdEllpack< Real, Device, Index >::getType()
{
    return String( "AdEllpack< ") +
           String( TNL::getType< Real >() ) +
           String( ", " ) +
           Device::getDeviceType() +
           String( " >" );
}

template< typename Real,
          typename Device,
          typename Index >
void
AdEllpack< Real, Device, Index >::
setCompressedRowLengths( ConstCompressedRowLengthsVectorView rowLengths )
{
    TNL_ASSERT( this->getRows() > 0, );
    TNL_ASSERT( this->getColumns() > 0, );
    if( std::is_same< DeviceType, Devices::Host >::value )
    {
         RealType average = 0.0;
         for( IndexType row = 0; row < this->getRows(); row++ )
            average += rowLengths.getElement( row );
         average /= ( RealType ) this->getRows();
         this->totalLoad = average;

        warpList* list = new warpList();

        if( !this->balanceLoad( average, rowLengths, list ) )
            throw 0; // TODO: Make better exception

        IndexType SMs = 15;
        IndexType threadsPerSM = 2048;

        this->computeWarps( SMs, threadsPerSM, list );

        if( !this->createArrays( list ) )
            throw 0; // TODO: Make better excpetion



        //this->performRowTest();
        //cout << "========================" << std::endl;
        //cout << "Testing row lengths" << std::endl;
        //cout << "========================" << std::endl;
        //this->performRowLengthsTest( rowLengths );
    }

    if( std::is_same< DeviceType, Devices::Cuda >::value )
    {
        AdEllpack< RealType, Devices::Host, IndexType > hostMatrix;
        hostMatrix.setDimensions( this->getRows(), this->getColumns() );
        Containers::Vector< IndexType, Devices::Host, IndexType > hostRowLengths;
        hostRowLengths.setLike( rowLengths );
        hostRowLengths = rowLengths;
        hostMatrix.setCompressedRowLengths( hostRowLengths );

        this->offset.setLike( hostMatrix.offset );
        this->offset = hostMatrix.offset;
        this->rowOffset.setLike( hostMatrix.rowOffset );
        this->rowOffset = hostMatrix.rowOffset;
        this->localLoad.setLike( hostMatrix.localLoad );
        this->localLoad = hostMatrix.localLoad;
        this->reduceMap.setLike( hostMatrix.reduceMap );
        this->reduceMap = hostMatrix.reduceMap;
         this->totalLoad = hostMatrix.getTotalLoad();

        this->allocateMatrixElements( this->offset.getElement( this->offset.getSize() - 1 ) );
    }
}

template< typename Real,
          typename Device,
          typename Index >
Index AdEllpack< Real, Device, Index >::getTotalLoad() const
{
    return this->totalLoad;
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::performRowLengthsTest( ConstCompressedRowLengthsVectorView rowLengths )
{
    bool found = false;
    for( IndexType row = 0; row < this->getRows(); row++ )
    {
        found = false;
        IndexType warp = this->getWarp( row );
        IndexType inWarpOffset = this->getInWarpOffset( row, warp );
        if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
            warp++;
        IndexType rowLength = this->localLoad.getElement( warp );
        while( !found )
        {
            if( ( inWarpOffset < this->warpSize - 1 ) &&
                ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
		{
                    inWarpOffset++;
		    rowLength += this->localLoad.getElement( warp );
		}
            else if( ( inWarpOffset == this->warpSize - 1 ) &&
                     ( this->rowOffset.getElement( warp + 1 ) == row ) )
            {
                warp++;
                inWarpOffset = 0;
		rowLength += localLoad.getElement( warp );
            }
            else
                found = true;
        }
	if( rowLength < rowLengths.getElement( row ) )
	    std::cout << "Row: " << row << " is short!";
    }
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::performRowTest()
{
    for( IndexType warp = 0; warp < this->localLoad.getSize(); warp++ )
    {
	IndexType row = this->rowOffset.getElement( warp );
        for( IndexType i = warp * this->warpSize + 1; i < ( warp + 1 ) * this->warpSize; i++ )
	{
	    if( this->reduceMap.getElement( i ) == 1 )
		row++;
	}
	if( row == this->rowOffset.getElement( warp + 1 ) || row + 1 == this->rowOffset.getElement( warp + 1 ) )
	    ;
	else 
        {
	    std::cout << "Error warp = " << warp << std::endl;
	    std::cout << "Row: " << row << ", Row offset: " << this->rowOffset.getElement( warp + 1 ) << std::endl;
        }
    }
}

template< typename Real,
          typename Device,
          typename Index >
Index AdEllpack< Real, Device, Index >::getWarp( const IndexType row ) const
{
    for( IndexType searchedWarp = 0; searchedWarp < this->getRows(); searchedWarp++ )
    {
        if( ( this->rowOffset.getElement( searchedWarp ) == row ) ||
            ( ( this->rowOffset.getElement( searchedWarp ) < row ) && ( this->rowOffset.getElement( searchedWarp + 1 ) >= row ) ) )
            return searchedWarp;
    }
    // FIXME: non-void function always has to return something sensible
    throw "bug - row was not found";
}

template< typename Real,
          typename Device,
          typename Index >
Index AdEllpack< Real, Device, Index >::getInWarpOffset( const IndexType row,
                                                                  const IndexType warp ) const
{
    IndexType inWarpOffset = warp * this->warpSize;
    while( ( inWarpOffset < ( warp + 1 ) * this->warpSize ) && ( this->reduceMap.getElement( inWarpOffset ) < row ) )
        inWarpOffset++;
    return ( inWarpOffset & ( this->warpSize - 1 ) );
}

template< typename Real,
          typename Device,
          typename Index >
Index AdEllpack< Real, Device, Index >::getRowLength( const IndexType row ) const
{
    IndexType warp = this->getWarp( row );
    IndexType inWarpOffset = this->getInWarpOffset( row, warp );
    if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
        warp++;

    bool found = false;
    IndexType length = 0;
    IndexType elementPtr;
    while( !found )
    {
        elementPtr = this->offset.getElement( warp ) + inWarpOffset;
        for( IndexType i = 0; i < this->localLoad.getElement( warp ); i++ )
        {
            if( this->columnIndexes.getElement( elementPtr ) != this->getPaddingIndex() )
                length++;
            elementPtr += this->warpSize;
        }
        if( ( inWarpOffset < this->warpSize - 1 ) &&
            ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
            inWarpOffset++;
        else if( ( inWarpOffset == this->warpSize - 1 ) &&
                 ( this->rowOffset.getElement( warp + 1 ) == row ) )
        {
            warp++;
            inWarpOffset = 0;
        }
        else
            found = true;
    }
    return length;
}

template< typename Real,
          typename Device,
          typename Index >
template< typename Real2,
          typename Device2,
          typename Index2 >
bool AdEllpack< Real, Device, Index >::setLike( const AdEllpack< Real2, Device2, Index2 >& matrix )
{
    if( !Sparse< Real, Device, Index >::setLike( matrix ) ||
        !this->offset.setLike( matrix.offset ) ||
        !this->rowOffset.setLike( matrix.rowOffset ) ||
        !this->localLoad.setLike( matrix.localLoad ) ||
        !this->reduceMap.setLike( matrix.reduceMap ) )
        return false;
    return true;
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::reset()
{
    Sparse< Real, Device, Index >::reset();
    this->offset.reset();
    this->rowOffset.reset();
    this->localLoad.reset();
    this->reduceMap.reset();
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::setElement( const IndexType row,
                                                            const IndexType column,
                                                            const RealType& value )
{
    return this->addElement( row, column, value, 0.0 );
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::addElement( const IndexType row,
                                                            const IndexType column,
                                                            const RealType& value,
                                                            const RealType& thisElementMultiplicator )
{
    IndexType warp = this->getWarp( row );
    IndexType inWarpOffset = this->getInWarpOffset( row, warp );
    if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
        warp++;

    bool found = false;
    IndexType elementPtr;
    IndexType iterator = 0;
    while( !found )
    {
        elementPtr = this->offset.getElement( warp ) + inWarpOffset;
        for( IndexType i = 0; i < this->localLoad.getElement( warp ); i++ )
        {
            iterator++;
            if( this->columnIndexes.getElement( elementPtr ) == column )
            {
                this->values.setElement( elementPtr, thisElementMultiplicator * this->values.getElement( elementPtr ) + value );
                return true;
            }
            if( this->columnIndexes.getElement( elementPtr ) == this->getPaddingIndex() )
            {
                this->columnIndexes.setElement( elementPtr, column );
                this->values.setElement( elementPtr, value );
                return true;
            }
            elementPtr += this->warpSize;
        }
        if( ( inWarpOffset < this->warpSize - 1 ) &&
            ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
	{
            inWarpOffset++;
	}
        else if( ( inWarpOffset == this->warpSize - 1 ) &&
                 ( this->rowOffset.getElement( warp + 1 ) == row ) )
        {
            warp++;
            inWarpOffset = 0;
        }
        else
            found = true;
    }
    return false;
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::setRow( const IndexType row,
                                                        const IndexType* columnIndexes,
                                                        const RealType* values,
                                                        const IndexType elements )
{
    IndexType warp = this->getWarp( row );
    IndexType inWarpOffset = this->getInWarpOffset( row, warp );
    if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
        warp++;

    bool found = false;
    IndexType elementPtr;
    IndexType elPtr = 0;
    while( ( !found ) && ( elPtr < elements ) )
    {
        elementPtr = this->offset.getElement( warp ) + inWarpOffset;
        for( IndexType i = 0; ( i < this->localLoad.getElement( warp ) ) && ( elPtr < elements ); i++ )
        {
            this->values.setElement( elementPtr, values[ elPtr ] );
            this->columnIndexes.setElement( elementPtr, columnIndexes[ elPtr ] );
            elPtr++;
            elementPtr += this->warpSize;
        }
        if( ( inWarpOffset < this->warpSize - 1 ) &&
            ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
            inWarpOffset++;
        else if( ( inWarpOffset == this->warpSize - 1 ) &&
                 ( this->rowOffset.getElement( warp + 1 ) == row ) )
        {
            warp++;
            inWarpOffset = 0;
        }
        else
            found = true;
    }
    return true;
}

template< typename Real,
          typename Device,

          typename Index >
bool AdEllpack< Real, Device, Index >::addRow( const IndexType row,
                                                        const IndexType* columnIndexes,
                                                        const RealType* values,
                                                        const IndexType elements,
                                                        const RealType& thisElementMultiplicator )
{
    return false;
}

template< typename Real,
          typename Device,
          typename Index >
Real AdEllpack< Real, Device, Index >::getElement( const IndexType row,
                                                            const IndexType column ) const
{
    IndexType warp = this->getWarp( row );
    IndexType inWarpOffset = this->getInWarpOffset( row, warp );
    if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
        warp++;

    bool found = false;
    IndexType elementPtr;
    while( !found )
    {
        elementPtr = this->offset.getElement( warp ) + inWarpOffset;
        for( IndexType i = 0; i < this->localLoad.getElement( warp ); i++ )
        {
            if( this->columnIndexes.getElement( elementPtr ) == column )
                return this->values.getElement( elementPtr );
            elementPtr += this->warpSize;
        }
        if( ( inWarpOffset < this->warpSize - 1 ) &&
            ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
            inWarpOffset++;
        else if( ( inWarpOffset == this->warpSize - 1 ) &&
                 ( this->rowOffset.getElement( warp + 1 ) == row ) )
        {
            warp++;
            inWarpOffset = 0;
        }
        else
            found = true;
    }
    return 0.0;
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::getRow( const IndexType row,
                                                        IndexType* columns,
                                                        RealType* values ) const
{
    IndexType warp = this->getWarp( row );
    IndexType inWarpOffset = this->getInWarpOffset( row, warp );
    if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
        warp++;

    bool found = false;
    IndexType arrayElementPtr = 0;
    IndexType elementPtr;
    while( !found )
    {
        elementPtr = this->offset.getElement( warp ) + inWarpOffset;
        for( IndexType i = 0; i < this->localLoad.getElement( warp ); i++ )
        {
            if( this->columnIndexes.getElement( elementPtr ) != this->getPaddingIndex() )
            {
                columns[ arrayElementPtr ] = this->columnIndexes.getElement( elementPtr );
                values[ arrayElementPtr ] = this->values.getElement( elementPtr );
                arrayElementPtr++;
            }
            elementPtr += this->warpSize;
        }
        if( ( inWarpOffset < this->warpSize - 1 ) &&
            ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
            inWarpOffset++;
        else if( ( inWarpOffset == this->warpSize - 1 ) &&
                 ( this->rowOffset.getElement( warp + 1 ) == row ) )
        {
            warp++;
            inWarpOffset = 0;
        }
        else
            found = true;
    }
}

template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
void AdEllpack< Real, Device, Index >::vectorProduct( const InVector& inVector,
                                                               OutVector& outVector ) const
{
    DeviceDependentCode::vectorProduct( *this, inVector, outVector );
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::save( File& file ) const
{
    if( Sparse< Real, Device, Index >::save( file ) ||
        this->offset.save( file ) ||
        this->rowOffset.save( file ) ||
        this->localLoad.save( file ) ||
        this->reduceMap.save( file ) )
        return false;
    return true;
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::load( File& file )
{
    if( Sparse< Real, Device, Index >::load( file ) ||
        this->offset.load( file ) ||
        this->rowOffset.load( file ) ||
        this->localLoad.load( file ) ||
        this->reduceMap.load( file ) )
        return false;
    return true;
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::save( const String& fileName ) const
{
    return Object::save( fileName );
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::load( const String& fileName )
{
    return Object::load( fileName );
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::print( std::ostream& str ) const
{
    for( IndexType row = 0; row < this->getRows(); row++ )
    {
        str  << "Row: " << row << " -> \t";

        IndexType warp = this->getWarp( row );
        IndexType inWarpOffset = this->getInWarpOffset( row, warp );
        if( inWarpOffset == 0 && rowOffset.getElement( warp ) != row )
            warp++;

        bool found = false;
        IndexType elementPtr;
        while( !found )
        {
            elementPtr = this->offset.getElement( warp ) + inWarpOffset;
            for( IndexType i = 0; i < this->localLoad.getElement( warp ); i++ )
            {
                if( this->columnIndexes.getElement( elementPtr ) != this->getPaddingIndex() )
                    str << " column: " << this->columnIndexes.getElement( elementPtr ) << " -> "
                        << " value: " << this->values.getElement( elementPtr ) << std::endl;
                elementPtr += this->warpSize;
            }
            if( ( inWarpOffset < this->warpSize - 1 ) &&
                ( this->reduceMap.getElement( this->warpSize * warp + inWarpOffset + 1 ) == row ) )
                inWarpOffset++;
            else if( ( inWarpOffset == this->warpSize - 1 ) &&
                     ( this->rowOffset.getElement( warp + 1 ) == row ) )
            {
                warp++;
                inWarpOffset = 0;
            }
            else
                found = true;
        }
    }
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::balanceLoad( const RealType average,
                                                    ConstCompressedRowLengthsVectorView rowLengths,
                                                    warpList* list )
{
    IndexType offset, rowOffset, localLoad, reduceMap[ 32 ];

    IndexType numberOfThreads = 0;
    IndexType ave = average * 1.1;
    offset = 0;
    rowOffset = 0;
    localLoad = 0;
    bool addedWarp = false;
    IndexType warpId = 0;

    for( IndexType row = 0; row < this->getRows(); row++ )
    {
        addedWarp = false;
        // takes every row
        IndexType rowLength = rowLengths.getElement( row );
        // counts necessary number of threads
        IndexType threadsPerRow = ( IndexType ) ( rowLength / ave ) + ( rowLength % ave == 0 ? 0 : 1 );
        // if this number added to number of total threads in this warp is higher than 32
        // then spreads them in as many as necessary

        while( numberOfThreads + threadsPerRow >= this->warpSize )
        {
            // counts number of usable threads
            IndexType usedThreads = this->warpSize - numberOfThreads;

            // if this condition applies we can finish the warp
            if( usedThreads == threadsPerRow )
            {
                localLoad = max( localLoad, roundUpDivision( rowLength, threadsPerRow ) );

                reduceMap[ numberOfThreads ] = 1;
                for( IndexType i = numberOfThreads + 1; i < this->warpSize; i++ )
                    reduceMap[ i ] = 0;

                if( !list->addWarp( offset, rowOffset, localLoad, reduceMap ) )
                    return false;

                offset += this->warpSize * localLoad;
                rowOffset = row + 1;

                // IMPORTANT TO RESET VARIABLES
                localLoad = 0;
                for( IndexType i = 0; i < this->warpSize; i++ )
                    reduceMap[ i ] = 0;
                numberOfThreads = 0;
                threadsPerRow = 0;
                addedWarp = true;
		warpId++;
            }

            // if it doesnt apply and if local load isn't equal 0 it will use old local load ( for better balance )
            else if( localLoad != 0 )
            {
                // subtract unmber of used elements and number of used threads
                rowLength -= localLoad * usedThreads;

                threadsPerRow = ( IndexType ) ( rowLength / ave ) + ( rowLength % ave == 0 ? 0 : 1 );

                // fill the reduction map
                reduceMap[ numberOfThreads ] = 1;
                for( IndexType i = numberOfThreads + 1; i < this->warpSize; i++ )
                    reduceMap[ i ] = 0;

                // count new offsets, add new warp and reset variables
                if( !list->addWarp( offset, rowOffset, localLoad, reduceMap ) )
                    return false;
                offset += this->warpSize * localLoad;
                rowOffset = row;

                // RESET VARIABLES
                localLoad = 0;
                for( IndexType i = 0; i < this->warpSize; i++ )
                    reduceMap[ i ] = 0;
                numberOfThreads = 0;
                addedWarp = true;
		warpId++;
            }
            // otherwise we are starting a new warp
            else
            {
                threadsPerRow = ( IndexType ) ( rowLength / ave ) + ( rowLength % ave == 0 ? 0 : 1 );
                if( threadsPerRow < this->warpSize )                
                    break;

                localLoad = ave;
            }
        }
	if( threadsPerRow <= 0 )
	{
	    threadsPerRow = 1;
	    continue;
	}
        localLoad = max( localLoad, roundUpDivision( rowLength, threadsPerRow ) );
        reduceMap[ numberOfThreads ] = 1;
        for( IndexType i = numberOfThreads + 1; i < numberOfThreads + threadsPerRow; i++ )
            reduceMap[ i ] = 0;

        numberOfThreads += threadsPerRow;

        // if last row doesnt fill the whole warp or it fills more warps and threads still remain
        if( ( ( row == this->getRows() - 1 ) && !addedWarp ) ||
            ( ( row == this->getRows() - 1 ) && ( threadsPerRow == numberOfThreads ) && ( numberOfThreads > 0 ) ) )
        {
            list->addWarp( offset, rowOffset, localLoad, reduceMap );
        }
    }
    return true;
}

template< typename Real,
          typename Device,
          typename Index >
void AdEllpack< Real, Device, Index >::computeWarps( const IndexType SMs,
                                                              const IndexType threadsPerSM,
                                                              warpList* list )
{
    IndexType averageLoad = 0;
    warpInfo* temp = list->getHead()->next;
    while( temp->next != list->getTail() )
    {
        averageLoad += temp->localLoad;
        temp = temp->next;
    }
    averageLoad /= list->getNumberOfWarps();

    IndexType totalWarps = SMs * ( threadsPerSM / this->warpSize );
    IndexType remainingThreads = list->getNumberOfWarps();
    bool warpsToSplit = true;

    while( remainingThreads < ( totalWarps / 2 ) && warpsToSplit )
    {
        warpsToSplit = false;
        temp = list->getHead()->next;
        while( temp != list->getTail() )
        {
            if( temp->localLoad > averageLoad )
            {
                temp = list->splitInHalf( temp );
                warpsToSplit = true;
		
            }
            temp = temp->next;
        }
	remainingThreads = list->getNumberOfWarps();
    }
}

template< typename Real,
          typename Device,
          typename Index >
bool AdEllpack< Real, Device, Index >::createArrays( warpList* list )
{
    IndexType length = list->getNumberOfWarps();

    this->offset.setSize( length + 1 );
    this->rowOffset.setSize( length + 1 );
    this->localLoad.setSize( length );
    this->reduceMap.setSize( length * this->warpSize );

    IndexType iteration = 0;
    warpInfo* warp = list->getHead()->next;
    while( warp != list->getTail() )
    {
        this->offset.setElement( iteration, warp->offset );
        this->rowOffset.setElement( iteration, warp->rowOffset );
        this->localLoad.setElement( iteration, warp->localLoad );
	IndexType row = this->rowOffset.getElement( iteration );
	this->reduceMap.setElement( iteration * this->warpSize, row );
        for( int i = iteration * this->warpSize + 1; i < ( iteration + 1 ) * this->warpSize; i++ )
	{
	    if( warp->reduceMap[ i & ( this->warpSize - 1 ) ] == 1 )
		row++;
            this->reduceMap.setElement( i, row );
	}
	if( localLoad > this->totalLoad )
	    std::cout << "Error localLoad!!" << std::endl;
        iteration++;
        warp = warp->next;
    }
    this->rowOffset.setElement( length, this->getRows() );
    this->offset.setElement( length, this->offset.getElement( length - 1 ) + this->warpSize * this->localLoad.getElement( length - 1 ) );
    this->allocateMatrixElements( this->offset.getElement( length ) );

    return true;
}

template<>
class AdEllpackDeviceDependentCode< Devices::Host >
{
public:

    typedef Devices::Host Device;

    template< typename Real,
              typename Index,
              typename InVector,
              typename OutVector >
    static void vectorProduct( const AdEllpack< Real, Device, Index >& matrix,
                               const InVector& inVector,
                               OutVector& outVector )
    {
	// parallel vector product simulation
	const Index blockSize = 256; 
	const Index blocks = ( Index ) ( matrix.reduceMap.getSize() / blockSize ) + ( matrix.reduceMap.getSize() % blockSize != 0 );
	for( Index block = 0; block < blocks; block++ )
	{
	    Containers::Vector< Real, Device, Index > temp, reduceMap;
	    temp.setSize( blockSize );
	    reduceMap.setSize( blockSize );
	    for( Index threadIdx = 0; threadIdx < blockSize; threadIdx++ )
	    {
		Index globalIdx = block * blockSize + threadIdx;
		Index warpIdx = globalIdx >> 5;
		Index inWarpIdx = globalIdx & ( matrix.warpSize - 1 );
		if( globalIdx >= matrix.reduceMap.getSize() )
		{
		    for( Index i = globalIdx % blockSize; i < blockSize; i++ )
			reduceMap.setElement( i, -1 );
		    break;
		}
		temp.setElement( threadIdx, 0.0 );
		reduceMap.setElement( threadIdx, matrix.reduceMap.getElement( globalIdx ) );
		Index elementPtr = matrix.offset.getElement( warpIdx ) + inWarpIdx;
		for( Index i = 0; i < matrix.localLoad.getElement( warpIdx ); i++ )
		{
		    if( matrix.columnIndexes.getElement( elementPtr ) != matrix.getPaddingIndex() )
		        temp.setElement( threadIdx, temp.getElement( threadIdx ) + matrix.values.getElement( elementPtr ) * inVector[ matrix.columnIndexes.getElement( elementPtr ) ] );
		    elementPtr += matrix.warpSize;
		}
	    }
	    for( Index threadIdx = 0; threadIdx < blockSize; threadIdx++ )
	    {
		Index globalIdx = block * blockSize + threadIdx;
		Index warpIdx = globalIdx >> 5;
		Index inWarpIdx = globalIdx & ( matrix.warpSize - 1 );
		if( globalIdx >= matrix.reduceMap.getSize() )
		    break;
		if( ( inWarpIdx == 0 ) || ( reduceMap.getElement( threadIdx - 1 ) < reduceMap.getElement( threadIdx ) ) )
		{
		    Index end = ( warpIdx + 1 ) << 5;
		    Index elementPtr = threadIdx + 1;
		    globalIdx++;
		    while( globalIdx < end && reduceMap.getElement( elementPtr ) == reduceMap.getElement( threadIdx ) )
		    {
			temp.setElement( threadIdx, temp.getElement( elementPtr ) + temp.getElement( threadIdx ) );
			elementPtr++;
			globalIdx++;
		    }
		    outVector[ reduceMap.getElement( threadIdx ) ] += temp.getElement( threadIdx );
		}
	    }
	}
    }

};

#ifdef HAVE_CUDA
template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
__device__
void AdEllpack< Real, Device, Index >::spmvCuda2( const InVector& inVector,
                                                  OutVector& outVector,
                                                  const int gridIdx ) const
{
    IndexType globalIdx = ( gridIdx * Devices::Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
    IndexType warpIdx = globalIdx >> 5;
    IndexType inWarpIdx = globalIdx & ( this->warpSize - 1 );
    if( globalIdx >= this->reduceMap.getSize() )
	return;

    const int blockSize = 256;
    Real* temp = Devices::Cuda::getSharedMemory< Real >();
    __shared__ IndexType reduceMap[ blockSize ];
    reduceMap[ threadIdx.x ] = this->reduceMap[ globalIdx ];
    temp[ threadIdx.x ] = 0.0;

        IndexType i = 0;
        IndexType elementPtr = this->offset[ warpIdx ] + inWarpIdx;
        for( ; i < this->localLoad[ warpIdx ]; i++ )
        {
	    if( this->columnIndexes[ elementPtr ] < this->getColumns() )
            {
	        temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                elementPtr += this->warpSize;
	    }
        }
    if( ( inWarpIdx == 0 ) || ( reduceMap[ threadIdx.x ] > reduceMap[ threadIdx.x - 1 ] ) )
    {
	IndexType elementPtr = threadIdx.x + 1;
        globalIdx++;
        while( //elementPtr < this->reduceMap.getSize() && 
	       globalIdx < ( ( warpIdx + 1 ) << 5 ) && 
               reduceMap[ elementPtr ] == reduceMap[ threadIdx.x ] )
        {
		temp[ threadIdx.x ] += temp[ elementPtr ];
                elementPtr++;
                globalIdx++;
        }
        outVector[ reduceMap[ threadIdx.x] ] += temp[ threadIdx.x ];
    }
} 

template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
__device__
void AdEllpack< Real, Device, Index >::spmvCuda4( const InVector& inVector,
                                                           OutVector& outVector,
                                                           const int gridIdx ) const
{
    IndexType globalIdx = ( gridIdx * Devices::Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
    IndexType warpIdx = globalIdx >> 5;
    IndexType inWarpIdx = globalIdx & ( this->warpSize - 1 );
    if( globalIdx >= this->reduceMap.getSize() )
	return;

    const int blockSize = 192;
    Real* temp = Devices::Cuda::getSharedMemory< Real >();
    __shared__ IndexType reduceMap[ blockSize ];
    reduceMap[ threadIdx.x ] = this->reduceMap[ globalIdx ];
    temp[ threadIdx.x ] = 0.0;

        IndexType i = 0;
        IndexType elementPtr = this->offset[ warpIdx ] + inWarpIdx;

        if( ( this->localLoad[ warpIdx ] & 1 ) == 1 )	
	    if( this->columnIndexes[ elementPtr ] < this->getColumns() )
	    {
	        temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                elementPtr += this->warpSize;
                i++;
	    }
        for( ; i < this->localLoad[ warpIdx ]; i += 2 )
        {
	    #pragma unroll
            for( IndexType j = 0; j < 2; j++ )
	        if( this->columnIndexes[ elementPtr ] < this->getColumns() )
                {
	            temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                    elementPtr += this->warpSize;
	        }
        }

    if( ( inWarpIdx == 0 ) || ( reduceMap[ threadIdx.x ] > reduceMap[ threadIdx.x - 1 ] ) )
    {
	IndexType elementPtr = threadIdx.x + 1;
        globalIdx++;
        while( //elementPtr < this->reduceMap.getSize() && 
	       globalIdx < ( ( warpIdx + 1 ) << 5 ) && 
               reduceMap[ elementPtr ] == reduceMap[ threadIdx.x ] )
        {
		temp[ threadIdx.x ] += temp[ elementPtr ];
                elementPtr++;
                globalIdx++;
        }
        outVector[ reduceMap[ threadIdx.x] ] += temp[ threadIdx.x ];
    }
}

template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
__device__
void AdEllpack< Real, Device, Index >::spmvCuda8( const InVector& inVector,
                                                           OutVector& outVector,
                                                           const int gridIdx ) const
{
    IndexType globalIdx = ( gridIdx * Devices::Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
    IndexType warpIdx = globalIdx >> 5;
    IndexType inWarpIdx = globalIdx & ( this->warpSize - 1 );
    if( globalIdx >= this->reduceMap.getSize() )
	return;

    const int blockSize = 128;
    Real* temp = Devices::Cuda::getSharedMemory< Real >();
    __shared__ IndexType reduceMap[ blockSize ];
    reduceMap[ threadIdx.x ] = this->reduceMap[ globalIdx ];
    temp[ threadIdx.x ] = 0.0;

        IndexType i = 0;
        IndexType elementPtr = this->offset[ warpIdx ] + inWarpIdx;

        while( ( this->localLoad[ warpIdx ] & 7 ) != 0 )	
	    if( this->columnIndexes[ elementPtr ] < this->getColumns() )
	    {
	        temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                elementPtr += this->warpSize;
                i++;
	    }
        for( ; i < this->localLoad[ warpIdx ]; i += 4 )
        {
	    #pragma unroll
            for( IndexType j = 0; j < 4; j++ )
	        if( this->columnIndexes[ elementPtr ] < this->getColumns() )
                {
	            temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                    elementPtr += this->warpSize;
	        }
        }
    if( ( inWarpIdx == 0 ) || ( reduceMap[ threadIdx.x ] > reduceMap[ threadIdx.x - 1 ] ) )
    {
	IndexType elementPtr = threadIdx.x + 1;
        globalIdx++;
        while( //elementPtr < this->reduceMap.getSize() && 
	       globalIdx < ( ( warpIdx + 1 ) << 5 ) && 
               reduceMap[ elementPtr ] == reduceMap[ threadIdx.x ] )
        {
		temp[ threadIdx.x ] += temp[ elementPtr ];
                elementPtr++;
                globalIdx++;
        }
        outVector[ reduceMap[ threadIdx.x] ] += temp[ threadIdx.x ];
    }
}

template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
__device__
void AdEllpack< Real, Device, Index >::spmvCuda16( const InVector& inVector,
                                                            OutVector& outVector,
                                                            const int gridIdx ) const
{
    IndexType globalIdx = ( gridIdx * Devices::Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
    IndexType warpIdx = globalIdx >> 5;
    IndexType inWarpIdx = globalIdx & ( this->warpSize - 1 );
    if( globalIdx >= this->reduceMap.getSize() )
	return;

    const int blockSize = 128;
    Real* temp = Devices::Cuda::getSharedMemory< Real >();
    __shared__ IndexType reduceMap[ blockSize ];
    reduceMap[ threadIdx.x ] = this->reduceMap[ globalIdx ];
    temp[ threadIdx.x ] = 0.0;

        IndexType i = 0;
        IndexType elementPtr = this->offset[ warpIdx ] + inWarpIdx;

        while( ( this->localLoad[ warpIdx ] & 15 ) != 0 )	
	    if( this->columnIndexes[ elementPtr ] < this->getColumns() )
	    {
	        temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                elementPtr += this->warpSize;
                i++;
	    }
        for( ; i < this->localLoad[ warpIdx ]; i += 8 )
        {
	    #pragma unroll
            for( IndexType j = 0; j < 8; j++ )
	        if( this->columnIndexes[ elementPtr ] < this->getColumns() )
                {
	            temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                    elementPtr += this->warpSize;
	        }
        }
    if( ( inWarpIdx == 0 ) || ( reduceMap[ threadIdx.x ] > reduceMap[ threadIdx.x - 1 ] ) )
    {
	IndexType elementPtr = threadIdx.x + 1;
        globalIdx++;
        while( //elementPtr < this->reduceMap.getSize() && 
	       globalIdx < ( ( warpIdx + 1 ) << 5 ) && 
               reduceMap[ elementPtr ] == reduceMap[ threadIdx.x ] )
        {
		temp[ threadIdx.x ] += temp[ elementPtr ];
                elementPtr++;
                globalIdx++;
        }
        outVector[ reduceMap[ threadIdx.x] ] += temp[ threadIdx.x ];
    }
}

template< typename Real,
          typename Device,
          typename Index >
template< typename InVector,
          typename OutVector >
__device__
void AdEllpack< Real, Device, Index >::spmvCuda32( const InVector& inVector,
                                                            OutVector& outVector,
                                                            const int gridIdx ) const
{
    IndexType globalIdx = ( gridIdx * Devices::Cuda::getMaxGridSize() + blockIdx.x ) * blockDim.x + threadIdx.x;
    IndexType warpIdx = globalIdx >> 5;
    IndexType inWarpIdx = globalIdx & ( this->warpSize - 1 );
    if( globalIdx >= this->reduceMap.getSize() )
	return;

    const int blockSize = 96;
    Real* temp = Devices::Cuda::getSharedMemory< Real >();
    __shared__ IndexType reduceMap[ blockSize ];
    reduceMap[ threadIdx.x ] = this->reduceMap[ globalIdx ];
    temp[ threadIdx.x ] = 0.0;
        IndexType i = 0;
        IndexType elementPtr = this->offset[ warpIdx ] + inWarpIdx;

        while( ( this->localLoad[ warpIdx ] & 31 ) != 0 )	
	    if( this->columnIndexes[ elementPtr ] < this->getColumns() )
	    {
	        temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                elementPtr += this->warpSize;
                i++;
	    }
        for( ; i < this->localLoad[ warpIdx ]; i += 16 )
        {
	    #pragma unroll
            for( IndexType j = 0; j < 16; j++ )
	        if( this->columnIndexes[ elementPtr ] < this->getColumns() )
                {
	            temp[ threadIdx.x ] += inVector[ this->columnIndexes[ elementPtr ] ] * this->values[ elementPtr ];
                    elementPtr += this->warpSize;
	        }
        }
    if( ( inWarpIdx == 0 ) || ( reduceMap[ threadIdx.x ] > reduceMap[ threadIdx.x - 1 ] ) )
    {
	IndexType elementPtr = threadIdx.x + 1;
        globalIdx++;
        while( //elementPtr < this->reduceMap.getSize() && 
	       globalIdx < ( ( warpIdx + 1 ) << 5 ) && 
               reduceMap[ elementPtr ] == reduceMap[ threadIdx.x ] )
        {
		temp[ threadIdx.x ] += temp[ elementPtr ];
                elementPtr++;
                globalIdx++;
        }
        outVector[ reduceMap[ threadIdx.x] ] += temp[ threadIdx.x ];
    }
}
#endif

#ifdef HAVE_CUDA
template< typename Real,
          typename Index,
          typename InVector,
          typename OutVector >
__global__
void AdEllpackVectorProductCuda2( const AdEllpack< Real, Devices::Cuda, Index >* matrix,
                                           const InVector* inVector,
                                           OutVector* outVector,
                                           const int gridIdx )
{
    matrix->spmvCuda2( *inVector, *outVector, gridIdx );
}

template< typename Real,
          typename Index,
          typename InVector,
          typename OutVector >
__global__
void AdEllpackVectorProductCuda4( const AdEllpack< Real, Devices::Cuda, Index >* matrix,
                                           const InVector* inVector,
                                           OutVector* outVector,
                                           const int gridIdx )
{
    matrix->spmvCuda4( *inVector, *outVector, gridIdx );
}

template< typename Real,
          typename Index,
          typename InVector,
          typename OutVector >
__global__
void AdEllpackVectorProductCuda8( const AdEllpack< Real, Devices::Cuda, Index >* matrix,
                                           const InVector* inVector,
                                           OutVector* outVector,
                                           const int gridIdx )
{
    matrix->spmvCuda8( *inVector, *outVector, gridIdx );
}

template< typename Real,
          typename Index,
          typename InVector,
          typename OutVector >
__global__
void AdEllpackVectorProductCuda16( const AdEllpack< Real, Devices::Cuda, Index >* matrix,
                                            const InVector* inVector,
                                            OutVector* outVector,
                                            const int gridIdx )
{
    matrix->spmvCuda16( *inVector, *outVector, gridIdx );
}

template< typename Real,
          typename Index,
          typename InVector,
          typename OutVector >
__global__
void AdEllpackVectorProductCuda32( const AdEllpack< Real, Devices::Cuda, Index >* matrix,
                                            const InVector* inVector,
                                            OutVector* outVector,
                                            const int gridIdx )
{
    matrix->spmvCuda32( *inVector, *outVector, gridIdx );
}
#endif

#ifdef HAVE_CUDA
template<>
class AdEllpackDeviceDependentCode< Devices::Cuda >
{
public:

    typedef Devices::Cuda Device;

    template< typename Real,
              typename Index,
              typename InVector,
              typename OutVector >
    static void vectorProduct( const AdEllpack< Real, Device, Index >& matrix,
                               const InVector& inVector,
                               OutVector& outVector )
    {
        typedef AdEllpack< Real, Devices::Cuda, Index > Matrix;
	typedef typename Matrix::IndexType IndexType;
	Matrix* kernel_this = Devices::Cuda::passToDevice( matrix );
	InVector* kernel_inVector = Devices::Cuda::passToDevice( inVector );
	OutVector* kernel_outVector = Devices::Cuda::passToDevice( outVector );
	if( matrix.totalLoad < 2 )
	{
	    dim3 blockSize( 256 ), cudaGridSize( Devices::Cuda::getMaxGridSize() );
	    IndexType cudaBlocks = roundUpDivision( matrix.reduceMap.getSize(), blockSize.x );
	    IndexType cudaGrids = roundUpDivision( cudaBlocks, Devices::Cuda::getMaxGridSize() );
	    for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
	    {
	        if( gridIdx == cudaGrids - 1 )
		    cudaGridSize.x = cudaBlocks % Devices::Cuda::getMaxGridSize();
	        const int sharedMemory = blockSize.x * sizeof( Real );
	        AdEllpackVectorProductCuda2< Real, Index, InVector, OutVector >
                                                    <<< cudaGridSize, blockSize, sharedMemory >>>
                                                    ( kernel_this,
                                                      kernel_inVector,
                                                      kernel_outVector,
                                                      gridIdx );
	    }
	    TNL_CHECK_CUDA_DEVICE;
	    Devices::Cuda::freeFromDevice( kernel_this );
	    Devices::Cuda::freeFromDevice( kernel_inVector );
	    Devices::Cuda::freeFromDevice( kernel_outVector );
	    TNL_CHECK_CUDA_DEVICE;
	}
	else if( matrix.totalLoad < 4 )
	{
	    dim3 blockSize( 192 ), cudaGridSize( Devices::Cuda::getMaxGridSize() );
	    IndexType cudaBlocks = roundUpDivision( matrix.reduceMap.getSize(), blockSize.x );
	    IndexType cudaGrids = roundUpDivision( cudaBlocks, Devices::Cuda::getMaxGridSize() );
	    for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
	    {
	        if( gridIdx == cudaGrids - 1 )
		    cudaGridSize.x = cudaBlocks % Devices::Cuda::getMaxGridSize();
	        const int sharedMemory = blockSize.x * sizeof( Real );
	        AdEllpackVectorProductCuda4< Real, Index, InVector, OutVector >
                                                    <<< cudaGridSize, blockSize, sharedMemory >>>
                                                    ( kernel_this,
                                                      kernel_inVector,
                                                      kernel_outVector,
                                                      gridIdx );
	    }
	    TNL_CHECK_CUDA_DEVICE;
	    Devices::Cuda::freeFromDevice( kernel_this );
	    Devices::Cuda::freeFromDevice( kernel_inVector );
	    Devices::Cuda::freeFromDevice( kernel_outVector );
	    TNL_CHECK_CUDA_DEVICE;
	}
	else if( matrix.totalLoad < 8 )
	{
	    dim3 blockSize( 128 ), cudaGridSize( Devices::Cuda::getMaxGridSize() );
	    IndexType cudaBlocks = roundUpDivision( matrix.reduceMap.getSize(), blockSize.x );
	    IndexType cudaGrids = roundUpDivision( cudaBlocks, Devices::Cuda::getMaxGridSize() );
	    for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
	    {
	        if( gridIdx == cudaGrids - 1 )
		    cudaGridSize.x = cudaBlocks % Devices::Cuda::getMaxGridSize();
	        const int sharedMemory = blockSize.x * sizeof( Real );
	        AdEllpackVectorProductCuda8< Real, Index, InVector, OutVector >
                                                    <<< cudaGridSize, blockSize, sharedMemory >>>
                                                    ( kernel_this,
                                                      kernel_inVector,
                                                      kernel_outVector,
                                                      gridIdx );
	    }
	    TNL_CHECK_CUDA_DEVICE;
	    Devices::Cuda::freeFromDevice( kernel_this );
	    Devices::Cuda::freeFromDevice( kernel_inVector );
	    Devices::Cuda::freeFromDevice( kernel_outVector );
	    TNL_CHECK_CUDA_DEVICE;
	}
	else if( matrix.totalLoad < 16 )
	{
	    dim3 blockSize( 128 ), cudaGridSize( Devices::Cuda::getMaxGridSize() );
	    IndexType cudaBlocks = roundUpDivision( matrix.reduceMap.getSize(), blockSize.x );
	    IndexType cudaGrids = roundUpDivision( cudaBlocks, Devices::Cuda::getMaxGridSize() );
	    for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
	    {
	        if( gridIdx == cudaGrids - 1 )
		    cudaGridSize.x = cudaBlocks % Devices::Cuda::getMaxGridSize();
	        const int sharedMemory = blockSize.x * sizeof( Real );
	        AdEllpackVectorProductCuda16< Real, Index, InVector, OutVector >
                                                     <<< cudaGridSize, blockSize, sharedMemory >>>
                                                     ( kernel_this,
                                                       kernel_inVector,
                                                       kernel_outVector,
                                                       gridIdx );
	    }
	    TNL_CHECK_CUDA_DEVICE;
	    Devices::Cuda::freeFromDevice( kernel_this );
	    Devices::Cuda::freeFromDevice( kernel_inVector );
	    Devices::Cuda::freeFromDevice( kernel_outVector );
	    TNL_CHECK_CUDA_DEVICE;
	}
	else
	{
	    dim3 blockSize( 96 ), cudaGridSize( Devices::Cuda::getMaxGridSize() );
	    IndexType cudaBlocks = roundUpDivision( matrix.reduceMap.getSize(), blockSize.x );
	    IndexType cudaGrids = roundUpDivision( cudaBlocks, Devices::Cuda::getMaxGridSize() );
	    for( IndexType gridIdx = 0; gridIdx < cudaGrids; gridIdx++ )
	    {
	        if( gridIdx == cudaGrids - 1 )
		    cudaGridSize.x = cudaBlocks % Devices::Cuda::getMaxGridSize();
	        const int sharedMemory = blockSize.x * sizeof( Real );
	        AdEllpackVectorProductCuda32< Real, Index, InVector, OutVector >
                                                     <<< cudaGridSize, blockSize, sharedMemory >>>
                                                     ( kernel_this,
                                                       kernel_inVector,
                                                       kernel_outVector,
                                                       gridIdx );
	    }
	    TNL_CHECK_CUDA_DEVICE;
	    Devices::Cuda::freeFromDevice( kernel_this );
	    Devices::Cuda::freeFromDevice( kernel_inVector );
	    Devices::Cuda::freeFromDevice( kernel_outVector );
	    TNL_CHECK_CUDA_DEVICE;
	}
    }

};

#endif

} // namespace Matrices
} // namespace TNL
