#include <SparseMatrixFwd.hpp>

Collaboration diagram for Nektar::StorageSmvBsr< T >:

Classes
class	const_iterator

Public Types
typedef T	DataType

typedef Array< OneD, DataType >	DataVectorType

typedef Array< OneD, const DataType >	ConstDataVectorType

typedef Array< OneD, IndexType >	IndexVectorType

typedef void(*	MultiplyKernel )(const double , const double , double *)

Public Member Functions
	StorageSmvBsr (const IndexType blkRows, const IndexType blkCols, const IndexType blkDim, const BCOMatType &bcoMat, const MatrixStorage matType=eFULL)

	StorageSmvBsr (const StorageSmvBsr &src)

	~StorageSmvBsr ()

const IndexType	GetRows () const

const IndexType	GetColumns () const

const IndexType	GetNumNonZeroEntries () const

const IndexType	GetNumStoredDoubles () const

const IndexType	GetBlkSize () const

const DataType	GetFillInRatio () const

const size_t	GetMemoryUsage (IndexType nnz, IndexType nRows) const

const_iterator	begin () const

const_iterator	end () const

const boost::call_traits < DataType >::const_reference	GetValue (IndexType row, IndexType column) const

void	Multiply (const DataType in, DataType out)

void	Multiply (const DataVectorType &in, DataVectorType &out)

void	MultiplyLight (const DataVectorType &in, DataVectorType &out)

Protected Member Functions
void	processBcoInput (const IndexType blkRows, const IndexType blkColumns, const IndexType blkDim, const BCOMatType &bcoMat)

void	Multiply_1x1 (const int mb, const int kb, const double val, const int bindx, const int bpntrb, const int bpntre, const double b, double c)
	Zero-based CSR multiply. Essentially this is slightly modified copy-paste from NIST Sparse Blas 0.9b routine CSR_VecMult_CAB_double() More...

void	Multiply_2x2 (const int mb, const int kb, const double val, const int bindx, const int bpntrb, const int bpntre, const double b, double c)
	Zero-based BSR multiply unrolled for 2x2 blocks. Essentially this is slightly optimised copy-paste from NIST Sparse Blas 0.9b routine BSR_VecMult_BAB_double() More...

void	Multiply_3x3 (const int mb, const int kb, const double val, const int bindx, const int bpntrb, const int bpntre, const double b, double c)
	Zero-based BSR multiply unrolled for 3x3 blocks. More...

void	Multiply_4x4 (const int mb, const int kb, const double val, const int bindx, const int bpntrb, const int bpntre, const double b, double c)
	Zero-based BSR multiply unrolled for 4x4 blocks. More...

void	Multiply_generic (const int mb, const int kb, const double val, const int bindx, const int bpntrb, const int bpntre, const double b, double c)
	Generic zero-based BSR multiply for higher matrix ranks. More...

Protected Attributes
MultiplyKernel	m_mvKernel

MatrixStorage	m_matType

IndexType	m_blkRows

IndexType	m_blkCols

IndexType	m_blkDim

IndexType	m_bnnz

IndexType	m_nnz

DataVectorType	m_val

IndexVectorType	m_indx

IndexVectorType	m_pntr

Detailed Description

template<typename T>
class Nektar::StorageSmvBsr< T >

Definition at line 71 of file SparseMatrixFwd.hpp.

Member Typedef Documentation

template<typename T >

typedef Array<OneD, const DataType> Nektar::StorageSmvBsr< T >::ConstDataVectorType

Definition at line 76 of file StorageSmvBsr.hpp.

template<typename T >

typedef T Nektar::StorageSmvBsr< T >::DataType

Definition at line 74 of file StorageSmvBsr.hpp.

template<typename T >

typedef Array<OneD, DataType> Nektar::StorageSmvBsr< T >::DataVectorType

Definition at line 75 of file StorageSmvBsr.hpp.

template<typename T >

typedef Array<OneD, IndexType> Nektar::StorageSmvBsr< T >::IndexVectorType

Definition at line 77 of file StorageSmvBsr.hpp.

template<typename T >

typedef void(* Nektar::StorageSmvBsr< T >::MultiplyKernel)(const double *, const double *, double *)

Definition at line 79 of file StorageSmvBsr.hpp.

Constructor & Destructor Documentation

template<typename DataType >

Nektar::StorageSmvBsr< DataType >::StorageSmvBsr	(	const IndexType	blkRows,
		const IndexType	blkCols,
		const IndexType	blkDim,
		const BCOMatType &	bcoMat,
		const MatrixStorage	matType = `eFULL`
	)

Todo:

: - iterators over strictly lower-triangular part

number off-diagonal elements calculation
clear distinction between stored and factual nonzeros (row density)

Definition at line 220 of file StorageSmvBsr.cpp.

References Nektar::eFULL, and Nektar::StorageSmvBsr< T >::processBcoInput().

                                                 :
         m_matType (matType),
         m_blkRows (blkRows),
         m_blkCols (blkCols),
         m_blkDim  (blkDim),
         m_bnnz    (bcoMat.size()),
         m_nnz     (0),
         m_val     (m_bnnz * blkDim*blkDim),
         m_indx    (m_bnnz+1),
         m_pntr    (blkRows+1)
     {
         if (matType != Nektar::eFULL)
         {
             /// \todo: - iterators over strictly lower-triangular part
             ///        - number off-diagonal elements calculation
             ///        - clear distinction between stored and factual nonzeros
             ///          (row density)
             std::cout << "matrix type not implemented" << std::endl;
             throw 1;
         }
 
 #ifdef NEKTAR_USING_SMV
         // Set pointer to rank-specific matrix-vector multiply kernel.
         // Number of ranks is defined by LibSMV library
         switch (blkDim)
         {
             #define BOOST_PP_LOCAL_MACRO(n)   case n: m_mvKernel = Smv::F77NAME(smv_##n); break;
             #define BOOST_PP_LOCAL_LIMITS     (1, LIBSMV_MAX_RANK)
             #include BOOST_PP_LOCAL_ITERATE()
         }
 #endif
 
         processBcoInput(blkRows,blkCols,blkDim,bcoMat);
     }

template<typename DataType >

Nektar::StorageSmvBsr< DataType >::StorageSmvBsr ( const StorageSmvBsr< T > & src )

Definition at line 262 of file StorageSmvBsr.cpp.

                                                                   :
         m_matType(src.m_matType),
         m_blkRows (src.m_blkRows),
         m_blkCols (src.m_blkCols),
         m_blkDim(src.m_blkDim),
         m_bnnz(src.m_bnnz),
         m_nnz(src.m_nnz),
         m_val(src.m_val),
         m_indx(src.m_indx),
         m_pntr(src.m_pntr)
     {
     }

template<typename DataType >

Nektar::StorageSmvBsr< DataType >::~StorageSmvBsr ( )

Definition at line 276 of file StorageSmvBsr.cpp.

277 {

278 }

Member Function Documentation

template<typename DataType >

StorageSmvBsr< DataType >::const_iterator Nektar::StorageSmvBsr< DataType >::begin ( ) const

Definition at line 358 of file StorageSmvBsr.cpp.

Referenced by Nektar::StorageSmvBsr< T >::const_iterator::const_iterator().

     {
         return const_iterator(m_matType, 0, m_nnz, m_blkDim, m_val, m_indx, m_pntr);
     }

template<typename DataType >

StorageSmvBsr< DataType >::const_iterator Nektar::StorageSmvBsr< DataType >::end ( ) const

Definition at line 364 of file StorageSmvBsr.cpp.

     {
         return const_iterator(m_matType, m_nnz, m_nnz, m_blkDim, m_val, m_indx, m_pntr);
     }

template<typename DataType >

const IndexType Nektar::StorageSmvBsr< DataType >::GetBlkSize ( ) const

Definition at line 300 of file StorageSmvBsr.cpp.

     {
         return m_blkDim;
     }

template<typename DataType >

const IndexType Nektar::StorageSmvBsr< DataType >::GetColumns ( ) const

Definition at line 288 of file StorageSmvBsr.cpp.

     {
         return m_blkCols*m_blkDim;
     }

template<typename DataType >

const DataType Nektar::StorageSmvBsr< DataType >::GetFillInRatio ( ) const

Definition at line 313 of file StorageSmvBsr.cpp.

     {
         return (DataType)(m_bnnz*m_blkDim*m_blkDim)/(DataType)m_nnz;
     }

template<typename DataType >

const size_t Nektar::StorageSmvBsr< DataType >::GetMemoryUsage	(	IndexType	nnz,
		IndexType	nRows
	)		const

Definition at line 320 of file StorageSmvBsr.cpp.

     {
         return sizeof(DataType) *m_val.capacity()   +
                sizeof(IndexType)*m_indx.capacity() +
                sizeof(IndexType)*m_pntr.capacity()  +
                sizeof(IndexType)*5   + //< blkRows + blkCols + blkDim + nnz + bnnz
                sizeof(MatrixStorage);
     }

template<typename DataType >

const IndexType Nektar::StorageSmvBsr< DataType >::GetNumNonZeroEntries ( ) const

Definition at line 294 of file StorageSmvBsr.cpp.

     {
         return m_nnz;
     }

template<typename DataType >

const IndexType Nektar::StorageSmvBsr< DataType >::GetNumStoredDoubles ( ) const

Definition at line 307 of file StorageSmvBsr.cpp.

     {
         return m_bnnz*m_blkDim*m_blkDim;
     }

template<typename DataType >

const IndexType Nektar::StorageSmvBsr< DataType >::GetRows ( ) const

Definition at line 282 of file StorageSmvBsr.cpp.

     {
         return m_blkRows*m_blkDim;
     }

template<typename DataType >

const boost::call_traits< DataType >::const_reference Nektar::StorageSmvBsr< DataType >::GetValue	(	IndexType	row,
		IndexType	column
	)		const

Definition at line 331 of file StorageSmvBsr.cpp.

     {
         IndexType  brow = grow    / m_blkDim;
         IndexType  bcol = gcolumn / m_blkDim;
         IndexType  lrow = grow    % m_blkDim;
         IndexType  lcol = gcolumn % m_blkDim;
 
         // rewind to the first entry of the first
         // block in the current block row
         IndexType  offset = m_pntr[brow]*m_blkDim*m_blkDim;
 
         IndexType i;
         static DataType defaultReturnValue;
         for( i = m_pntr[brow]; i < m_pntr[brow+1]; i++)
         {
             if(bcol == m_indx[i])
             {
                 return m_val[offset+lrow + lcol*m_blkDim];
             }
             offset += m_blkDim*m_blkDim;
         }
 
         return defaultReturnValue;
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply	(	const DataType *	in,
		DataType *	out
	)

Definition at line 412 of file StorageSmvBsr.cpp.

     {
         const double* b = &in[0];
               double* c = &out[0];
         const double* val = &m_val[0];
         const int* bindx  = (int*)&m_indx[0];
         const int* bpntrb = (int*)&m_pntr[0];
         const int* bpntre = (int*)&m_pntr[0]+1;
         const int  mb = m_blkRows;
         const int  kb = m_blkCols;
 
         switch(m_blkDim)
         {
         case 1:  Multiply_1x1(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 2:  Multiply_2x2(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #ifndef NEKTAR_USING_SMV
         case 3:  Multiply_3x3(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 4:  Multiply_4x4(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #endif
         default:
 #ifdef NEKTAR_USING_SMV
             if (m_blkDim <= LIBSMV_MAX_RANK)
             {
                 Multiply_libsmv(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
             else
 #endif
             {
                 Multiply_generic(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply	(	const DataVectorType &	in,
		DataVectorType &	out
	)

Definition at line 375 of file StorageSmvBsr.cpp.

     {
         const double* b = &in[0];
               double* c = &out[0];
         const double* val = &m_val[0];
         const int* bindx  = (int*)&m_indx[0];
         const int* bpntrb = (int*)&m_pntr[0];
         const int* bpntre = (int*)&m_pntr[0]+1;
         const int  mb = m_blkRows;
         const int  kb = m_blkCols;
 
         switch(m_blkDim)
         {
         case 1:  Multiply_1x1(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 2:  Multiply_2x2(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #ifndef NEKTAR_USING_SMV
         case 3:  Multiply_3x3(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 4:  Multiply_4x4(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #endif
         default:
 #ifdef NEKTAR_USING_SMV
             if (m_blkDim <= LIBSMV_MAX_RANK)
             {
                 Multiply_libsmv(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
             else
 #endif
             {
                 Multiply_generic(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_1x1	(	const int	mb,
		const int	kb,
		const double *	val,
		const int *	bindx,
		const int *	bpntrb,
		const int *	bpntre,
		const double *	b,
		double *	c
	)

protected

Zero-based CSR multiply. Essentially this is slightly modified copy-paste from NIST Sparse Blas 0.9b routine CSR_VecMult_CAB_double()

Definition at line 490 of file StorageSmvBsr.cpp.

     {
         for (int i=0;i!=mb;i++)
         {
             double t = 0;
             int jb = bpntrb[i];
             int je = bpntre[i];
             for (int j=jb;j!=je;j++)
             {
                 t += b[bindx[j]] * (*val++);
             }
             c[i] = t;
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_2x2	(	const int	mb,
		const int	kb,
		const double *	val,
		const int *	bindx,
		const int *	bpntrb,
		const int *	bpntre,
		const double *	b,
		double *	c
	)

protected

Zero-based BSR multiply unrolled for 2x2 blocks. Essentially this is slightly optimised copy-paste from NIST Sparse Blas 0.9b routine BSR_VecMult_BAB_double()

Definition at line 517 of file StorageSmvBsr.cpp.

     {
         const int lb = 2;
 
         const double *pval = val;
         double *pc=c;
 
         for (int i=0;i!=mb;i++)
         {
             int jb = bpntrb[i];
             int je = bpntre[i];
             pc[0] = 0;
             pc[1] = 0;
             for (int j=jb;j!=je;j++)
             {
                 int bs=bindx[j]*lb;
                 const double *pb = &b[bs];
 
                 pc[0] += pb[0] * pval[0] + pb[1] * pval[2];
                 pc[1] += pb[0] * pval[1] + pb[1] * pval[3];
                 pval += 4;
             }
             pc += 2;
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_3x3	(	const int	mb,
		const int	kb,
		const double *	val,
		const int *	bindx,
		const int *	bpntrb,
		const int *	bpntre,
		const double *	b,
		double *	c
	)

protected

Zero-based BSR multiply unrolled for 3x3 blocks.

Definition at line 553 of file StorageSmvBsr.cpp.

     {
         const int lb = 3;
 
         const double *pval = val;
         double *pc=c;
 
         for (int i=0;i!=mb;i++)
         {
             int jb = bpntrb[i];
             int je = bpntre[i];
             pc[0] = 0;
             pc[1] = 0;
             pc[2] = 0;
             for (int j=jb;j!=je;j++)
             {
                 int bs=bindx[j]*lb;
                 const double *pb = &b[bs];
 
                 pc[0] += pb[0] * pval[0] + pb[1] * pval[3] + pb[2] * pval[6];
                 pc[1] += pb[0] * pval[1] + pb[1] * pval[4] + pb[2] * pval[7];
                 pc[2] += pb[0] * pval[2] + pb[1] * pval[5] + pb[2] * pval[8];
                 pval += 9;
             }
             pc += 3;
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_4x4	(	const int	mb,
		const int	kb,
		const double *	val,
		const int *	bindx,
		const int *	bpntrb,
		const int *	bpntre,
		const double *	b,
		double *	c
	)

protected

Zero-based BSR multiply unrolled for 4x4 blocks.

Definition at line 591 of file StorageSmvBsr.cpp.

     {
         const int lb = 4;
 
         const double *pval = val;
         double *pc=c;
 
         for (int i=0;i!=mb;i++)
         {
             int jb = bpntrb[i];
             int je = bpntre[i];
             pc[0] = 0;
             pc[1] = 0;
             pc[2] = 0;
             pc[3] = 0;
             for (int j=jb;j!=je;j++)
             {
                 int bs=bindx[j]*lb;
                 const double *pb = &b[bs];
 
                 pc[0] += pb[0] * pval[0] + pb[1] * pval[4] + pb[2] * pval[ 8] + pb[3] * pval[12];
                 pc[1] += pb[0] * pval[1] + pb[1] * pval[5] + pb[2] * pval[ 9] + pb[3] * pval[13];
                 pc[2] += pb[0] * pval[2] + pb[1] * pval[6] + pb[2] * pval[10] + pb[3] * pval[14];
                 pc[3] += pb[0] * pval[3] + pb[1] * pval[7] + pb[2] * pval[11] + pb[3] * pval[15];
                 pval += 16;
             }
             pc += 4;
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_generic	(	const int	mb,
		const int	kb,
		const double *	val,
		const int *	bindx,
		const int *	bpntrb,
		const int *	bpntre,
		const double *	b,
		double *	c
	)

protected

Generic zero-based BSR multiply for higher matrix ranks.

Definition at line 667 of file StorageSmvBsr.cpp.

     {
         const int lb = m_blkDim;
         const double *pval = val;
         const int mm=lb*lb;
         double *pc=c;
         for (int i=0;i!=mb*lb;i++) *pc++ = 0;
 
         pc=c;
         for (int i=0;i!=mb;i++)
         {
             int jb = bpntrb[i];
             int je = bpntre[i];
             for (int j=jb;j!=je;j++)
             {
                 Blas::Dgemv('N',lb,lb,1.0,pval,lb,&b[bindx[j]*lb],1,1.0,pc,1);
                 pval+=mm;
             }
             pc += lb;
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::MultiplyLight	(	const DataVectorType &	in,
		DataVectorType &	out
	)

Definition at line 451 of file StorageSmvBsr.cpp.

     {
         const double* b = &in[0];
               double* c = &out[0];
         const double* val = &m_val[0];
         const int* bindx  = (int*)&m_indx[0];
         const int* bpntrb = (int*)&m_pntr[0];
         const int* bpntre = (int*)&m_pntr[0]+1;
         const int  mb = m_blkRows;
         const int  kb = m_blkCols;
 
         switch(m_blkDim)
         {
         case 1:  Multiply_1x1(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 2:  Multiply_2x2(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #ifndef NEKTAR_USING_SMV
         case 3:  Multiply_3x3(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
         case 4:  Multiply_4x4(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
 #endif
         default:
 #ifdef NEKTAR_USING_SMV
             if (m_blkDim <= LIBSMV_MAX_RANK)
             {
                 Multiply_libsmv(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
             else
 #endif
             {
                 Multiply_generic(mb,kb,val,bindx,bpntrb,bpntre,b,c); return;
             }
         }
     }

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::processBcoInput	(	const IndexType	blkRows,
		const IndexType	blkColumns,
		const IndexType	blkDim,
		const BCOMatType &	bcoMat
	)

protected

Definition at line 702 of file StorageSmvBsr.cpp.

References CellMLToNektar.pycml::copy(), and Nektar::NekConstants::kNekSparseNonZeroTol.

Referenced by Nektar::StorageSmvBsr< T >::StorageSmvBsr().

     {
         IndexType i;
         BCOMatTypeConstIt  entry;
         IndexType rowcoord;
         IndexType colcoord;
         BCOEntryType       value;
 
         std::vector<IndexType> tmp(blkRows+1,0);
 
         // calculate the number of entries on each row
         // and store the result in tmp
         for(entry = bcoMat.begin(); entry != bcoMat.end(); entry++)
         {
             rowcoord = (entry->first).first;
             tmp[rowcoord]++;
         }
         // Based upon this information, fill the array m_pntr
         // which basically contains the offset of each row's 
         // first entry in the other arrays m_val and m_indx
         m_pntr[0] = 0;
         for(i = 0; i < blkRows; i++)
         {
             m_pntr[i+1] = m_pntr[i] + tmp[i];
         }
 
         // Copy the values of m_pntr into tmp as this will be needed 
         // in the following step
         std::copy(&m_pntr[0],&m_pntr[0]+blkRows+1,&tmp[0]);
 
         // Now, fill in index and value entries.
         for(entry = bcoMat.begin(); entry != bcoMat.end(); entry++)
         {
             rowcoord = (entry->first).first;
             colcoord = (entry->first).second;
             value    =  entry->second;
             int blkSize = blkDim*blkDim;
 
             for (i = 0; i < blkSize; i++)
             {
                 m_val [ blkSize*(tmp[rowcoord]) + i ] = value[i];
                 if (std::abs(value[i]) > NekConstants::kNekSparseNonZeroTol) m_nnz++;
             }
 
             m_indx[ tmp[rowcoord] ] = colcoord;
             tmp[rowcoord]++;
         }
     }