UT/UT_SparseMatrix.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      David Wong
 *      Side Effects
 *      477 Richmond Street West
 *      Toronto, Ontario
 *      Canada   M5V 3E7
 *      416-504-9876
 *
 * NAME:        Sparse Matrix Class (C++)
 *
 * COMMENTS:
 *
 * Storage format: Co-ordinate (COO). See
 *     Y. Saad. SPARSKIT: A basic toolkit for sparse matrix computations.
 *         Research Institute for Advanced Computer Science, RIACS-90-20,
 *         1990.
 * for details.
 */

#ifndef __UT_SparseMatrix_H__ 
#define __UT_SparseMatrix_H__

#include "UT_API.h"

#include "UT_ThreadedAlgorithm.h"
#include "UT_Vector.h"

#include <SYS/SYS_Types.h>

#include <iterator>

class UT_IStream;

template <typename T> class UT_SparseMatrixRowT;

template <typename T, bool IsPaged>
class UT_API UT_SparseMatrixT
{
    class ut_MatrixCell
    {
    public:
        int myRow;
        int myCol;
        T myValue;

        inline bool operator<(const ut_MatrixCell &o) const
        {
            if (myRow < o.myRow)
                return true;
            if (myRow == o.myRow && myCol < o.myCol)
                return true;
            return false;
        }
    };

    // For paged matrices, store the size of each page
    // We always allocate full pages.
    static const int CELL_PAGESIZE = 1024;
    static const int CELL_PAGEMASK = 1023;
    static const int CELL_PAGEBITS = 10;

    class ut_CellIterator : public std::iterator<random_access_iterator_tag, ut_MatrixCell>
    {
    public:
        ut_CellIterator(const UT_SparseMatrixT<T, IsPaged> *matrix, int idx)
        { myMatrix = matrix; myPos = idx; repage(); }

        ut_CellIterator operator+(ptrdiff_t n) { return ut_CellIterator(myMatrix, myPos+n); }
        ut_CellIterator operator-(ptrdiff_t n) { return ut_CellIterator(myMatrix, myPos-n); }
        ut_CellIterator &operator+=(ptrdiff_t n) { myPos += n; repage(); return *this; }
        ut_CellIterator &operator-=(ptrdiff_t n) { myPos -= n; repage(); return *this; }
        ut_CellIterator &operator++() { myPos++;  myOffset++;  myData++; if (myOffset >= CELL_PAGESIZE) repage(); return *this; }
        ut_CellIterator &operator--() { myPos--;  myOffset--;  myData--; if (myOffset < 0) repage(); return *this; }
        ut_CellIterator operator++(int) { ut_CellIterator result = *this; myPos++;  myOffset++;  myData++; if (myOffset >= CELL_PAGESIZE) repage(); return result; }
        ut_CellIterator operator--(int) { ut_CellIterator result = *this; myPos--;  myOffset--;  myData--; if (myOffset < 0) repage(); return result; }
        int             operator-(ut_CellIterator b) { return myPos - b.myPos; }
        ut_MatrixCell   &operator[](ptrdiff_t idx) { return myMatrix->getCell(idx + myPos); }

        bool            operator<(ut_CellIterator b) { return myPos < b.myPos; }
        bool            operator==(ut_CellIterator b) { return myPos == b.myPos; }
        bool            operator!=(ut_CellIterator b) { return myPos != b.myPos; }

        ut_MatrixCell   &operator*() const { return *myData; }
        ut_MatrixCell   *operator->() const { return myData; }

    protected:

        void            repage()
        { myPage = myPos >> CELL_PAGEBITS; myOffset = myPos & CELL_PAGEMASK;
          myData = 0;
          if (myPage < myMatrix->myCellPages.entries())
              myData = &myMatrix->myCellPages(myPage)[myOffset];
        }

        const UT_SparseMatrixT<T, IsPaged>      *myMatrix;
        ut_MatrixCell                   *myData;
        ptrdiff_t                        myPos;
        int                              myPage, myOffset;
    };

    // For compact compiled matrix cells, specify the number of values to
    // pack into a single cell.
    static const int CELLBITS = 2;
    static const int CELLSIZE = 1 << CELLBITS;
    static const int CELLMASK = CELLSIZE-1;

    class ut_4MatrixCell
    {
    public:
        T myValue[CELLSIZE];
        int myRow;
        int myCol;
    } SYS_ALIGN16;

public:
    UT_SparseMatrixT();

    // creates a valid 'rows' x 'cols' zero matrix
    UT_SparseMatrixT(int rows, int cols);
    UT_SparseMatrixT(const UT_SparseMatrixT<T, IsPaged> &m);
    ~UT_SparseMatrixT();

    int getNumRows() const { return myNumRows; }
    int getNumCols() const { return myNumCols; }

    /// Return the amount of memory used by this array.
    int64        getMemoryUsage() const;
    int64        getIdealMemoryUsage() const;

    /// Ensures at least this number of cells are available, useful
    /// if you can predict your size.
    void         reserve(int numcells);

    /// Shrinks the array to exactly fit the number of cells present
    void         shrinkToFit();

    // sets equal to a 'rows' x 'cols' zero matrix
    void init(int rows, int cols);

    // Sets all entries of the matrix to zero.
    void zero();

    // Determines if we are large enough to justify mulithreading.
    // The 5000 was determined from experiments with UT_Vector.
    bool         shouldMultiThread() const
                 {
#ifdef CELLBE
                     return false;
#else
                     return getNumRows() > 5000;
#endif
                 }

    // does M(row,col) += value
    bool addToElement(int row, int col, T value);

    // Returns the index of the cell for which cell->row >= row is true
    // first.  This could be myCount if row is too high.
    // myCells[result].myRow may be strictly greater than row if it isn't
    // present.
    int findCellFromRow(int row) const;

    // computes M * v
    THREADED_METHOD2_CONST(UT_SparseMatrixT, shouldMultiThread(), multVec,
                    const UT_VectorT<T> &, v,
                    UT_VectorT<T> &, result)
    void multVecPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
                        const UT_JobInfo &info) const;

    /// Do result = result - (M * v)
    THREADED_METHOD2_CONST(UT_SparseMatrixT, shouldMultiThread(), 
                           subtractMultVec,
                           const UT_VectorT<T>&, v,
                           UT_VectorT<T>&, result)
    void subtractMultVecPartial(const UT_VectorT<T>& v, UT_VectorT<T>& result,
                                const UT_JobInfo& info) const;

    // computes M^T * v
    void transposeMultVec(const UT_VectorT<T> &v, UT_VectorT<T> &result) const;

    /// Square of L2-norm of all columns
    /// Use to form the diagonal part of A^tA of an overdetermined system A.
    /// Needed for both Jacobi and Gauss-Seidel iterations.
    void allColNorm2(UT_VectorT<T> &result) const;

    /// Initializes out to a matrix where the (i,j) elements are
    /// (rowstart+i, colstart+j) elements of this.  The rowend and
    /// colend represent the exclusive end - [rowstart,rowend) will
    /// be extracted for numrows = rowend - rowstart
    /// This function will compile *this, if it wasn't before the call.
    /// The resulting matrix will be compiled (by construction).
    void extractSubMatrix(UT_SparseMatrixT<T, IsPaged> &out,
                            int rowstart, int rowend,
                            int colstart, int colend) const;
    /// Extract submatrix without compiling *this.
    /// The resulting submatrix won't be compiled either.
    void extractSubMatrixUncompiled(UT_SparseMatrixT<T, IsPaged> &out,
                            int rowstart, int rowend,
                            int colstart, int colend) const;

    /// Extract the diagonal vector from a matrix. This does not
    /// require the matrix to be compiled, nor does it try to compile it.
    /// For entries (i, j), this method will extract the entries where
    /// (i-j == idx).
    void extractDiagonal(UT_VectorT<T> &out, int idx = 0) const;

    /// Extract everything but the diagonal:
    /// out(i, j) = 0,           if i == j,
    //              *this(i, j), otherwise
    void extractNondiagonal(UT_SparseMatrixT<T, IsPaged> &out) const;

    /// Incomplete Cholesky Factorization.
    /// Does a Cholesky Factorization but does not make any elements non-zero
    /// The result of this operation is an UPPER triangular matrix G
    /// such that A = Gt G.
    /// The input must be a symmetric matrix.
    /// Note that this factorization is not always stable.
    /// Returns 0 on success, -1 if semi-definite (diagonal contained zeros
    /// within the given tolerance) and -2 if not positive definite (diagonal
    /// contained negative numbers)
    int         incompleteCholeskyFactorization(T tol=1e-5);

    /// Modified Incomplete Cholesky
    /// Same as incomplte cholesky, except attempt to account for the
    /// discarded entries by adjusting the diagonal.
    /// tau is a tuning constant.
    int         modifiedIncompleteCholesky(T tau = 0.97, 
                                            T mindiagratio = 0.25,
                                            T tol=1e-5);

    /// Assumes this is a lower triangular matrix.  Solves the equation
    /// A x = b
    /// If the diagonal of A is zero within tolerance, the corresponding
    /// x coordinate is zero.
    /// Returned is the number of artifical zeros places into x.  0 means
    /// the solution encountered no singularities, 10 would mean 10 
    /// singularities.
    int         solveLowerTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
                                    T tol=1e-5) const;
    int         solveUpperTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
                                    T tol=1e-5) const;
    /// Given an upper triangular matrix, solves the lower triangular
    /// transposed of it and negates the result.
    int         solveLowerTriangularTransposeNegate(UT_VectorT<T> &x, 
                                    const UT_VectorT<T> &b,
                                    T tol=1e-5) const;

    // Use the conjugate gradient method to solve Ax=b.
    // NB: 'x' should be initialized with the initial estimate
    bool        solveConjugateGradient(UT_VectorT<T> &x, const UT_VectorT<T> &b,
                        bool (*callback_func)(void *) = 0,
                        void *callback_data = 0, T tol=1e-5,
                        int max_iters = -1) const;

    /// Transposes this matrix.
    THREADED_METHOD(UT_SparseMatrixT, shouldMultiThread(), transpose)
    void        transposePartial(const UT_JobInfo &info);

    /// Makes this a transposed copy of source.  By not working
    /// in place we can avoid a sort if the source is compiled.
    /// If source is not compiled, we compile it.
    void        transposeCompiled(const UT_SparseMatrixT<T, IsPaged> &src);

    /// *this = -*this.
    THREADED_METHOD(UT_SparseMatrixT, shouldMultiThread(), negate)
    void        negatePartial(const UT_JobInfo &info);

    UT_SparseMatrixT<T, IsPaged> &operator=(const UT_SparseMatrixT<T, IsPaged> &m);
    UT_SparseMatrixT<T, IsPaged> &operator*=(T scalar);
    UT_SparseMatrixT<T, IsPaged> &operator+=(const UT_SparseMatrixT<T, IsPaged> &m);

    // Writes out the full matrix as an Row X Column array.  Likely
    // not too useful with big matrices...
    void                printFull(ostream &os) const;

    // Writes out all of the non-zero elements.
    void                printSparse(ostream &os) const;

    // Binary save and load of a sparse matrix
    void                save(ostream &os) const;
    void                load(UT_IStream &is);
    
    /// Reorders the cell array to be sorted by row than column.
    /// Consolidates any duplicate entries by adding them together.
    /// Culls any zero entries.
    /// compile will be a no-op if the matrix's structure has not changed
    /// since the last invocation.
    /// While compilation is expensive, it also performs the collapse/cull
    /// operation, so you may want to invoke it explicitly if you are
    /// repeatedly gathering similar matrices, say with +=.
    void                         compile() const;
    bool                         isCompiled() const { return myCompiledFlag; }
    bool                         isStillSorted() const { return myStillSortedFlag; }

private:
    /// Always recompiles the matrix regardless of the compile flag setting.
    void                         forceCompile() const;

    /// Extract a cell by its index.
    inline ut_MatrixCell        &getCell(int idx) const 
    { 
        if (IsPaged)
            return myCellPages(idx >> CELL_PAGEBITS)[idx & CELL_PAGEMASK]; 
        else
            return myCells[idx]; 
    }


    /// Compile into 4-value cell data structure
    void                         compile4() const;

    int                          myNumRows;
    int                          myNumCols;
    mutable UT_PtrArray<ut_MatrixCell *>        myCellPages;
    mutable ut_MatrixCell       *myCells;
    mutable int                  myCount;
    mutable int                 *my4RowOffsets;
    mutable ut_4MatrixCell      *my4Cells;
    mutable int                  my4Count;
    int                          myMaxSize;
    mutable bool                 myCompiledFlag;
    mutable bool                 myStillSortedFlag;
    mutable bool                 my4CompiledFlag;

    friend class UT_SparseMatrixRowT<T>;
};

typedef UT_SparseMatrixT<fpreal32, false>       UT_SparseMatrixF;
typedef UT_SparseMatrixT<fpreal64, false>       UT_SparseMatrixD;
typedef UT_SparseMatrixT<fpreal64, false>       UT_SparseMatrix;


///
/// This is a highly specialized varient of the SparseMatrix
/// which does not let you change the number or position of cells
/// It has, however, compiled fixed row offsets to make it fast
/// to look up specific rows in the matrix.
///
template <typename T>
class UT_API UT_SparseMatrixRowT
{
    class ut_MatrixCell
    {
    public:
        int myCol;
        T myValue;
    };

public:
    UT_SparseMatrixRowT();
    ~UT_SparseMatrixRowT();

    // Can't mix paged and non paged, and this is a non-paged algorithm
    // Destroys the source matrix in the process as we reuse its
    // myCells array.
    void         buildFrom(UT_SparseMatrixT<T, false> &m,
                            bool invertdiag = false,
                            T tol=1e-5f);

    int getNumRows() const { return myNumRows; }
    int getNumCols() const { return myNumCols; }

    /// Return the amount of memory used by this array.
    int64        getMemoryUsage() const;

    // The 5000 was determined from experiments with UT_Vector.
    bool         shouldMultiThread() const
                 {
                     return getNumRows() > 5000;
                 }

    // Returns the index of the cell for which cell->row >= row is true
    // first.  This could be myCount if row is too high.
    // myCells[result].myRow may be strictly greater than row if it isn't
    // present.
    int         findCellFromRow(int row) const
                { return myRowOffsets(row); }

    // computes M * v
    THREADED_METHOD2_CONST(UT_SparseMatrixRowT, shouldMultiThread(), multVec,
                    const UT_VectorT<T> &, v,
                    UT_VectorT<T> &, result)
    void multVecPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
                        const UT_JobInfo &info) const;

    // computes M * v and accumulates the dot product of v & result
    THREADED_METHOD3_CONST(UT_SparseMatrixRowT, shouldMultiThread(), 
                    multVecAndDot,
                    const UT_VectorT<T> &, v,
                    UT_VectorT<T> &, result,
                    fpreal64 *, dotpq)
    void multVecAndDotPartial(const UT_VectorT<T> &v, UT_VectorT<T> &result,
                        fpreal64 *dotpq,
                        const UT_JobInfo &info) const;

    /// Assumes this is a lower triangular matrix.  Solves the equation
    /// A x = b
    /// If the diagonal of A is zero within tolerance, the corresponding
    /// x coordinate is zero.
    /// Returned is the number of artifical zeros places into x.  0 means
    /// the solution encountered no singularities, 10 would mean 10 
    /// singularities.
    int         solveUpperTriangular(UT_VectorT<T> &x, const UT_VectorT<T> &b,
                                    T tol=1e-5) const;
    /// Given an upper triangular matrix, solves the lower triangular
    /// transposed of it and negates the result.
    int         solveLowerTriangularTransposeNegate(UT_VectorT<T> &x, 
                                    const UT_VectorT<T> &b,
                                    T tol=1e-5) const;

    /// Solves conjugate gradient using our specialized functions.
    /// These allow us to perform some normal and dot operations
    /// while the cache is still hot.
    /// This matrix is the matrix to solve.  The provided GT matrix
    /// is the upper triangular result of cholesky factoriziation.
    /// Norm is hardcoded to 2.
    float       solveConjugateGradient(UT_VectorT<T> &x, const UT_VectorT<T> &b,
                        UT_SparseMatrixRowT<T> *GT, T tol=1e-5,
                        int max_iters = -1) const;
private:
    /// Extract a cell by its index.
    inline ut_MatrixCell        &getCell(int idx) const 
    { 
        return myCells[idx]; 
    }

    int                          myNumRows;
    int                          myNumCols;
    ut_MatrixCell               *myCells;
    int                          myCount;
    UT_IntArray                  myRowOffsets;
    UT_VectorT<T>                myDiagonal;
};

typedef UT_SparseMatrixRowT<fpreal32>   UT_SparseMatrixRowF;
typedef UT_SparseMatrixRowT<fpreal64>   UT_SparseMatrixRowD;

#endif

---