Commit b6f94aa5 authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Adding tnlSlicedEllpackMatrix.

parent 65593480

src/core/mfuncs.h

+5 −0
Original line number Diff line number Diff line
@@ -74,6 +74,11 @@ bool isSmall( const Real& v, 
{

   return ( -tolerance <= v && v <= tolerance );

}



int roundUpDivison( const int num, const int div )

{

   return num / div + ( num % div != 0 );

}

/*template< typename T >

void swap( T& a, T& b)

{

src/implementation/matrices/CMakeLists.txt

+2 −1
Original line number Diff line number Diff line 
SET( headers tnlDenseMatrix_impl.h

             tnlTridiagonalMatrix_impl.h

             tnlMultidiagonalMatrix_impl.h 

             tnlEllpackMatrix_impl.h )

             tnlEllpackMatrix_impl.h

             tnlSlicedEllpackMatrix_impl.h )



SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/implementation/matrices )

set( common_SOURCES

src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h

0 → 100644
+396 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          tnlSlicedSlicedEllpackMatrix_impl.h  -  description

                             -------------------

    begin                : Dec 8, 2013

    copyright            : (C) 2013 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/***************************************************************************

 *                                                                         *

 *   This program is free software; you can redistribute it and/or modify  *

 *   it under the terms of the GNU General Public License as published by  *

 *   the Free Software Foundation; either version 2 of the License, or     *

 *   (at your option) any later version.                                   *

 *                                                                         *

 ***************************************************************************/



#ifndef TNLSLICEDELLPACKMATRIX_IMPL_H_

#define TNLSLICEDELLPACKMATRIX_IMPL_H_



#include <matrices/tnlSlicedEllpackMatrix.h>

#include <core/vectors/tnlVector.h>

#include <core/mfuncs.h>



template< typename Real,

          typename Device,

          typename Index >

tnlSlicedEllpackMatrix< Real, Device, Index > :: tnlSlicedEllpackMatrix()

: rows( 0 ),

  columns( 0 ),

  rowLengths( 0 )

{

};



template< typename Real,

          typename Device,

          typename Index >

tnlString tnlSlicedEllpackMatrix< Real, Device, Index > :: getType()

{

   return tnlString( "tnlSlicedEllpackMatrix< ") +

          tnlString( GetParameterType( Real( 0.0 ) ) ) +

          tnlString( ", " ) +

          Device :: getDeviceType() +

          tnlString( " >" );

}



template< typename Real,

          typename Device,

          typename Index >

tnlString tnlSlicedEllpackMatrix< Real, Device, Index >::getTypeVirtual() const

{

   return this->getType();

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: setDimensions( const IndexType rows,

                                                                     const IndexType columns )

{

   tnlAssert( rows > 0 && columns > 0,

              cerr << "rows = " << rows

                   << " columns = " << columns << endl );

   this->rows = rows;

   this->columns = columns;

   return allocateElements();

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Vector >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: setRowLengths( const Vector& rowLengths )

{

   const IndexType slices = roundUpDivision( this->rows, SliceSize );

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Real2,

             typename Device2,

             typename Index2 >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: setLike( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix )

{

   this->rowLengths = 0;

   this->setDimensions( matrix.getRows(), matrix.getColumns() );

   if( ! this->setConstantRowLengths( matrix.rowLengths ) )

      return false;

   return true;

}



template< typename Real,

          typename Device,

          typename Index >

Index tnlSlicedEllpackMatrix< Real, Device, Index > :: getNumberOfAllocatedElements() const

{

   return this->values.getSize();

}



template< typename Real,

          typename Device,

          typename Index >

void tnlSlicedEllpackMatrix< Real, Device, Index > :: reset()

{

   this->rows = 0;

   this->columns = 0;

   this->rowLengths = 0;

   this->values.reset();

   this->columnIndexes.reset();

}



template< typename Real,

          typename Device,

          typename Index >

Index tnlSlicedEllpackMatrix< Real, Device, Index >::getRows() const

{

   return this->rows;

}



template< typename Real,

          typename Device,

          typename Index >

Index tnlSlicedEllpackMatrix< Real, Device, Index >::getColumns() const

{

   return this->columns;

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Real2,

             typename Device2,

             typename Index2 >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::operator == ( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const

{

   tnlAssert( this->getRows() == matrix.getRows() &&

              this->getColumns() == matrix.getColumns(),

              cerr << "this->getRows() = " << this->getRows()

                   << " matrix.getRows() = " << matrix.getRows()

                   << " this->getColumns() = " << this->getColumns()

                   << " matrix.getColumns() = " << matrix.getColumns()

                   << " this->getName() = " << this->getName()

                   << " matrix.getName() = " << matrix.getName() );

   // TODO: implement this

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Real2,

             typename Device2,

             typename Index2 >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::operator != ( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const

{

   return ! ( ( *this ) == matrix );

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: setElement( const IndexType row,

                                                            const IndexType column,

                                                            const Real& value )

{

   return this->addToElement( row, column, value, 0.0 );

}



template< typename Real,

          typename Device,

          typename Index >

Real tnlSlicedEllpackMatrix< Real, Device, Index >::getElement( const IndexType row,

                                                          const IndexType column ) const

{

   IndexType i( row * this->rowLengths );

   const IndexType rowEnd( i + this->rowLengths );

   while( i < rowEnd && this->columnIndexes[ i ] < column ) i++;

   if( i == rowEnd || this->columnIndexes[ i ] != column )

      return 0.0;

   return this->values.getElement( i );

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: addToElement( const IndexType row,

                                                              const IndexType column,

                                                              const RealType& value,

                                                              const RealType& thisElementMultiplicator )

{

   tnlAssert( row >= 0 && row < this->rows &&

              column >= 0 && column <= this->rows,

              cerr << " row = " << row

                   << " column = " << column

                   << " this->rows = " << this->rows

                   << " this->columns = " << this-> columns );

   IndexType i( row * this->rowLengths );

   const IndexType rowEnd( i + this->rowLengths );

   while( i < rowEnd && this->columnIndexes[ i ] < column ) i++;

   if( i == rowEnd )

      return false;

   if( this->columnIndexes[ i ] == column )

   {

      this->values[ i ] = thisElementMultiplicator * this->values[ i ] + value;

      return true;

   }

   else

      if( this->columnIndexes[ i ] == this->columns )

      {

         this->columnIndexes[ i ] = column;

         this->values[ i ] = value;

         return true;

      }

      else

      {

         IndexType j = rowEnd - 1;

         while( j > i )

         {

            this->columnIndexes[ j ] = this->columnIndexes[ j - 1 ];

            this->values[ j ] = this->values[ j - 1 ];

            j--;

         }

         this->columnIndexes[ i ] = column;

         this->values[ i ] = value;

         return true;

      }

   return false;

}





template< typename Real,

          typename Device,

          typename Index >

   template< typename Vector >

void tnlSlicedEllpackMatrix< Real, Device, Index >::vectorProduct( const Vector& inVector,

                                                                   Vector& outVector ) const

{

   for( Index row = 0; row < this->getRows(); row ++ )

   {

      Real result = 0.0;

      IndexType i( row * this->rowLengths );

      const IndexType rowEnd( i + this->rowLengths );

      while( i < rowEnd && this->columnIndexes[ i ] < this->columns )

      {

         const Index column = this->columnIndexes.getElement( i );

         result += this->values.getElement( i++ ) * inVector.getElement( column );

      }

      outVector.setElement( row, result );

   }

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Real2,

             typename Index2 >

void tnlSlicedEllpackMatrix< Real, Device, Index > :: addMatrix( const tnlSlicedEllpackMatrix< Real2, Device, Index2 >& matrix,

                                                                 const RealType& matrixMultiplicator,

                                                                 const RealType& thisMatrixMultiplicator )

{

   tnlAssert( false, cerr << "TODO: implement" );

   // TODO: implement

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Real2,

             typename Index2 >

void tnlSlicedEllpackMatrix< Real, Device, Index >::getTransposition( const tnlSlicedEllpackMatrix< Real2, Device, Index2 >& matrix,

                                                                      const RealType& matrixMultiplicator )

{

   tnlAssert( false, cerr << "TODO: implement" );

   // TODO: implement

}



template< typename Real,

          typename Device,

          typename Index >

   template< typename Vector >

bool tnlSlicedEllpackMatrix< Real, Device, Index > :: performSORIteration( const Vector& b,

                                                                           const IndexType row,

                                                                           Vector& x,

                                                                           const RealType& omega ) const

{

   tnlAssert( row >=0 && row < this->getRows(),

              cerr << "row = " << row

                   << " this->getRows() = " << this->getRows()

                   << " this->getName() = " << this->getName() << endl );



   RealType diagonalValue( 0.0 );

   RealType sum( 0.0 );



   IndexType i( row * this->rowLengths );

   const IndexType rowEnd( i + this->rowLengths );

   IndexType column;

   while( i < rowEnd && ( column = this->columnIndexes[ i ] ) < this->columns )

   {

      if( column == row )

         diagonalValue = this->values.getElement( i );

      else

         sum += this->values.getElement( row * this->diagonalsShift.getSize() + i ) * x. getElement( column );

      i++;

   }

   if( diagonalValue == ( Real ) 0.0 )

   {

      cerr << "There is zero on the diagonal in " << row << "-th row of thge matrix " << this->getName() << ". I cannot perform SOR iteration." << endl;

      return false;

   }

   x. setElement( row, x[ row ] + omega / diagonalValue * ( b[ row ] - sum ) );

   return true;

}





template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::save( tnlFile& file ) const

{

   if( ! file.write( &this->rows ) ) return false;

   if( ! file.write( &this->columns ) ) return false;

   if( ! file.write( &this->rowLengths ) ) return false;

   if( ! this->values.save( file ) ) return false;

   if( ! this->columnIndexes.save( file ) ) return false;

   return true;

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::load( tnlFile& file )

{

   if( ! file.read( &this->rows ) ) return false;

   if( ! file.read( &this->columns ) ) return false;

   if( ! file.read( &this->rowLengths ) ) return false;

   if( ! this->values.load( file ) ) return false;

   if( ! this->columnIndexes.load( file ) ) return false;

   return true;

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::save( const tnlString& fileName ) const

{

   return tnlObject::save( fileName );

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::load( const tnlString& fileName )

{

   return tnlObject::load( fileName );

}



template< typename Real,

          typename Device,

          typename Index >

void tnlSlicedEllpackMatrix< Real, Device, Index >::print( ostream& str ) const

{

   for( IndexType row = 0; row < this->getRows(); row++ )

   {

      str <<"Row: " << row << " -> ";

      IndexType i( row * this->rowLengths );

      const IndexType rowEnd( i + this->rowLengths );

      while( i < rowEnd && this->columnIndexes[ i ] < this->columns )

      {

         const Index column = this->columnIndexes.getElement( i );

         str << " Col:" << column << "->" << this->values.getElement( i ) << "\t";

         i++;

      }

      str << endl;

   }

}



template< typename Real,

          typename Device,

          typename Index >

bool tnlSlicedEllpackMatrix< Real, Device, Index >::allocateElements()

{

   if( ! this->values.setSize( this->rows * this->rowLengths ) ||

       ! this->columnIndexes.setSize( this->rows * this->rowLengths ) )

      return false;



   /****

    * Setting a column index to this->columns means that the

    * index is undefined.

    */

   this->columnIndexes.setValue( this->columns );

   return true;

}







#endif /* TNLSLICEDELLPACKMATRIX_IMPL_H_ */

src/matrices/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -2,6 +2,7 @@ SET( headers tnlDenseMatrix.h 
             tnlTridiagonalMatrix.h

             tnlMultidiagonalMatrix.h

             tnlEllpackMatrix.h

             tnlSlicedEllpackMatrix.h



             tnlAdaptiveRgCSRMatrix.h

             tnlCSRMatrix.h

src/matrices/tnlSlicedEllpackMatrix.h

0 → 100644
+122 −0
Original line number Diff line number Diff line 
/***************************************************************************

                          tnlSlicedSlicedEllpackMatrix.h  -  description

                             -------------------

    begin                : Dec 8, 2013

    copyright            : (C) 2013 by Tomas Oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

 ***************************************************************************/



/***************************************************************************

 *                                                                         *

 *   This program is free software; you can redistribute it and/or modify  *

 *   it under the terms of the GNU General Public License as published by  *

 *   the Free Software Foundation; either version 2 of the License, or     *

 *   (at your option) any later version.                                   *

 *                                                                         *

 ***************************************************************************/



#ifndef TNLSLICEDELLPACKMATRIX_H_

#define TNLSLICEDELLPACKMATRIX_H_



#include <core/vectors/tnlVector.h>



template< typename Real, typename Device = tnlHost, typename Index = int, int SliceSize = 32 >

class tnlSlicedEllpackMatrix : public tnlObject

{

   public:



   typedef Real RealType;

   typedef Device DeviceType;

   typedef Index IndexType;



   tnlSlicedEllpackMatrix();



   static tnlString getType();



   tnlString getTypeVirtual() const;



   bool setDimensions( const IndexType rows,

                       const IndexType columns );



   template< typename Vector >

   bool setRowLengths( const Vector& rowLengths );



   template< typename Real2, typename Device2, typename Index2 >

   bool setLike( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix );



   IndexType getNumberOfAllocatedElements() const;



   void reset();



   IndexType getRows() const;



   IndexType getColumns() const;



   template< typename Real2, typename Device2, typename Index2 >

   bool operator == ( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const;



   template< typename Real2, typename Device2, typename Index2 >

   bool operator != ( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const;



   bool setElement( const IndexType row,

                    const IndexType column,

                    const RealType& value );



   bool setRow( const IndexType row,

                const IndexType* columnIndexes,

                const RealType* values,

                const IndexType elements );



   RealType getElement( const IndexType row,

                        const IndexType column ) const;



   bool addToElement( const IndexType row,

                      const IndexType column,

                      const RealType& value,

                      const RealType& thisElementMultiplicator = 1.0 );



   template< typename Vector >

   void vectorProduct( const Vector& inVector,

                       Vector& outVector ) const;



   template< typename Real2, typename Index2 >

   void addMatrix( const tnlSlicedEllpackMatrix< Real2, Device, Index2 >& matrix,

                   const RealType& matrixMultiplicator = 1.0,

                   const RealType& thisMatrixMultiplicator = 1.0 );



   template< typename Real2, typename Index2 >

   void getTransposition( const tnlSlicedEllpackMatrix< Real2, Device, Index2 >& matrix,

                          const RealType& matrixMultiplicator = 1.0 );



   template< typename Vector >

   bool performSORIteration( const Vector& b,

                             const IndexType row,

                             Vector& x,

                             const RealType& omega = 1.0 ) const;



   bool save( tnlFile& file ) const;



   bool load( tnlFile& file );



   bool save( const tnlString& fileName ) const;



   bool load( const tnlString& fileName );



   void print( ostream& str ) const;



   protected:



   bool allocateElements();



   IndexType rows, columns;



   tnlVector< Real, Device, Index > values;



   tnlVector< Index, Device, Index > columnIndexes, slicePointers, sliceRowLentghs;



};



#include <implementation/matrices/tnlSlicedEllpackMatrix_impl.h>





#endif /* TNLSLICEDELLPACKMATRIX_H_ */