MatrixRange.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_DENSE_MATRIXRANGE_H
#define MYRAMATH_DENSE_MATRIXRANGE_H

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

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

namespace myra {

// Forward declarations.
template<int Arity, class Number> class Expression;
template<class Number> class MatrixRange;
template<class Number> class CMatrixRange;
template<class Number> class VectorRange;
template<class Number> class CVectorRange;
template<class T> class Array1;
template<class T> class Array2;
class OutputStream;
class intCRange;

template<class Number> class LIBPUBLIC MatrixRange
  {
  public:

    Number* begin;        // Pointer to A(0,0)
    int I;                // Number of rows in A
    int J;                // Number of columns in A
    int LD;               // Leading dimension of A

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

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

    MatrixRange();

    MatrixRange(Number* in_begin, int in_I, int in_J, int in_LD);

    MatrixRange(Number* in_begin, int in_I, int in_J);

    void write(OutputStream& out) const;

    // ------------------------------------ Size queries.

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

    uint64_t n_words() const;
    
    // ------------------------------------ General slicing.

    CMatrixRange<Number> add_const() const;

    operator CMatrixRange<Number> () const;

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

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

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

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

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

    MatrixRange<Number> row(int i) const;

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

    MatrixRange<Number> top(int i) const;

    MatrixRange<Number> cut_top(int i) const;

    Pair split_top(int i) const;

    MatrixRange<Number> bottom(int i) const;

    MatrixRange<Number> cut_bottom(int i) const;

    Pair split_bottom(int i) const;

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

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

    MatrixRange<Number> column(int j) const;

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

    MatrixRange<Number> left(int j) const;

    MatrixRange<Number> cut_left(int j) const;

    Pair split_left(int j) const;

    MatrixRange<Number> right(int j) const;

    MatrixRange<Number> cut_right(int j) const;

    Pair split_right(int j) const;

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

    VectorRange<Number> vector(int j) const;

    // ------------------------------------ Mutating underlying contents.

    void swap(const MatrixRange<Number>& that) const;

    // 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
      {
      for (int j = 0; j < J; ++j)
        for (int i = 0; i < I; ++i)
          *pointer(i,j) = f( *pointer(i,j) );
      }
    
    // 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
      {
      if (uplo == 'U')
        {
        for (int j = 0; j < J; ++j)
          for (int i = 0; i <= j && i < I; ++i)
            *pointer(i,j) = f( *pointer(i,j) );
        }
      else if (uplo == 'L')
        {
        for (int j = 0; j < J; ++j)
          for (int i = j; i < I; ++i)
            *pointer(i,j) = f( *pointer(i,j) );
        }
      }

    // 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 N = I < J ? I : J;
      for (int n = 0; n < N; ++n)
        *pointer(n,n) = f( *pointer(n,n) );
      }

    void assign(const CMatrixRange<Number>& that) const;
    void assign(const CMatrixRange<Number>& that, char op) const;
    void assign(const Expression<2,Number>& that) const;

    void operator += (const CMatrixRange<Number>& that) const;
    void operator += (const Expression<2,Number>& that) const;

    void operator -= (const CMatrixRange<Number>& that) const;
    void operator -= (const Expression<2,Number>& that) const;
    
    void operator *= (Number alpha) const;

    void operator /= (Number alpha) const;

    void zero() const;

  private:

    // Returns pointer to (i,j)'th entry.
    Number* pointer(int i, int j) const;

    // Throws if IJ is not the same size as *this.
    void check_size(const std::pair<int,int>& IJ) const;
  
  };

template<class Number> class LIBPUBLIC CMatrixRange
  {
  public:

    const Number* begin;  // Pointer to A(0,0)
    int I;                // Number of rows in A
    int J;                // Number of columns in A
    int LD;               // Leading dimension of A

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

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

    CMatrixRange();

    CMatrixRange(const Number* in_begin, int in_I, int in_J, int in_LD);

    CMatrixRange(const Number* in_begin, int in_I, int in_J);

    void write(OutputStream& out) const;

    // ------------------------------------ Size queries.

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

    uint64_t n_words() const;

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

    MatrixRange<Number> remove_const() const;

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

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

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

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

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

    CMatrixRange<Number> row(int i) const;

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

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

    CMatrixRange<Number> top(int i) const;

    CMatrixRange<Number> cut_top(int i) const;

    Pair split_top(int i) const;

    CMatrixRange<Number> bottom(int i) const;

    CMatrixRange<Number> cut_bottom(int i) const;

    Pair split_bottom(int i) const;

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

    CMatrixRange<Number> column(int j) const;

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

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

    CMatrixRange<Number> left(int j) const;

    CMatrixRange<Number> cut_left(int j) const;

    Pair split_left(int j) const;

    CMatrixRange<Number> right(int j) const;

    CMatrixRange<Number> cut_right(int j) const;

    Pair split_right(int j) const;

    CVectorRange<Number> vector(int j) const;

  private:
    
    // Returns pointer to (i,j)'th entry.
    const Number* pointer(int i, int j) const;
  };

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

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

LIBPUBLIC Array1<CMatrixRange<NumberS> > add_const(const Array1<MatrixRange<NumberS> >& ranges);
LIBPUBLIC Array1<CMatrixRange<NumberD> > add_const(const Array1<MatrixRange<NumberD> >& ranges);
LIBPUBLIC Array1<CMatrixRange<NumberC> > add_const(const Array1<MatrixRange<NumberC> >& ranges);
LIBPUBLIC Array1<CMatrixRange<NumberZ> > add_const(const Array1<MatrixRange<NumberZ> >& ranges);

LIBPUBLIC Array1<MatrixRange<NumberS> > remove_const(const Array1<CMatrixRange<NumberS> >& ranges);
LIBPUBLIC Array1<MatrixRange<NumberD> > remove_const(const Array1<CMatrixRange<NumberD> >& ranges);
LIBPUBLIC Array1<MatrixRange<NumberC> > remove_const(const Array1<CMatrixRange<NumberC> >& ranges);
LIBPUBLIC Array1<MatrixRange<NumberZ> > remove_const(const Array1<CMatrixRange<NumberZ> >& ranges);

LIBPUBLIC Array2<CMatrixRange<NumberS> > add_const(const Array2<MatrixRange<NumberS> >& ranges);
LIBPUBLIC Array2<CMatrixRange<NumberD> > add_const(const Array2<MatrixRange<NumberD> >& ranges);
LIBPUBLIC Array2<CMatrixRange<NumberC> > add_const(const Array2<MatrixRange<NumberC> >& ranges);
LIBPUBLIC Array2<CMatrixRange<NumberZ> > add_const(const Array2<MatrixRange<NumberZ> >& ranges);

LIBPUBLIC Array2<MatrixRange<NumberS> > remove_const(const Array2<CMatrixRange<NumberS> >& ranges);
LIBPUBLIC Array2<MatrixRange<NumberD> > remove_const(const Array2<CMatrixRange<NumberD> >& ranges);
LIBPUBLIC Array2<MatrixRange<NumberC> > remove_const(const Array2<CMatrixRange<NumberC> >& ranges);
LIBPUBLIC Array2<MatrixRange<NumberZ> > remove_const(const Array2<CMatrixRange<NumberZ> >& ranges);

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

} // namespace myra

#endif