SparseMatrix.h

Go to the documentation of this file.

// ========================================================================= //
// This file is part of MyraMath, copyright (c) 2014-2019 by Ryan A Chilton  //
// and distributed by MyraCore, LLC. See LICENSE.txt for license terms.      //
// ========================================================================= //

#ifndef MYRAMATH_SPARSE_SPARSEMATRIX_H
#define MYRAMATH_SPARSE_SPARSEMATRIX_H

#include <myramath/utility/detail/LIBPUBLIC.h>
#include <myramath/utility/Number.h>

#include <myramath/sparse/SparseMatrixRange.h>

#include <utility>
#include <vector>
#include <iosfwd>

namespace myra {

// Forward declarations.
class InputStream;
class OutputStream;
template<int Arity, class Number> class Expression;
template<class Number> class Matrix;
template<class Number> class MatrixRange;
template<class Number> class SparseMatrix;
template<class Number> class SparseMatrixBuilder;
class Pattern;
class PatternRange;
template<class T> class Array1;
template<class T> class Array2;

namespace detail {

// Some of these functions can be implemented more efficiently with access to internal 
// data structures, they're forward declared here so they can be granted friendship later.
class permute_wrapper;
class horzcat_wrapper;
class vertcat_wrapper;
class bothcat_wrapper;
class diagcat_wrapper;
class transpose_wrapper;
class multiply_wrapper;
class flip_wrapper;

} // namespace detail

template<class Number> class LIBPUBLIC SparseMatrix
  {
  public:

    typedef std::pair< CSparseMatrixRange<Number>, CSparseMatrixRange<Number> > CPair;
    typedef std::pair<  SparseMatrixRange<Number>,  SparseMatrixRange<Number> >  Pair;

    // ------------------------------------ Construction, serialization, value semantics.

    SparseMatrix();

    explicit SparseMatrix(const PatternRange& P);

    SparseMatrix(const SparseMatrix& that);

    void swap(SparseMatrix& that);
    void swap(Pattern& that);

#ifdef MYRAMATH_ENABLE_CPP11
    SparseMatrix(SparseMatrix&& that);
    SparseMatrix(Pattern&& that);
#endif

    SparseMatrix& operator = (SparseMatrix that);

    explicit SparseMatrix(const CSparseMatrixRange<Number>& that);

    explicit SparseMatrix(InputStream& in);

    void write(OutputStream& out) const;

    const SparseMatrix<Number>& add_const() const;

   ~SparseMatrix();

    // ------------------------------------ Size and pattern queries.

    std::pair<int,int> size() const;

    int n_nonzeros() const;

    int n_nonzeros(int j) const;

    bool test(int i, int j) const;

    PatternRange pattern() const; 

    std::vector<int> pattern(int j) const;

    // ------------------------------------ General slicing.

    CSparseMatrixRange<Number> range() const;
     SparseMatrixRange<Number> range()      ;

    operator CSparseMatrixRange<Number> () const;
    operator  SparseMatrixRange<Number> ()      ;

    CSparseMatrixRange<Number> window(int i0, int i1, int j0, int j1) const;
     SparseMatrixRange<Number> window(int i0, int i1, int j0, int j1)      ;

    Array2<CSparseMatrixRange<Number> > windows(const intCRange& i, const intCRange& j) const;
    Array2< SparseMatrixRange<Number> > windows(const intCRange& i, const intCRange& j)      ;

    const Number* pointer(int i, int j) const;
          Number* pointer(int i, int j)      ;

    const Number& operator() (int i, int j) const;
          Number& operator() (int i, int j)      ;

    const Number& at(int i, int j) const;
          Number& at(int i, int j)      ;

    // ------------------------------------ Iterating.

    CSparseMatrixIterator<Number> begin() const;
    CSparseMatrixIterator<Number> end  () const;
     SparseMatrixIterator<Number> begin()      ;
     SparseMatrixIterator<Number> end  ()      ;
    
    CSparseMatrixIterator<Number> begin(int j) const;
    CSparseMatrixIterator<Number> end  (int j) const;
     SparseMatrixIterator<Number> begin(int j)      ;
     SparseMatrixIterator<Number> end  (int j)      ;

    // ------------------------------------ Row slicing.

    CSparseMatrixRange<Number> row(int i) const;
     SparseMatrixRange<Number> row(int i)      ;

    CSparseMatrixRange<Number> rows(int i0, int i1) const;
     SparseMatrixRange<Number> rows(int i0, int i1)      ;

    CSparseMatrixRange<Number> top(int i) const;
     SparseMatrixRange<Number> top(int i)      ;

    CSparseMatrixRange<Number> cut_top(int i) const;
     SparseMatrixRange<Number> cut_top(int i)      ;

    CPair split_top(int i) const;
     Pair split_top(int i)      ;

    CSparseMatrixRange<Number> bottom(int i) const;
     SparseMatrixRange<Number> bottom(int i)      ;

    CSparseMatrixRange<Number> cut_bottom(int i) const;
     SparseMatrixRange<Number> cut_bottom(int i)      ;

    CPair split_bottom(int i) const;
     Pair split_bottom(int i)      ;

    Array1<CSparseMatrixRange<Number> > rows(const intCRange& i) const;
    Array1< SparseMatrixRange<Number> > rows(const intCRange& i)      ;

    // ------------------------------------ Column slicing.

    CSparseMatrixRange<Number> column(int j) const;
     SparseMatrixRange<Number> column(int i)      ;

    CSparseMatrixRange<Number> columns(int j0, int j1) const;
     SparseMatrixRange<Number> columns(int j0, int j1)      ;

    CSparseMatrixRange<Number> left(int j) const;
     SparseMatrixRange<Number> left(int j)      ;

    CSparseMatrixRange<Number> cut_left(int j) const;
     SparseMatrixRange<Number> cut_left(int j)      ;

    CPair split_left(int j) const;
     Pair split_left(int j)      ;

    CSparseMatrixRange<Number> right(int j) const;
     SparseMatrixRange<Number> right(int j)      ;    

    CSparseMatrixRange<Number> cut_right(int j) const;
     SparseMatrixRange<Number> cut_right(int j)      ;    

    CPair split_right(int j) const;
     Pair split_right(int j)      ;

    Array1<CSparseMatrixRange<Number> > columns(const intCRange& j) const;
    Array1< SparseMatrixRange<Number> > columns(const intCRange& j)      ;

    // ------------------------------------ Manipulating numeric contents.

    // Functor f should be a have an operator(), that takes a Number and returns a Number.
    template<class Functor> void transform(const Functor& f)
      { this->range().transform(f); }

    // Functor f should be a have an operator(), that takes a Number and returns a Number.
    template<class Functor> void transform_triangle(const Functor& f, char uplo)
      { this->range().transform_triangle(f,uplo); }

    // Functor f should be a have an operator(), that takes a Number and returns a Number.
    template<class Functor> void transform_diagonal(const Functor& f)
      { this->range().transform_diagonal(f); }

    SparseMatrix& operator = (const CSparseMatrixRange<Number>& that);

    void operator *= (Number alpha);

    void operator /= (Number alpha);
    
    SparseMatrix<Number> operator -();

    // ------------------------------------ Static generators.

    static SparseMatrix<Number> identity(int IJ);

    static SparseMatrix<Number> random(int I, int J, int N);

    static SparseMatrix<Number> random(const PatternRange& P);

    static SparseMatrix<Number> zeros(int I, int J);

    static SparseMatrix<Number> zeros(const PatternRange& P);

    static SparseMatrix<Number> ones(const PatternRange& P);

    static SparseMatrix<Number> fill(const PatternRange& P, Number c);

    static SparseMatrix<Number> fill(const intCRange& i, const intCRange& j, const std::vector<Number>& v);

    static SparseMatrix<Number> swap(int I, int J, std::vector<int>& Ao, std::vector<int>& Ai, std::vector<Number>& Av);

    static SparseMatrix<Number> evaluate(const PatternRange& P, const Expression<2,Number>& e);

    Matrix<Number> make_Matrix() const;
    void make_Matrix(const MatrixRange<Number>& A) const;

    void debug(std::ostream& out) const;

  private:

    // Returns offset to *this A(i,j). Returns -1 if (i,j) is a structural zero. Throws if (i,j) are out of (I,J) bounds.
    int offset(int i, int j) const;

    // Stores size.
    int I;
    int J;

    // Stores offsets for each column of A.
    std::vector<int> Ao;

    // Holds row indices (i). The sorted i-pattern of column A(:,j) is given by Ai[Ao[j]:Ao[j+1]]
    std::vector<int> Ai;

    // Holds numeric values (v). The numeric values of column A(:,j) are given by Av[Ao[j]:Ao[j+1]]
    std::vector<Number> Av;

    // Nonmember friend functions that are granted access to internal data structures.
    friend class SparseMatrixBuilder<Number>;
    friend class detail::permute_wrapper;
    friend class detail::horzcat_wrapper;
    friend class detail::vertcat_wrapper;
    friend class detail::bothcat_wrapper;
    friend class detail::diagcat_wrapper;
    friend class detail::transpose_wrapper;
    friend class detail::multiply_wrapper;
    friend class detail::flip_wrapper;
  };

template<class Number> class ReflectNumber< SparseMatrix<Number> >
  { public: typedef Number type; };

LIBPUBLIC SparseMatrix<NumberS> operator + (const CSparseMatrixRange<NumberS>& A, const CSparseMatrixRange<NumberS>& B);
LIBPUBLIC SparseMatrix<NumberD> operator + (const CSparseMatrixRange<NumberD>& A, const CSparseMatrixRange<NumberD>& B);
LIBPUBLIC SparseMatrix<NumberC> operator + (const CSparseMatrixRange<NumberC>& A, const CSparseMatrixRange<NumberC>& B);
LIBPUBLIC SparseMatrix<NumberZ> operator + (const CSparseMatrixRange<NumberZ>& A, const CSparseMatrixRange<NumberZ>& B);

LIBPUBLIC SparseMatrix<NumberS> operator - (const CSparseMatrixRange<NumberS>& A, const CSparseMatrixRange<NumberS>& B);
LIBPUBLIC SparseMatrix<NumberD> operator - (const CSparseMatrixRange<NumberD>& A, const CSparseMatrixRange<NumberD>& B);
LIBPUBLIC SparseMatrix<NumberC> operator - (const CSparseMatrixRange<NumberC>& A, const CSparseMatrixRange<NumberC>& B);
LIBPUBLIC SparseMatrix<NumberZ> operator - (const CSparseMatrixRange<NumberZ>& A, const CSparseMatrixRange<NumberZ>& B);

LIBPUBLIC SparseMatrix<NumberS> operator * (NumberS alpha, const CSparseMatrixRange<NumberS>& A);
LIBPUBLIC SparseMatrix<NumberD> operator * (NumberD alpha, const CSparseMatrixRange<NumberD>& A);
LIBPUBLIC SparseMatrix<NumberC> operator * (NumberC alpha, const CSparseMatrixRange<NumberC>& A);
LIBPUBLIC SparseMatrix<NumberZ> operator * (NumberZ alpha, const CSparseMatrixRange<NumberZ>& A);

LIBPUBLIC SparseMatrix<NumberS> operator * (const CSparseMatrixRange<NumberS>& A, NumberS alpha);
LIBPUBLIC SparseMatrix<NumberD> operator * (const CSparseMatrixRange<NumberD>& A, NumberD alpha);
LIBPUBLIC SparseMatrix<NumberC> operator * (const CSparseMatrixRange<NumberC>& A, NumberC alpha);
LIBPUBLIC SparseMatrix<NumberZ> operator * (const CSparseMatrixRange<NumberZ>& A, NumberZ alpha);

LIBPUBLIC SparseMatrix<NumberC> make_complex(const CSparseMatrixRange<NumberS>& A);
LIBPUBLIC SparseMatrix<NumberZ> make_complex(const CSparseMatrixRange<NumberD>& A);

LIBPUBLIC SparseMatrix<NumberC> make_complex(const CSparseMatrixRange<NumberS>& R, const CSparseMatrixRange<NumberS>& I);
LIBPUBLIC SparseMatrix<NumberZ> make_complex(const CSparseMatrixRange<NumberD>& R, const CSparseMatrixRange<NumberD>& I);

LIBPUBLIC SparseMatrix<NumberS> realpart(const CSparseMatrixRange<NumberC>& A);
LIBPUBLIC SparseMatrix<NumberD> realpart(const CSparseMatrixRange<NumberZ>& A);

LIBPUBLIC SparseMatrix<NumberS> imagpart(const CSparseMatrixRange<NumberC>& A);
LIBPUBLIC SparseMatrix<NumberD> imagpart(const CSparseMatrixRange<NumberZ>& A);

} // namespace myra

#endif