ListSparseMatrix.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __solve_SparseMatrix_h__
#define __solve_SparseMatrix_h__

namespace fe
{
namespace ext
{
/**************************************************************************//**
    @brief Dense array of sparse rows

    SparseMatrix has some similar functionailty to Matrix through
    overloaded non-member operations.

    Only a small number of operations are implemented for SparseMatrix.
    Others may be added as needed.

*//***************************************************************************/
template<class T>
class SparseMatrix
{
    public:
                            SparseMatrix(void);
                            SparseMatrix(U32 rows,U32 columns);
                            SparseMatrix(const SparseMatrix<T>& rhs);
                            ~SparseMatrix(void);

        SparseMatrix<T>&    operator=(const SparseMatrix<T>& rhs);

        T                   operator()(U32 i, U32 j) const;
        T&                  operator()(U32 i, U32 j);

        void                reset(U32 rows,U32 columns);
        void                clear(void);
        U32                 rows(void) const        { return m_rows; }
        U32                 columns(void) const     { return m_columns; }

        void                setTranspose(const SparseMatrix<T>& rhs);
        void                setSum(const SparseMatrix<T> &rhs);
        void                setDifference(const SparseMatrix<T> &rhs);

        void                premultiplyDiagonal(const SparseMatrix<T> &diag,
                                    SparseMatrix<T> &b) const;

        void                scale(const F32 scalar,
                                    SparseMatrix<T>& result) const;

        void                transform(const DenseVector<T> &x,
                                    DenseVector<T> &b) const;
        void                transposeTransform(const DenseVector<T> &x,
                                    DenseVector<T> &b) const;

    protected:
        BWORD               seek(U32 i, U32 j) const;

        U32                 m_rows;
        U32                 m_columns;

    class Entry
    {
        public:
                Entry(U32 column):
                        m_column(column),
                        m_value(T(0))                                       {}
            U32 m_column;
            T   m_value;
    };
    class Row: public List<Entry*>
    {
        public:
                Row(void)   { ListCore::setAutoDestruct(TRUE); }
    };

        Row*                m_pRow;
mutable Row*                m_pCurrent;
mutable ListCore::Context   m_context;
};

template<class T>
inline SparseMatrix<T>::SparseMatrix(void):
    m_rows(1),
    m_columns(1),
    m_pCurrent(NULL)
{
    m_pRow=new Row[m_rows];
}

template<class T>
inline SparseMatrix<T>::SparseMatrix(U32 rows,U32 columns):
    m_rows(rows),
    m_columns(columns),
    m_pCurrent(NULL)
{
    m_pRow=new Row[m_rows];
}

template<class T>
inline SparseMatrix<T>::SparseMatrix(const SparseMatrix<T>& rhs):
    m_rows(rhs.m_rows),
    m_columns(rhs.m_columns),
    m_pRow(NULL),
    m_pCurrent(NULL)
{
    operator=(rhs);
}

template<class T>
inline SparseMatrix<T>::~SparseMatrix(void)
{
    delete[] m_pRow;
}

template<class T>
inline void SparseMatrix<T>::reset(U32 rows,U32 columns)
{
    delete[] m_pRow;
    m_rows=rows;
    m_columns=columns;
    m_pCurrent=NULL;
    m_pRow=new Row[m_rows];
}
template<class T>
inline void SparseMatrix<T>::clear(void)
{
    delete[] m_pRow;
    m_pCurrent=NULL;
    m_pRow=new Row[m_rows];
}

template<class T>
inline BWORD SparseMatrix<T>::seek(U32 i, U32 j) const
{
#if FE_BOUNDSCHECK
    if(i >= m_rows || j >= m_columns) { feX(e_invalidRange); }
#endif

    Row& row=m_pRow[i];
    if(m_pCurrent!=&row)
    {
        m_pCurrent=&row;
        row.toHead(m_context);
    }

//  feLog("seek(%d,%d)\n",i,j);
    I32 column= -1;
    if(row.isAtHeadNull(m_context))
    {
//      feLog("%d++\n",column);
        row.postIncrement(m_context);
    }
    while(!row.isAtHeadNull(m_context) && !row.isAtTailNull(m_context) &&
            (column=I32(row.current(m_context)->m_column))<I32(j))
    {
//      feLog("%d++\n",column);
        row.postIncrement(m_context);
    }
    while(!row.isAtHeadNull(m_context) && (row.isAtTailNull(m_context) ||
            (column=I32(row.current(m_context)->m_column))>I32(j)))
    {
//      feLog("%d--\n",column);
        row.postDecrement(m_context);
    }

    if(!row.isAtHeadNull(m_context) && row.isAtTailNull(m_context))
    {
        column=I32(row.current(m_context)->m_column);
    }
//  feLog("at %d\n",column);
    return column==I32(j);
}

template<class T>
inline T SparseMatrix<T>::operator()(U32 i, U32 j) const
{
    if(seek(i,j))
    {
        return m_pRow[i].current(m_context)->m_value;
    }

    return T(0);
}

template<class T>
inline T& SparseMatrix<T>::operator()(U32 i, U32 j)
{
    if(seek(i,j))
    {
        return m_pRow[i].current(m_context)->m_value;
    }

    Entry** ppEntry=m_pRow[i].insertAfter(m_context,new Entry(j));
    return (*ppEntry)->m_value;
}

/** Copy existing matrix
    @relates SparseMatrix
    */
template<class T>
inline SparseMatrix<T>& SparseMatrix<T>::operator=(const SparseMatrix<T>& rhs)
{
    if(this != &rhs)
    {
        m_rows=rhs.m_rows;
        m_columns=rhs.m_columns;

        clear();

        for(U32 i = 0; i<m_rows; i++)
        {
            Row& row=rhs.m_pRow[i];
            rhs.m_pCurrent=&row;
            row.toHead(rhs.m_context);

            Entry* pNode;
            while((pNode=row.postIncrement(rhs.m_context)) != NULL)
            {
                //* effectively  this(j,i) = rhs(i,j)
                operator()(i,pNode->m_column)=pNode->m_value;
            }
        }
    }
    return *this;
}

#if FALSE
/** Sets to the transpose of the argument
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::setTranspose(const SparseMatrix<T>& rhs)
{
    //* WARNING This will thrash the output context.  Is there a faster way?

    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=rhs.m_pRow[i];
        rhs.m_pCurrent=&row;
        row.toHead(rhs.m_context);

        Entry* pNode;
        while((pNode=row.postIncrement(rhs.m_context)) != NULL)
        {
            //* effectively  this(j,i) = rhs(i,j)
            operator()(pNode->m_column,i)=pNode->m_value;
        }
    }
}
#endif

/** add to SparseMatrix in place
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::setSum(const SparseMatrix<T> &rhs)
{
    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=rhs.m_pRow[i];
        rhs.m_pCurrent=&row;
        row.toHead(rhs.m_context);

        Entry* pNode;
        while((pNode=row.postIncrement(rhs.m_context)) != NULL)
        {
            //* effectively  this(i,j) += rhs(i,j)
            operator()(i,pNode->m_column)+=pNode->m_value;
        }
    }
}

/** subtract from SparseMatrix in place
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::setDifference(const SparseMatrix<T> &rhs)
{
    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=rhs.m_pRow[i];
        rhs.m_pCurrent=&row;
        row.toHead(rhs.m_context);

        Entry* pNode;
        while((pNode=row.postIncrement(rhs.m_context)) != NULL)
        {
            //* effectively  this(i,j) -= rhs(i,j)
            operator()(i,pNode->m_column)-=pNode->m_value;
        }
    }
}

/** SparseMatrix-Scalar multiply
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::scale(const F32 scalar,SparseMatrix<T> &result) const
{
    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=m_pRow[i];
        m_pCurrent=&row;
        row.toHead(m_context);

        Entry* pNode;
        while((pNode=row.postIncrement(m_context)) != NULL)
        {
            //* effectively  result(i,j) = A(i,j) * scalar
            result(i,pNode->m_column)=pNode->m_value*scalar;
        }
    }
}

/** SparseMatrix-Vector multiply
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::transform(const DenseVector<T> &x,DenseVector<T> &b) const
{
    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=m_pRow[i];
        m_pCurrent=&row;
        row.toHead(m_context);

        T sum=T(0);
        Entry* pNode;
        while((pNode=row.postIncrement(m_context)) != NULL)
        {
            //* effectively  b[i] += A(i,j) * x[j]
            sum+=pNode->m_value * x[pNode->m_column];
        }
        setAt(b,i,sum);
    }
}

/** SparseMatrix-Vector multiply with matrix transposed
    @relates SparseMatrix
    */
template<class T>
void SparseMatrix<T>::transposeTransform(const DenseVector<T> &x,
        DenseVector<T> &b) const
{
    set(b);

    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=m_pRow[i];
        m_pCurrent=&row;
        row.toHead(m_context);

        Entry* pNode;
        while((pNode=row.postIncrement(m_context)) != NULL)
        {
            //* effectively  b[j] += A(i,j) * x[i]
            b[pNode->m_column]+=pNode->m_value * x[i];
        }
    }
}

/** SparseMatrix-Diagonal multiply
    @relates SparseMatrix

    This amounts to a scale of each row by the corresponding diagonal element.
    */
template<class T>
void SparseMatrix<T>::premultiplyDiagonal(const SparseMatrix<T> &diag,
        SparseMatrix<T> &result) const
{
    for(U32 i = 0; i<m_rows; i++)
    {
        Row& row=m_pRow[i];
        m_pCurrent=&row;
        row.toHead(m_context);

        T d=diag(i,i);

        Entry* pNode;
        while((pNode=row.postIncrement(m_context)) != NULL)
        {
            //* effectively  result(i,j) = A(i,j) * diag(i,i)
            result(i,pNode->m_column)=pNode->m_value*d;
        }
    }
}

/* non member operations */

template<class T>
inline SparseMatrix<T>& set(SparseMatrix<T> &matrix)
{
    matrix.clear();
    return matrix;
}

/** Return true is matrix is square, otherwise return false.
    @relates SparseMatrix
    */
template<class T>
inline BWORD isSquare(const SparseMatrix<T> &matrix)
{
    return(matrix.rows()==matrix.columns());
}

/** Set matrix to identity matrix
    @relates SparseMatrix
    */
template<class T>
inline SparseMatrix<T> &setIdentity(SparseMatrix<T> &matrix)
{
    U32 limit=minimum(matrix.columns(),matrix.rows());
    matrix.clear();
    for(U32 i = 0;i < limit;i++)
        matrix(i,i) = T(1);
    return matrix;
}

/** Return transpose of matrix
    @relates SparseMatrix
    */
template<class T>
inline SparseMatrix<T> transpose(const SparseMatrix<T> &matrix)
{
    SparseMatrix<T> A(matrix.rows(),matrix.columns());
    A.setTranspose(matrix);

    return A;
}

/** SparseMatrix add
    @relates SparseMatrix
    */
template<class T>
SparseMatrix<T> operator+(const SparseMatrix<T> &lhs,
        const SparseMatrix<T> &rhs)
{
    SparseMatrix<T> result=lhs;
    result.setSum(rhs);
    return result;
}

/** SparseMatrix subtract
    @relates SparseMatrix
    */
template<class T>
SparseMatrix<T> operator-(const SparseMatrix<T> &lhs,
        const SparseMatrix<T> &rhs)
{
    SparseMatrix<T> result=lhs;
    result.setDifference(rhs);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where first matrix is a diagonal
    @relates SparseMatrix
    */
template<class T>
SparseMatrix<T> premultiplyDiagonal(const SparseMatrix<T> &lhs,
        const SparseMatrix<T> &rhs)
{
    SparseMatrix<T> result(rhs.rows(),rhs.columns());
    rhs.premultiplyDiagonal(lhs,result);
    return result;
}

/** SparseMatrix-Vector multiply
    @relates SparseMatrix
    */
template<class T>
DenseVector<T> operator*(const SparseMatrix<T> &lhs, const DenseVector<T> &rhs)
{
    DenseVector<T> b(rhs.size());
    lhs.transform(rhs,b);
    return b;
}

/** SparseMatrix-Vector multiply with matrix transposed
    @relates SparseMatrix
    */
template<class T>
DenseVector<T> transposeMultiply(const SparseMatrix<T> &lhs,
        const DenseVector<T>& rhs)
{
    DenseVector<T> b(rhs.size());
    lhs.transposeTransform(rhs,b);
    return b;
}

/** SparseMatrix-Scalar post-multiply
    @relates SparseMatrix
    */
template<class T>
SparseMatrix<T> operator*(const SparseMatrix<T> &matrix, const F32 scalar)
{
    SparseMatrix<T> result(matrix.rows(),matrix.columns());
    matrix.scale(scalar,result);
    return result;
}

/** SparseMatrix-Scalar pre-multiply
    @relates SparseMatrix
    */
template<class T>
SparseMatrix<T> operator*(const F32 scalar, const SparseMatrix<T> &matrix)
{
    SparseMatrix<T> result(matrix.rows(),matrix.columns());
    matrix.scale(scalar,result);
    return result;
}

/** SparseMatrix print
    @relates SparseMatrix

    The output isn't sparse.
    */
template<class T>
String print(const SparseMatrix<T> &matrix)
{
    String s;
    for(U32 i = 0; i<matrix.rows(); i++)
    {
        for(U32 j = 0; j<matrix.columns(); j++)
        {
            s.catf("%6.3f ", matrix(i,j));
        }
        s.cat("\n");
    }
    return s;
}

} /* namespace ext */
} /* namespace fe */

#endif // __solve_SparseMatrix_h__