#include <StorageSmvBsr.hpp>

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

const_iterator	begin () const

const_iterator	end () const

const DataType &	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 blkDim, const BCOMatType &bcoMat)

void	Multiply_1x1 (const int mb, 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 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 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 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 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 62 of file StorageSmvBsr.hpp.

Member Typedef Documentation

◆ ConstDataVectorType

template<typename T >

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

Definition at line 68 of file StorageSmvBsr.hpp.

◆ DataType

template<typename T >

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

Definition at line 66 of file StorageSmvBsr.hpp.

◆ DataVectorType

template<typename T >

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

Definition at line 67 of file StorageSmvBsr.hpp.

◆ IndexVectorType

template<typename T >

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

Definition at line 69 of file StorageSmvBsr.hpp.

◆ MultiplyKernel

template<typename T >

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

Definition at line 71 of file StorageSmvBsr.hpp.

Constructor & Destructor Documentation

◆ StorageSmvBsr() [1/2]

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 206 of file StorageSmvBsr.cpp.

    : 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;
    }
 
    processBcoInput(blkRows, blkDim, bcoMat);
}

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

◆ StorageSmvBsr() [2/2]

template<typename DataType >

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

Definition at line 229 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)
{
}

◆ ~StorageSmvBsr()

template<typename DataType >

Nektar::StorageSmvBsr< DataType >::~StorageSmvBsr

Definition at line 236 of file StorageSmvBsr.cpp.

237{

238}

Member Function Documentation

◆ begin()

template<typename DataType >

StorageSmvBsr< DataType >::const_iterator Nektar::StorageSmvBsr< DataType >::begin

Definition at line 314 of file StorageSmvBsr.cpp.

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

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

◆ end()

template<typename DataType >

StorageSmvBsr< DataType >::const_iterator Nektar::StorageSmvBsr< DataType >::end

Definition at line 320 of file StorageSmvBsr.cpp.

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

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

◆ GetBlkSize()

template<typename DataType >

IndexType Nektar::StorageSmvBsr< DataType >::GetBlkSize

Definition at line 258 of file StorageSmvBsr.cpp.

{
    return m_blkDim;
}

◆ GetColumns()

template<typename DataType >

IndexType Nektar::StorageSmvBsr< DataType >::GetColumns

Definition at line 246 of file StorageSmvBsr.cpp.

{
    return m_blkCols * m_blkDim;
}

◆ GetFillInRatio()

template<typename DataType >

DataType Nektar::StorageSmvBsr< DataType >::GetFillInRatio

Definition at line 270 of file StorageSmvBsr.cpp.

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

◆ GetMemoryUsage()

template<typename DataType >

size_t Nektar::StorageSmvBsr< DataType >::GetMemoryUsage

Definition at line 276 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);
}

◆ GetNumNonZeroEntries()

template<typename DataType >

IndexType Nektar::StorageSmvBsr< DataType >::GetNumNonZeroEntries

Definition at line 252 of file StorageSmvBsr.cpp.

{
    return m_nnz;
}

◆ GetNumStoredDoubles()

template<typename DataType >

IndexType Nektar::StorageSmvBsr< DataType >::GetNumStoredDoubles

Definition at line 264 of file StorageSmvBsr.cpp.

{
    return m_bnnz * m_blkDim * m_blkDim;
}

◆ GetRows()

template<typename DataType >

IndexType Nektar::StorageSmvBsr< DataType >::GetRows

Definition at line 240 of file StorageSmvBsr.cpp.

{
    return m_blkRows * m_blkDim;
}

◆ GetValue()

template<typename DataType >

const DataType & Nektar::StorageSmvBsr< DataType >::GetValue	(	IndexType	row,
		IndexType	column
	)		const

Definition at line 286 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;
}

◆ Multiply() [1/2]

template<typename DataType >

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

Definition at line 359 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;
 
    switch (m_blkDim)
    {
        case 1:
            Multiply_1x1(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        case 2:
            Multiply_2x2(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        default:
        {
            Multiply_generic(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        }
    }
}

◆ Multiply() [2/2]

template<typename DataType >

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

Definition at line 331 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;
 
    switch (m_blkDim)
    {
        case 1:
            Multiply_1x1(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        case 2:
            Multiply_2x2(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        default:
        {
            Multiply_generic(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        }
    }
}

◆ Multiply_1x1()

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_1x1	(	const int	mb,
		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 419 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;
    }
}

◆ Multiply_2x2()

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_2x2	(	const int	mb,
		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 441 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;
    }
}

◆ Multiply_3x3()

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_3x3	(	const int	mb,
		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 472 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;
    }
}

◆ Multiply_4x4()

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_4x4	(	const int	mb,
		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 505 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;
    }
}

◆ Multiply_generic()

template<typename DataType >

void Nektar::StorageSmvBsr< DataType >::Multiply_generic	(	const int	mb,
		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 544 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;
    }
}

References Blas::Dgemv().

◆ MultiplyLight()

template<typename DataType >

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

Definition at line 388 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;
 
    switch (m_blkDim)
    {
        case 1:
            Multiply_1x1(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        case 2:
            Multiply_2x2(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        default:
        {
            Multiply_generic(mb, val, bindx, bpntrb, bpntre, b, c);
            return;
        }
    }
}

◆ processBcoInput()

template<typename DataType >

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

protected

Definition at line 577 of file StorageSmvBsr.cpp.

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