LowerMatrixRange.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_LOWERMATRIXRANGE_H
#define MYRAMATH_DENSE_LOWERMATRIXRANGE_H

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

#include <iosfwd>
#include <utility>

namespace myra {

// Forward declarations.
template<int Arity, class Number> class Expression;
template<class Number> class Matrix;
template<class Number> class MatrixRange;
template<class Number> class CMatrixRange;
template<class Number> class LowerMatrixRange;
template<class Number> class CLowerMatrixRange;
class OutputStream;

template<class Number> class LIBPUBLIC LowerMatrixRange
  {
  public:

    Number* begin;
    int N1;
    int N2;
    int LD;

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

    LowerMatrixRange();

    LowerMatrixRange(Number* in_begin, int N);

    void write(OutputStream& out) const;

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

    int size() const;

    uint64_t n_words() const;

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

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

    CLowerMatrixRange<Number> add_const() const;

    operator CLowerMatrixRange<Number> () const;

    MatrixRange<Number> rectangle() const;

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

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

    // Optionally, the 'T'ranspose or 'H'ermitian of *this can be copied into the upper triangle.
    Matrix<Number> make_Matrix(char op = 'N') const;
    void make_Matrix(const MatrixRange<Number>& A, char op = 'N') const;

    // ------------------------------------ Lower triangular slicing.

    MatrixRange<Number> range11 () const;
    MatrixRange<Number> range21 () const;
    MatrixRange<Number> range22t() 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
      {
      int N = N1+N2;
      for (int j = 0; j < N; ++j)
        for (int i = j; i < N; ++i)
          *pointer(i,j) = f( *pointer(i,j) );
      }

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

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

    void operator -= (const CLowerMatrixRange<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 a pointer to L(i,j)
    Number* pointer(int i, int j) const;

    // Throws if that is not the same size as *this.
    void check_size(const CLowerMatrixRange<Number>& that) const;
    
    // Throws if S is not the same size as *this.
    void check_size(int S) const;
  };

template<class Number> class LIBPUBLIC CLowerMatrixRange
  {
  public:

    const Number* begin;
    int N1;
    int N2;
    int LD;

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

    CLowerMatrixRange();

    CLowerMatrixRange(const Number* in_begin, int N);

    void write(OutputStream& out) const;

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

    int size() const;

    uint64_t n_words() const;

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

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

    LowerMatrixRange<Number> remove_const() const;

    CMatrixRange<Number> rectangle() const;

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

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

    // Optionally, the 'T'ranspose or 'H'ermitian of *this can be copied into the upper triangle.
    Matrix<Number> make_Matrix(char op = 'N') const;
    void make_Matrix(const MatrixRange<Number>& A, char op = 'N') const;

    // ------------------------------------ Lower triangular slicing.

    CMatrixRange<Number> range11 () const;
    CMatrixRange<Number> range21 () const;
    CMatrixRange<Number> range22t() const;

  private:

    // Returns a pointer to L(i,j)
    const Number* pointer(int i, int j) const;
    
  };

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

/*

  Stores a lower triangular matrix in lapack-friendly format, by chopping
  off the lower right corner, transposing it, and putting it "on top" of 
  the rest. Pictures are worth a thousand words, here's a few to cover it:
  
  Examples for even N:
  
    N=2
               C   N1 = 1                     2
    A          A   N2 = 1      0         ->   0
    B C   ->   B   LD = 3      1 2            1
  
  
    N=4
                C C                           7  8
    A           A C   N1 = 2   0              0  9 
    A A     ->  A A   N2 = 2   1 4       ->   1  4  
    B B C       B B   LD = 5   2 5 7          2  5
    B B C C     B B            3 6 8 9        3  6
              
    N=6

                    C C C                                 15 16 18
    A               A C C   N1 = 3    0                    0 17 19
    A A             A A C   N2 = 3    1  6                 1  6 20
    A A A      ->   A A A   LD = 7    2  7 11        ->    2  7 11
    B B B C         B B B             3  8 12 15           3  8 12
    B B B C C       B B B             4  9 13 16 18        4  9 13
    B B B C C C     B B B             5 10 14 17 19 20     5 10 14
  
  Examples for odd N:
  
    N=1
  
             N1 = 1
    A -> A   N2 = 0   0  ->  0
             LD = 1
  
    N=3
  
    A          A C   N1 = 2   0           0 5
    A A    ->  A A   N2 = 1   1 3    ->   1 3 
    B B C      B B   LD = 3   2 4 5       2 4
      
    N=5 
  
    A            A C C            1                  0 12 13
    A A          A A C   N1 = 3   2 6                1  5 14 
    A A A    ->  A A A   N2 = 2   3 7 10             2  6  9
    B B B C      B B B   LD = 5   4 8 11 13     ->   3  7 10
    B B B C C    B B B            5 9 12 14 15       4  8 11
  
    N=7 
  
    A              A C C C            0                      0 22 23 25
    A A            A A C C   N1 = 4   1  7                   1  7 24 26
    A A A    ->    A A A C   N2 = 3   2  8 13                2  8 13 27
    A A A A        A A A A   LD = 7   3  9 14 18        ->   3  9 14 18
    B B B B C      B B B B            4 10 15 19 22          4 10 15 19
    B B B B C C    B B B B            5 11 16 20 23 25       5 11 16 20
    B B B B C C C  B B B B            6 12 17 21 24 26 27    6 12 17 21

  MatrixRange's can be obtained over the the various pieces. For example (N=7):

      [ A C C C ]                    [ A C C C ]
      [ A A C C ]                    [ A A C C ] == .range11()
      [ A A A C ]                    [ A A A C ]
      [ A A A A ] == .rectangle()    [ A A A A ]
      [ B B B B ]                      B B B B
      [ B B B B ]                      B B B B
      [ B B B B ]                      B B B B

        A C C C                       A [C C C]
        A A C C                       A [A C C] == .range22t()
        A A A C                       A [A A C]
        A A A A                       A  A A A
      [ B B B B ]                     B  B B B
      [ B B B B ] == .range21()       B  B B B
      [ B B B B ]                     B  B B B

  This storage scheme allows for efficient "column major"/"BLAS3" arithmetic, 
  but maintains the factor of two space advantage that packed formats enjoy 
  compared to just ignoring/wasting the upper triangle of rectangular storage.

  */
  
} // namespace myra

#endif