#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 ()

IndexType	GetRows () const

IndexType	GetColumns () const

IndexType	GetNumNonZeroEntries () const

IndexType	GetNumStoredDoubles () const

IndexType	GetBlkSize () const

DataType	GetFillInRatio () 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 >

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

Definition at line 300 of file StorageSmvBsr.cpp.

     {
         return m_blkDim;
     }

template<typename DataType >

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

Definition at line 288 of file StorageSmvBsr.cpp.

     {
         return m_blkCols*m_blkDim;
     }

template<typename DataType >

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 >

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 >

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

Definition at line 294 of file StorageSmvBsr.cpp.

     {
         return m_nnz;
     }

template<typename DataType >

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

Definition at line 307 of file StorageSmvBsr.cpp.

     {
         return m_bnnz*m_blkDim*m_blkDim;
     }

template<typename DataType >

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]++;
         }
     }