SparseMatrix.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

    @ingroup solve

    SparseMatrix has some similar functionality 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 ROW=SparseArray<T> >
class SparseMatrix
{
    public:
                            SparseMatrix(void);
                            SparseMatrix(U32 rows,U32 columns);
                            SparseMatrix(const SparseMatrix<T,ROW>& rhs);
                            ~SparseMatrix(void);

        SparseMatrix<T,ROW>&    operator=(const SparseMatrix<T,ROW>& 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,ROW>& rhs);
        void                setSum(const SparseMatrix<T,ROW> &rhs);
        void                setDifference(const SparseMatrix<T,ROW> &rhs);
        void                setProduct(
                                    const SparseMatrix<T,ROW> &lhs,
                                    const SparseMatrix<T,ROW> &rhs);

        void                premultiplyDiagonal(
                                    const SparseMatrix<T,ROW> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                premultiplyDiagonal(
                                    const DenseVector<T> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                premultiplyInverseDiagonal(
                                    const SparseMatrix<T,ROW> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                premultiplyInverseDiagonal(
                                    const DenseVector<T> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                postmultiplyDiagonal(
                                    const DenseVector<T> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                postmultiplyInverseDiagonal(
                                    const SparseMatrix<T,ROW> &diag,
                                    SparseMatrix<T,ROW> &b) const;
        void                postmultiplyInverseDiagonal(
                                    const DenseVector<T> &diag,
                                    SparseMatrix<T,ROW> &b) const;

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

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

    protected:
        U32                 m_rows;
        U32                 m_columns;

        ROW*                m_pRow;
};

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

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

template<class T,class ROW>
inline SparseMatrix<T,ROW>::SparseMatrix(const SparseMatrix<T,ROW>& rhs):
    m_rows(rhs.m_rows),
    m_columns(rhs.m_columns)
{
    m_pRow=new ROW[m_rows];
    operator=(rhs);
}

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

template<class T,class ROW>
inline void SparseMatrix<T,ROW>::reset(U32 rows,U32 columns)
{
    delete[] m_pRow;
    m_rows=rows;
    m_columns=columns;
    m_pRow=new ROW[m_rows];
}
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::clear(void)
{
    for(U32 i=0;i<m_rows;i++)
    {
        m_pRow[i].clear();
    }
}

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

    const ROW& row=m_pRow[i];
    return row[j];
}

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

    return m_pRow[i][j];
}

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

        for(U32 i=0;i<m_rows;i++)
        {
            const ROW& row=rhs.m_pRow[i];
            m_pRow[i]=row;
        }
    }
    return *this;
}

/** set as transpose of argument
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::setTranspose(const SparseMatrix<T,ROW> &rhs)
{
    if(this != &rhs)
    {
        if(m_rows==rhs.m_rows)
        {
            m_columns=rhs.m_columns;
            clear();
        }
        else
        {
            reset(rhs.m_rows,rhs.m_columns);
        }

        for(U32 i=0;i<m_rows;i++)
        {
            const ROW& row=rhs.m_pRow[i];
            U32 entries=row.entries();
            for(U32 loc=0; loc<entries; loc++)
            {
                //* effectively  this(j,i) = rhs(i,j)
                operator()(row.index(loc),i)=row.entry(loc);
            }
        }
    }
}

/** add to SparseMatrix in place
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::setSum(const SparseMatrix<T,ROW> &rhs)
{
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=rhs.m_pRow[i];
        U32 entries=row.entries();
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  this(i,j) += rhs(i,j)
            operator()(i,row.index(loc))+=row.entry(loc);
        }
    }
}

/** subtract from SparseMatrix in place
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::setDifference(const SparseMatrix<T,ROW> &rhs)
{
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=rhs.m_pRow[i];
        U32 entries=row.entries();
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  this(i,j) -= rhs(i,j)
            operator()(i,row.index(loc))-=row.entry(loc);
        }
    }
}

/** SparseMatrix-SparseMatrix multiply
    @relates SparseMatrix

    Very slow.  For testing only.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::setProduct(
        const SparseMatrix<T,ROW> &lhs,
        const SparseMatrix<T,ROW> &rhs)
{
    FEASSERT(lhs.m_columns==rhs.m_rows);

    if(m_rows==lhs.m_rows)
    {
        m_columns=rhs.m_columns;
        clear();
    }
    else
    {
        reset(lhs.m_rows,rhs.m_columns);
    }

    for(U32 i=0; i<m_rows; i++)
    {
        for(U32 j=0; j<m_columns; j++)
        {
            T sum=T(0);
            for(U32 k=0; k<lhs.m_columns; k++)
            {
                sum += lhs(i,k)*rhs(k,j);
            }
            if(sum!=0.0)
            {
                operator()(i,j)=sum;
            }
        }
    }
}

/** SparseMatrix-Scalar multiply in place
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::scale(const F32 scalar)
{
    for(U32 i=0; i<m_rows; i++)
    {
        ROW& row=m_pRow[i];
        U32 entries=row.entries();
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  this(i,j) *= scalar
            row.entry(loc)*=scalar;
        }
    }
}

/** SparseMatrix-Scalar multiply
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::scale(const F32 scalar,
        SparseMatrix<T,ROW> &result) const
{
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * scalar
            result(i,row.index(loc))=row.entry(loc)*scalar;
        }
    }
}

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

        T sum=T(0);
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  b[i] += this(i,j) * x[j]
            sum+=row.entry(loc) * x[row.index(loc)];
        }
        setAt(b,i,sum);
    }
}

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

    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  b[j] += this(i,j) * x[i]
            b[row.index(loc)]+=row.entry(loc) * x[i];
        }
    }
}

/** Diagonal-SparseMatrix multiply
    @relates SparseMatrix

    This amounts to a scale of each row by the corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::premultiplyDiagonal(
        const SparseMatrix<T,ROW> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        T d=diag(i,i);
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(i,i)
            result(i,row.index(loc))=row.entry(loc)*d;
        }
    }
}

/** Diagonal-SparseMatrix multiply
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    This amounts to a scale of each row by the corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::premultiplyDiagonal(
        const DenseVector<T> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        T d=diag[i];
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(i,i)
            result(i,row.index(loc))=row.entry(loc)*d;
        }
    }
}

/** Inverse_Diagonal-SparseMatrix multiply
    @relates SparseMatrix

    This amounts to a scale of each row by the reciprical of the
    corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::premultiplyInverseDiagonal(
        const SparseMatrix<T,ROW> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        T d=T(1)/diag(i,i);
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(i,i)
            result(i,row.index(loc))=row.entry(loc)*d;
        }
    }
}

/** Inverse_Diagonal-SparseMatrix multiply
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    This amounts to a scale of each row by the reciprical of the
    corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::premultiplyInverseDiagonal(
        const DenseVector<T> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        T d=T(1)/diag[i];
        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(i)
            result(i,row.index(loc))=row.entry(loc)*d;
        }
    }
}

/** Diagonal-SparseMatrix multiply
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    This amounts to a scale of each column by the corresponding
    diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::postmultiplyDiagonal(
        const DenseVector<T> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(j,j)
            const U32 j=row.index(loc);
            result(i,j)=row.entry(loc)*diag[j];
        }
    }
}

/** SparseMatrix-Inverse_Diagonal multiply
    @relates SparseMatrix

    This amounts to a scale of each column by the reciprical of the
    corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::postmultiplyInverseDiagonal(
        const SparseMatrix<T,ROW> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    T scale[m_rows];
    for(U32 i=0; i<m_rows; i++)
    {
        scale[i]=T(1)/diag(i,i);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(j,j)
            const U32 j=row.index(loc);
            result(i,row.index(loc))=row.entry(loc)*scale[j];
        }
    }
}

/** SparseMatrix-Inverse_Diagonal multiply
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    This amounts to a scale of each column by the reciprical of the
    corresponding diagonal element.
    */
template<class T,class ROW>
inline void SparseMatrix<T,ROW>::postmultiplyInverseDiagonal(
        const DenseVector<T> &diag,SparseMatrix<T,ROW> &result) const
{
    if(result.rows()!=m_rows || result.columns()!=m_columns)
    {
        result.reset(m_rows,m_columns);
    }
    T scale[m_rows];
    for(U32 i=0; i<m_rows; i++)
    {
        scale[i]=T(1)/diag(i);
    }
    for(U32 i=0; i<m_rows; i++)
    {
        const ROW& row=m_pRow[i];
        U32 entries=row.entries();

        for(U32 loc=0; loc<entries; loc++)
        {
            //* effectively  result(i,j) = this(i,j) * diag(j)
            const U32 j=row.index(loc);
            result(i,row.index(loc))=row.entry(loc)*scale[j];
        }
    }
}

/* non member operations */

/** Return the horizonatal dimension
    @relates SparseMatrix
*/
template<class T,class ROW>
inline U32 width(const SparseMatrix<T,ROW>& matrix)
{
    return matrix.rows();
}

/** Return the vertical dimension
    @relates SparseMatrix
*/
template<class T,class ROW>
inline U32 height(const SparseMatrix<T,ROW>& matrix)
{
    return matrix.columns();
}

/** @brief Compute the per-element product of the vector and the
    diagonal entries

    @relates SparseMatrix

    Only the diagonal elements are used.  The non-diagonal values are
    not read and treated as zero.
*/
template<class T,class ROW>
inline void premultiplyDiagonal(
        DenseVector<T>& result,
        const SparseMatrix<T,ROW>& diagonal,
        const DenseVector<T>& vector)
{
    //* result[i] = vector[i] / diag[i][i]

    const U32 size=vector.size();
    if(result.size()!=size)
    {
        result.reset(size);
    }
    for(U32 m=0;m<size;m++)
    {
        result[m]=vector[m]*diagonal(m,m);
    }
}

/** @brief Compute the per-element product of the vector and the
    inverse diagonal entries

    @relates SparseMatrix

    Only the diagonal elements are used.  The non-diagonal values are
    not read and treated as zero.
*/
template<class T,class ROW>
inline void premultiplyInverseDiagonal(
        DenseVector<T>& result,
        const SparseMatrix<T,ROW>& diagonal,
        const DenseVector<T>& vector)
{
    //* result[i] = vector[i] / diag[i][i]

    const U32 size=vector.size();
    if(result.size()!=size)
    {
        result.reset(size);
    }
    for(U32 m=0;m<size;m++)
    {
        result[m]=vector[m]/diagonal(m,m);
    }
}

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

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

/** Set matrix to identity matrix
    @relates SparseMatrix
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW> &setIdentity(SparseMatrix<T,ROW> &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,class ROW>
inline SparseMatrix<T,ROW> transpose(const SparseMatrix<T,ROW> &matrix)
{
    SparseMatrix<T,ROW> A(matrix.rows(),matrix.columns());
    A.setTranspose(matrix);

    return A;
}

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

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

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

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

/** SparseMatrix-SparseMatrix multiply where first matrix is presumed diagonal
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& premultiplyDiagonal(SparseMatrix<T,ROW> &result,
        const DenseVector<T> &lhs,const SparseMatrix<T,ROW> &rhs)
{
    result.clear();
    rhs.premultiplyDiagonal(lhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where first matrix is presumed diagonal
    @relates SparseMatrix

    Non-diagonal elements of the first matrix are ignored and treated as zero.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& premultiplyDiagonal(SparseMatrix<T,ROW> &result,
        const SparseMatrix<T,ROW> &lhs,const SparseMatrix<T,ROW> &rhs)
{
    result.clear();
    rhs.premultiplyDiagonal(lhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where first matrix is presumed diagonal
    @relates SparseMatrix

    Diagonal elements of the first matrix are inverted during the multiply.
    Non-diagonal elements of the first matrix are ignored and treated as zero.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& premultiplyInverseDiagonal(
        SparseMatrix<T,ROW> &result,
        const SparseMatrix<T,ROW> &lhs,const SparseMatrix<T,ROW> &rhs)
{
    result.clear();
    rhs.premultiplyInverseDiagonal(lhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where first matrix is a diagonal
    represented with a vector

    @relates SparseMatrix

    Diagonal elements of the first matrix are inverted during the multiply.
    Non-diagonal elements of the first matrix are ignored and treated as zero.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& premultiplyInverseDiagonal(
        SparseMatrix<T,ROW> &result,
        const DenseVector<T> &lhs,const SparseMatrix<T,ROW> &rhs)
{
    result.clear();
    rhs.premultiplyInverseDiagonal(lhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where second matrix is a diagonal
    @relates SparseMatrix

    The diagonal matrix is represented as a vector.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& postmultiplyDiagonal(SparseMatrix<T,ROW> &result,
        const SparseMatrix<T,ROW> &lhs,const DenseVector<T> &rhs)
{
    result.clear();
    lhs.postmultiplyDiagonal(rhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where second matrix is presumed diagonal
    @relates SparseMatrix

    Diagonal elements of the first matrix are inverted during the multiply.
    Non-diagonal elements of the first matrix are ignored and treated as zero.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& postmultiplyInverseDiagonal(
        SparseMatrix<T,ROW> &result,
        const SparseMatrix<T,ROW> &lhs,const SparseMatrix<T,ROW> &rhs)
{
    result.clear();
    lhs.postmultiplyInverseDiagonal(rhs,result);
    return result;
}

/** SparseMatrix-SparseMatrix multiply where second matrix is a diagonal
    represented by a vector

    @relates SparseMatrix

    Diagonal elements of the first matrix are inverted during the multiply.
    Non-diagonal elements of the first matrix are ignored and treated as zero.

    The result argument can be used to avoid the expensive allocation of
    a temporary.  The same object is returned by reference for efficient
    operation in a compound expression.
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& postmultiplyInverseDiagonal(
        SparseMatrix<T,ROW> &result,
        const SparseMatrix<T,ROW> &lhs,const DenseVector<T> &rhs)
{
    result.clear();
    lhs.postmultiplyInverseDiagonal(rhs,result);
    return result;
}

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

/** SparseMatrix-Vector multiply with matrix
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void transformVector(const SparseMatrix<T,ROW> &lhs,
        const DenseVector<T>& in, DenseVector<T>& out)
{
    lhs.transform(in,out);
}

/** SparseMatrix-Vector multiply with matrix transposed
    @relates SparseMatrix
    */
template<class T,class ROW>
inline void transposeTransformVector(const SparseMatrix<T,ROW> &lhs,
        const DenseVector<T>& in, DenseVector<T>& out)
{
    lhs.transposeTransform(in,out);
}

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

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

/** SparseMatrix-Scalar multiply in place
    @relates SparseMatrix
    */
template<class T,class ROW>
inline SparseMatrix<T,ROW>& operator*=(SparseMatrix<T,ROW> &lhs,
        const F32 scalar)
{
    lhs.scale(scalar);
    return lhs;
}

// BCRS 1D, patterned after the 3D version

template <typename T>
class FE_DL_EXPORT SparseMatrix1 : public Component
{
    public:
typedef T                   t_real;

virtual ~SparseMatrix1(void) {}
virtual void multiply(std::vector<t_real> &a_b,
            const std::vector<t_real> &a_x) const                   = 0;
virtual void scale(t_real a_scale)                                  = 0;
};

template <typename T>
class FE_DL_EXPORT FullBCRS1 : public SparseMatrix1<T>
{
    public:
typedef T                   t_real;
        FullBCRS1(void);
virtual ~FullBCRS1(void);


virtual void multiply(std::vector<t_real> &a_b,
            const std::vector<t_real> &a_x) const;
virtual void scale(t_real a_scale);

        void                clear(void);
        void                startrow(void);
virtual void                done(void);
        t_real          &insert(unsigned int a_columnIndex,
                                const t_real &a_block);

        void                dump(void) const;

        void                write(const String &a_filename);
        void                read(const String &a_filename);

        void                writeDense(const String &a_filename);
        void                readDense(const String &a_filename, unsigned int a_M, unsigned int a_N);

        void                writeStructure(const String &a_filename) const;

        void                sparse(const sp<FullBCRS1> &a_other);
        void                project(const sp<FullBCRS1> &a_other);

        //void              write(CRS<T> &aCRS) const;
        //void              read(const CRS<T> &aCRS);


    public:
        std::vector<t_real> &blocks(void) { return m_blocks; }
        std::vector<unsigned int>   &colind(void) { return m_colind; }
        std::vector<unsigned int>   &rowptr(void) { return m_rowptr; }

    protected:
        std::vector<t_real>         m_blocks;
        std::vector<unsigned int>       m_colind;
        std::vector<unsigned int>       m_rowptr;
        unsigned int                    m_r;
        t_real                          m_zeroVector;
};

template <typename T>
class FE_DL_EXPORT UpperTriangularBCRS1 : public FullBCRS1<T>
{
    public:
typedef T                   t_real;
        UpperTriangularBCRS1(void){}
virtual ~UpperTriangularBCRS1(void){}

virtual void multiply(std::vector<t_real> &a_b,
            const std::vector<t_real> &a_x) const;

        void incompleteSqrt(void);
        void backSub(std::vector<t_real> &a_x, const std::vector<t_real> &a_b);
        void transposeForeSub(std::vector<t_real> &a_x, const std::vector<t_real> &a_b);
        void backSolve(std::vector<t_real> &a_x, const std::vector<t_real> &a_b);
virtual void done(void);

    private:
        std::vector<unsigned int>   m_ptr; // buffer for transposeForeSub
        std::vector<t_real>         m_y; // buffer for backSolve

};

template <typename T>
inline FullBCRS1<T>::FullBCRS1(void)
{
    m_r = 0;
    m_zeroVector = 0.0;
}

template <typename T>
inline FullBCRS1<T>::~FullBCRS1(void)
{
}

template <typename T>
inline void FullBCRS1<T>::clear(void)
{
    m_blocks.clear();
    m_colind.clear();
    m_rowptr.clear();
    m_r = 0;
}

template <typename T>
inline void FullBCRS1<T>::startrow(void)
{
    m_rowptr.push_back(m_r);
}

template <typename T>
inline void FullBCRS1<T>::done(void)
{
    m_rowptr.push_back(m_r);
}

template <typename T>
inline void UpperTriangularBCRS1<T>::done(void)
{
    m_y.resize(FullBCRS1<T>::m_rowptr.size());
    FullBCRS1<T>::m_rowptr.push_back(FullBCRS1<T>::m_r);
    m_ptr.resize(FullBCRS1<T>::m_rowptr.size());
}

template <typename T>
inline T &FullBCRS1<T>::insert(unsigned int a_columnIndex, const t_real &a_block)
{
    m_blocks.push_back(a_block);
    m_colind.push_back(a_columnIndex);
    m_r++;
    return m_blocks.back();
}

template <typename T>
inline void FullBCRS1<T>::multiply(std::vector<t_real> &a_b, const std::vector<t_real> &a_x) const
{
    unsigned int n = (unsigned int)(m_rowptr.size() - 1);
    for(unsigned int i = 0; i < n; i++)
    {
        a_b[i] = m_zeroVector;
        for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
        {
            a_b[i] += m_blocks[k] * a_x[m_colind[k]];
        }
    }
}

template <typename T>
inline void FullBCRS1<T>::dump(void) const
{
    t_real offdiagonal_sum;
    offdiagonal_sum = 0.0;
    unsigned int n = (unsigned int)(m_rowptr.size() - 1);
    for(unsigned int i = 0; i < n; i++)
    {
        fe_fprintf(stderr, " %d:  ", i);
        for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
        {
fe_fprintf(stderr, "%g ", m_blocks[k]);
            if(m_colind[k] != i) // diagonal
            {
                offdiagonal_sum = offdiagonal_sum + m_blocks[k];
            }
        }
        fe_fprintf(stderr, "\n");
    }

    feLogGroup("BCRS", "off diagonal\n%g\n", offdiagonal_sum);

fe_fprintf(stderr, "--------------------------\n");
    for(unsigned int i = 0; i < n; i++)
    {
        unsigned int j = 0;
        for(unsigned int k = m_rowptr[i]; k < m_rowptr[i+1]; k++)
        {
            while(j < m_colind[k] && j < n)
            {
                fe_fprintf(stderr, "%6.3g ", 0.0);
                j++;
            }
            fe_fprintf(stderr, "%6.3g ", m_blocks[k]);
            j++;
        }
        fe_fprintf(stderr, "\n");
    }
fe_fprintf(stderr, "--------------------------\n");
}

template <typename T>
inline void FullBCRS1<T>::scale(t_real a_scale)
{
    unsigned int n = (unsigned int)m_blocks.size();
    for(unsigned int i = 0; i < n; i++)
    {
        m_blocks[i] *= a_scale;
    }
}

template <typename T>
inline void UpperTriangularBCRS1<T>::multiply(std::vector<t_real> &a_b, const std::vector<t_real> &a_x) const
{
    unsigned int n = (unsigned int)(FullBCRS1<T>::m_rowptr.size() - 1);
    for(unsigned int i = 0; i < n; i++)
    {
        a_b[i] = FullBCRS1<T>::m_zeroVector;
    }
    for(unsigned int i = 0; i < n; i++)
    {
        unsigned int k = FullBCRS1<T>::m_rowptr[i];
        a_b[i] += FullBCRS1<T>::m_blocks[k] * a_x[i];
        for(k=k+1; k < FullBCRS1<T>::m_rowptr[i+1]; k++)
        {
            a_b[i] += FullBCRS1<T>::m_blocks[k] * a_x[FullBCRS1<T>::m_colind[k]];
            a_b[FullBCRS1<T>::m_colind[k]] += FullBCRS1<T>::m_blocks[k] * a_x[i];
        }
    }
}

template <typename T>
inline void FullBCRS1<T>::sparse(const sp<FullBCRS1> &a_other)
{
    m_blocks.clear();
    m_colind.clear();
    m_rowptr.clear();
    m_r = 0;
    unsigned int n = a_other->m_rowptr.size()-1;
    for(unsigned int i = 0; i < n; i++)
    {
        startrow();
        int d = a_other->m_rowptr[i];
        for(int k=d; k < a_other->m_rowptr[i+1]; k++)
        {
            int c = a_other->m_colind[k];
            if(!(0.0 == a_other->m_blocks[k]))
            {
                insert(c, a_other->m_blocks[k]);
            }
        }
    }
    done();
}

template <typename T>
inline void FullBCRS1<T>::project(const sp<FullBCRS1> &a_other)
{
    unsigned int n = a_other->m_rowptr.size()-1;
    for(unsigned int i = 0; i < n; i++)
    {
        int d = a_other->m_rowptr[i];
        int d_self = m_rowptr[i];
        for(int k=d; k < a_other->m_rowptr[i+1]; k++)
        {
            int c = a_other->m_colind[k];
            for(int k_self=d_self; k_self < m_rowptr[i+1]; k_self++)
            {
                int c_self = m_colind[k_self];
                if(c_self == c)
                {
                    m_blocks[k_self] = a_other->m_blocks[k];
                }
            }
        }
    }
}

template <typename T>
inline void UpperTriangularBCRS1<T>::backSub(std::vector<t_real> &a_x, const std::vector<t_real> &a_b)
{
    unsigned int n = (unsigned int)(FullBCRS1<T>::m_rowptr.size()-1);
    //fe::backSub(FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[n-1]], a_x[n-1], a_b[n-1]);
    a_x[n-1] = a_b[n-1]/FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[n-1]];
    t_real tmp;

    for(int i = n-2; i >= 0; i--)
    {
        a_x[i] = a_b[i];
        int d = FullBCRS1<T>::m_rowptr[i];
        for(unsigned int k=d+1; k < FullBCRS1<T>::m_rowptr[i+1]; k++)
        {
            int c = FullBCRS1<T>::m_colind[k];
            tmp = FullBCRS1<T>::m_blocks[k] * a_x[c];
            a_x[i] -= tmp;
        }
        //fe::backSub(FullBCRS1<T>::m_blocks[d], a_x[i], a_x[i]);
        a_x[i] = a_x[i]/FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[i]];
    }
}

template <typename T>
inline void UpperTriangularBCRS1<T>::transposeForeSub(std::vector<t_real> &a_x, const std::vector<t_real> &a_b)
{
    unsigned int n = (unsigned int)(FullBCRS1<T>::m_rowptr.size()-1);
    for(unsigned int i = 0 ; i < FullBCRS1<T>::m_rowptr.size(); i++)
    {
        m_ptr[i] = FullBCRS1<T>::m_rowptr[i];
    }
    //fe::transposeForeSub(FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[0]], a_x[0], a_b[0]);
    a_x[0] = a_b[0] / FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[0]];

    for(unsigned int i = 1; i < n; i++)
    {
        a_x[i] = a_b[i];
        for(unsigned int m = 0; m < i; m++)
        {
            unsigned int c = FullBCRS1<T>::m_colind[m_ptr[m]];
            while(c < i && m_ptr[m] < FullBCRS1<T>::m_rowptr[m+1]-1)
            {
                m_ptr[m]++;
                c = FullBCRS1<T>::m_colind[m_ptr[m]];
            }
            if(c == i)
            {
                //a_x[i] -= fe::transposeMultiply(FullBCRS1<T>::m_blocks[m_ptr[m]], a_x[m]);
                a_x[i] -= FullBCRS1<T>::m_blocks[m_ptr[m]] * a_x[m];
            }
        }
        //fe::transposeForeSub(FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[i]], a_x[i], a_x[i]);
        a_x[i] = a_x[i] / FullBCRS1<T>::m_blocks[FullBCRS1<T>::m_rowptr[i]];
    }
}

template <typename T>
inline void UpperTriangularBCRS1<T>::backSolve(std::vector<t_real> &a_x, const std::vector<t_real> &a_b)
{
    transposeForeSub(m_y, a_b);
    backSub(a_x, m_y);
}

template <typename T>
inline void UpperTriangularBCRS1<T>::incompleteSqrt(void)
{
    unsigned int d;
    unsigned int n = (unsigned int)(FullBCRS1<T>::m_rowptr.size()-1);
    t_real inv_z;
    t_real z;
    for(unsigned int k = 0; k < n-1; k++)
    {
        d = FullBCRS1<T>::m_rowptr[k];
        FullBCRS1<T>::m_blocks[d] = sqrt(FullBCRS1<T>::m_blocks[d]);
        z = FullBCRS1<T>::m_blocks[d];

        inv_z = 1.0 / z;

        for(unsigned int i=d+1; i < FullBCRS1<T>::m_rowptr[k+1]; i++)
        {
            z = FullBCRS1<T>::m_blocks[i] * inv_z;

            FullBCRS1<T>::m_blocks[i] = z;
        }

        for(unsigned int i=d+1; i < FullBCRS1<T>::m_rowptr[k+1]; i++)
        {
            z = FullBCRS1<T>::m_blocks[i];
            unsigned int h = FullBCRS1<T>::m_colind[i];
            unsigned int g = i;

            for(unsigned int j = FullBCRS1<T>::m_rowptr[h]; j < FullBCRS1<T>::m_rowptr[h+1]; j++)
            {
                for(; g < FullBCRS1<T>::m_rowptr[k+1] && FullBCRS1<T>::m_colind[g] <= FullBCRS1<T>::m_colind[j]; g++)
                {
                    if (FullBCRS1<T>::m_colind[g] == FullBCRS1<T>::m_colind[j])
                    {
                        t_real t;
                        t = z * FullBCRS1<T>::m_blocks[g];
                        FullBCRS1<T>::m_blocks[j] = FullBCRS1<T>::m_blocks[j] - t;
                    }
                }
            }
        }
    }
    d = FullBCRS1<T>::m_rowptr[n-1];
    t_real tmp;
    tmp = sqrt(FullBCRS1<T>::m_blocks[d]);
    FullBCRS1<T>::m_blocks[d] = tmp;
}


} /* namespace */

/** SparseMatrix print
    @relates SparseMatrix

    The output isn't sparse.
    */
template<class T,class ROW>
String print(const ext::SparseMatrix<T,ROW> &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 */

#endif // __solve_SparseMatrix_h__