SparseMatrixRange.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_SPARSEMATRIXRANGE_H
#define MYRAMATH_SPARSE_SPARSEMATRIXRANGE_H

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

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

namespace myra {

// Forward declarations.
template<int Arity, class Number> class Expression;
class intCRange;
template<class Number> class Matrix;
template<class Number> class MatrixRange;
class PatternRange;
template<class Number> class SparseMatrix;
template<class Number> class CSparseMatrixRange;
template<class Number> class  SparseMatrixRange;
template<class Number> class CSparseMatrixIterator;
template<class Number> class  SparseMatrixIterator;
template<class T> class Array1;
template<class T> class Array2;
class OutputStream;

template<class Number> class LIBPUBLIC SparseMatrixRange
  {
  public:

    // Data pointers for underlying SparseMatrix
    const int* Ao; // Offsets for each j into Ai[] and Av[] data.
    const int* Ai; // Row indices (i) in compressed-sparse-column format.
    Number* Av;    // Numeric values (v) in compressed-sparse-column format.
     
    // Extents of this Range within underlying SparseMatrix.
    int I0,I1;     // I-extents of Range
    int J0,J1;     // J-extents of Range
    bool sliced;   // If true, this SparseMatrixRange has been sliced.

    // Useful typedefs.
    typedef std::pair<SparseMatrixRange<Number>, SparseMatrixRange<Number> > Pair;

    // ------------------------------------ Construction and range semantics.

    SparseMatrixRange();

    SparseMatrixRange(const int* in_Ao, const int* in_Ai, Number* in_Av, int I, int J);

    void write(OutputStream& out) const;
    
    // ------------------------------------ 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;

    // ------------------------------------ Iterating.
    typedef SparseMatrixIterator<Number> Iterator;

    SparseMatrixIterator<Number> begin() const;
    SparseMatrixIterator<Number> end() const;

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

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

    CSparseMatrixRange<Number> add_const() const;

    operator CSparseMatrixRange<Number> () const;

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

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

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

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

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

    Matrix<Number> make_Matrix() const;
    void make_Matrix(const MatrixRange<Number>& A) const;
    
    // ------------------------------------ Row slicing.
    
    SparseMatrixRange<Number> row(int i) const;

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

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

    SparseMatrixRange<Number> cut_top(int i) const;

    Pair split_top(int i) const;

    SparseMatrixRange<Number> bottom(int i) const;

    SparseMatrixRange<Number> cut_bottom(int i) const;

    Pair split_bottom(int i) const;

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

    SparseMatrixRange<Number> column(int j) const;

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

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

    SparseMatrixRange<Number> left(int j) const;

    SparseMatrixRange<Number> cut_left(int j) const;

    Pair split_left(int j) const;

    SparseMatrixRange<Number> right(int j) const;

    SparseMatrixRange<Number> cut_right(int j) const;

    Pair split_right(int j) const;

    // ------------------------------------ Mutating underlying 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) const
      { 
      typedef SparseMatrixIterator<Number> Iterator;
      Iterator b = this->begin();
      Iterator e = this->end();
      for (Iterator aij = b; aij != e; ++aij) 
        aij.value() = f(aij.value()); 
      }

    // 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) const
      {
      typedef SparseMatrixIterator<Number> Iterator;
      Iterator b = this->begin();
      Iterator e = this->end();
      if (uplo == 'U')
        {
        for (Iterator aij = b; aij != e; ++aij) 
          if (aij.i() <= aij.j())
            aij.value() = f(aij.value());
        }
      else if (uplo == 'L')
        {
        for (Iterator aij = b; aij != e; ++aij) 
          if (aij.i() >= aij.j())
            aij.value() = f(aij.value());
        }
      }

    // 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) const
      {
      int II = this->I();
      int JJ = this->J();
      int N = II < JJ ? II : JJ;
      for (int n = 0; n < N; ++n)
        {
        Number* p = this->pointer(n,n);
        if (p) 
          *p = f(*p);
        }
      }

    void operator *= (Number alpha) const;
    
    void operator /= (Number alpha) const;

  private:


    // Implementation details, helps make subranges.
    int I() const;
    int J() const;

    // Implementation detail for iterator compare, verifies both iterator's are over same range.
    bool same (const SparseMatrixRange<Number>& that) const;

    // Implementation detail for iterator compare, verifies both iterator's are over same range.
    bool notsame (const SparseMatrixRange<Number>& that) const;

    // A special memory location so that a default-constructed 0x0 SparseMatrixRange can point to valid Ai[] data
    static int Ai_zero;
  
    friend class SparseMatrixIterator<Number>;
    friend class CSparseMatrixRange<Number>;
    friend class SparseMatrix<Number>;

  };

// Initializes SparseMatrixRange<Number>::Ai_zero.
template<class Number> int SparseMatrixRange<Number>::Ai_zero = 0;

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

template<class Number> class LIBPUBLIC SparseMatrixIterator
  {
  public:

    void operator ++ ();
    bool operator == (const SparseMatrixIterator<Number>& that) const;
    bool operator != (const SparseMatrixIterator<Number>& that) const;

    int i() const;
    int j() const;
    Number& value() const;

  private:

    // Implementation detail of increment()
    void seek();

    // Should call this when jay is modified, sets up traversal A(I0:I1,jay)
    void j_initialize(int jj);

    // Implementation detail of initialize(), returns offset to A(ii,jay)
    int i_initialize(int ii);
  
    // Copy of underlying range data.
    SparseMatrixRange<Number> r;

    // State members of iterator.
    int jay; // Column iterator jay
    int b;   // Offset to A(I0,jay), the begin of column jay
    int e;   // Offset to A(I1,jay), the end of column jay
    int o;   // Offset to current nonzero.

    friend class SparseMatrixRange<Number>;
    
  };

template<class Number> class LIBPUBLIC CSparseMatrixRange
  {
  public:

    // Data pointers for underlying SparseMatrix
    const int* Ao;    // Offsets for each j into Ai[] and Av[] data.
    const int* Ai;    // Row indices (i) in compressed-sparse-column format.
    const Number* Av; // Numeric values (v) in compressed-sparse-column format.
     
    // Extents of this Range within underlying SparseMatrix.
    int I0,I1;        // I-extents of Range
    int J0,J1;        // J-extents of Range
    bool sliced;      // If true, this SparseMatrixRange has been sliced.

    // Useful typedefs.
    typedef std::pair<CSparseMatrixRange<Number>, CSparseMatrixRange<Number> > Pair;

    CSparseMatrixRange();

    CSparseMatrixRange(const int* in_Ao, const int* in_Ai, const Number* in_Av, int I, int J);

    void write(OutputStream& out) const;

    // ------------------------------------ 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;

    // ------------------------------------ Iterating.
    typedef CSparseMatrixIterator<Number> Iterator;

    CSparseMatrixIterator<Number> begin() const;
    CSparseMatrixIterator<Number> end() const;

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

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

    SparseMatrixRange<Number> remove_const() const;

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

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

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

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

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

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

    // ------------------------------------ Row slicing.
    
    CSparseMatrixRange<Number> row(int i) const;

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

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

    CSparseMatrixRange<Number> cut_top(int i) const;

    Pair split_top(int i) const;

    CSparseMatrixRange<Number> bottom(int i) const;

    CSparseMatrixRange<Number> cut_bottom(int i) const;

    Pair split_bottom(int i) const;

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

    CSparseMatrixRange<Number> column(int j) const;

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

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

    CSparseMatrixRange<Number> left(int j) const;

    CSparseMatrixRange<Number> cut_left(int j) const;

    Pair split_left(int j) const;

    CSparseMatrixRange<Number> right(int j) const;

    CSparseMatrixRange<Number> cut_right(int j) const;

    Pair split_right(int j) const;

  private:

    // Implementation details, helps make subranges.
    int I() const;
    int J() const;

    // Implementation detail for iterator compare, verifies both iterators are over same range.
    bool operator == (const CSparseMatrixRange<Number>& that) const;

    // Implementation detail for iterator compare, verifies both iterators are over same range.
    bool operator != (const CSparseMatrixRange<Number>& that) const;

    // A special memory location so that a default-constructed 0x0 CSparseMatrixRange can point to valid Ai[] data
    static int Ai_zero;

    friend class CSparseMatrixIterator<Number>;
    friend class SparseMatrixRange<Number>;    
    friend class SparseMatrix<Number>;
    
  };

// Initializes CSparseMatrixRange<Number>::Ai_zero.
template<class Number> int CSparseMatrixRange<Number>::Ai_zero = 0;

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

template<class Number> class LIBPUBLIC CSparseMatrixIterator
  {
  public:

    void operator ++ ();
    bool operator == (const CSparseMatrixIterator<Number>& that) const;
    bool operator != (const CSparseMatrixIterator<Number>& that) const;

    int i() const;
    int j() const;
    const Number& value() const;
    
  private:

    // Implementation detail of increment()
    void seek();

    // Should call this when jay is modified, sets up traversal A(I0:I1,jay)
    void j_initialize(int jj);

    // Implementation detail of initialize(), returns offset to A(ii,jay)
    int i_initialize(int ii);
  
    // Copy of underlying range data.
    CSparseMatrixRange<Number> r;

    // State members of iterator.
    int jay; // Column iterator jay
    int b;   // Offset to A(I0,jay), the begin of column jay
    int e;   // Offset to A(I1,jay), the end of column jay
    int o;   // Offset to current nonzero.

    friend class CSparseMatrixRange<Number>;
    
  };

LIBPUBLIC std::ostream& operator << (std::ostream& out, const CSparseMatrixRange<NumberS>& A);
LIBPUBLIC std::ostream& operator << (std::ostream& out, const CSparseMatrixRange<NumberD>& A);
LIBPUBLIC std::ostream& operator << (std::ostream& out, const CSparseMatrixRange<NumberC>& A);
LIBPUBLIC std::ostream& operator << (std::ostream& out, const CSparseMatrixRange<NumberZ>& A);

} // namespace

#endif