SchurTree.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_MULTIFRONTAL_SYMBOLIC_SCHURTREE_H
#define MYRAMATH_MULTIFRONTAL_SYMBOLIC_SCHURTREE_H

#include <myramath/MYRAMATH_EXPORT.h>
#include <myramath/utility/Number.h>
#include <myramath/utility/Array1.h>
#include <myramath/utility/Array2.h>

// Options pack.
#include <myramath/multifrontal/Options.h>

// Required for AssemblyTree:Node
#include <myramath/multifrontal/symbolic/AssemblyTree.h>

// Required std::libraries.
#include <map>
#include <vector>
#include <iosfwd>
#include <stdint.h>

namespace myra {

// Forward declarations.
class intCRange;
class PatternRange;
class InputStream;
class OutputStream;
class AssemblyTree;


class MYRAMATH_EXPORT SchurTree
  {
  public:

    typedef std::pair<const int*,const int*> Range;
    typedef multifrontal::Options Options;

    // Inner class detail, storage type within WriteMap
    class WriteMapItem
     {
     public:
       int node;
       int column;
       int block;
       int offset;
     };

    // Inner class detail, data structure that records which Node's write to column j.
    class WriteMap2
      {
      public:

        // Default constructor, makes empty container.
        WriteMap2();

        // Constructs from a SchurTree.
        WriteMap2(const SchurTree& stree);

        // Copy constructor.
        WriteMap2(const WriteMap2& that);

        // Presents a range of Node's that write to column j.
        const SchurTree::WriteMapItem* begin(int j) const;
        const SchurTree::WriteMapItem* end(int j) const;

            // Inspection routine, prints debugging information on out.
            void inspect(std::ostream& out) const;

        // Destructor.
       ~WriteMap2();
      
      private:

        // Internal representation.
        std::vector<WriteMapItem> items;
        std::vector<int> strides;

      };

    // Inner class detail, data structure that records which Node's write to column j.
    class WriteMap
      {
      public:

        // Copy constructor.
        WriteMap(const WriteMap& that);

        // Presents a range of Node's that write to column j.
        const int* begin(int j) const;
        const int* end(int j) const;

            // Inspection routine, prints debugging information on out.
            void inspect(std::ostream& out) const;

        // Destructor.
       ~WriteMap();
      
      private:

        // Private constructor, SchurTree manipulates contents directly.
        WriteMap();

        // Internal representation.
        std::vector<int> nodes;
        std::vector<int> strides;
        
        friend class SchurTree;
      };

    // Inner class detail, single Node of the SchurTree
    class Node
      {
      public:
        
        // Unique identifier for this node.
        int id;

        // Number of pivots (p), nonpivots (q), and columns (j)
        int p;
        int q;
        int j;
        std::vector<int> i_pattern;
        std::vector<int> j_pattern;

        // Number of pivot blocks (P) and nonpivot blocks (Q) and column blocks (J)
        int P;
        int Q;
        int J;
        std::vector<int> i_blocking;
        std::vector<int> j_blocking;

        // Striding vector, cumsum(blocking)
        std::vector<int> i_striding;
        std::vector<int> j_striding;

        // Parent/child links.
        Node* parent;
        std::vector<Node*> children;

        // Default constructor, empty Node
        Node();

        // Constructor, requires matching AssemblyTree::Node
        Node(const AssemblyTree::Node& anode);

        // InputStream constructor.
        Node(InputStream& in);

        // Writes to OutputStream.
        void write(OutputStream& out);

        // Returns id of parent (or -1 if *this is a root)
        int parent_id() const;

        // Returns id's of all children.
        std::vector<int> children_ids() const;

            // Inspection routine, prints debugging information on out.
            void inspect(std::ostream& out) const;

        // ------------------------------------ Pattern interrogation.
        
        // Returns pattern of pivots.
        Range p_pattern() const
          { return Range(i_pattern.data(), i_pattern.data()+p); }

        // Returns pattern of nonpivots.
        Range q_pattern() const
          { return Range(i_pattern.data()+p, i_pattern.data()+p+q); }

        // Returns blocking of pivots.
        Range p_blocking() const
          { return Range(i_blocking.data(), i_blocking.data()+P); }

        // Returns blocking of nonpivots.
        Range q_blocking() const
          { return Range(i_blocking.data()+P, i_blocking.data()+P+Q); }

        // Returns pattern of block ib, for ib in [0,P+Q]
        Range iblock_pattern(int ib) const
          { return Range(i_pattern.data()+i_striding[ib], i_pattern.data()+i_striding[ib+1]); }

        // Returns pattern of block jb, for jb in [0,J]
        Range jblock_pattern(int jb) const
          { return Range(j_pattern.data()+j_striding[jb], j_pattern.data()+j_striding[jb+1]); }

        // Returns true iff two SchurTree::Node's j_patterns have common entries.
        static bool j_intersects(const Node& n1, const Node& n2);

        // Finds the (block,offset) to a given column j.
        std::pair<int,int> j_find(int j) const;     

        // ------------------------------------ Space/time queries.

        // Returns (permanent,temporary) storage requirements.
        std::pair<uint64_t,uint64_t> n_words() const;        

        // Returns nodal contribution n_work_solve()
        uint64_t n_work_solve() const;

        // Returns nodal contribution to n_work_syrk()
        uint64_t n_work_rankk() const;

        // Returns nodal contributions to n_work_gemm()
        static uint64_t n_work_gemm(const Node& x, const Node& y);

        // ------------------------------------ Dependency tracking along i-axis     
        
        // Given a contribution block ib in [P,P+Q], returns the set of (parent,iblock)'s that it pushes contributions into.
        Range iblock_up(int ib) const;

        // Given a block ib, returns the set of children[c] iblocks it pulls contributions from.
        Range iblock_down(int ib, int c) const;
        
        // ------------------------------------ Dependency tracking along j-axis

        // Given a block jb in [0,J], returns the set of (parent,jblock)'s that it pushes contributions into.
        Range jblock_up(int jb) const;

        // Given a block jb, returns the set of children[c] jblocks it pulls contributions from.
        Range jblock_down(int jb, int c) const;

      private:
      
        // Caches the child to parent relationships, exposed by iblock_up()
        std::vector<int> iup_data;
        Array1<int> iup_begin;
        Array1<int> iup_end;
        
        // Caches the parent to child relationships, exposed by block_down()
        std::vector<int> idown_data;
        Array2<int> idown_begin;
        Array2<int> idown_end;
        
        // Given a contribution block b in [P,P+Q], returns the set of (parent,block)'s that it pushes contributions into.
        // Implementation detail of finalize()
        std::vector<int> iblock_up_finalize(int b, int blocksize) const;

        // Given a block b, returns the set of children[c] blocks it pulls contributions from.
        // Implementation detail of finalize()
        std::vector<int> iblock_down_finalize(int b, int c, int blocksize) const;

        // Given a sorted range of indices, finds the match for each one in this->pattern. (or -1, if !found)
        // Implementation detail of block_up_finalize() and block_down_finalize()
        std::vector<int> imap(Range range) const;

        // Populates the caches above.
        void ifinalize(int blocksize);

        // Caches the child to parent relationships, exposed by iblock_up()
        std::vector<int> jup_data;
        Array1<int> jup_begin;
        Array1<int> jup_end;
        
        // Caches the parent to child relationships, exposed by block_down()
        std::vector<int> jdown_data;
        Array2<int> jdown_begin;
        Array2<int> jdown_end;

        // Given a block jb in [0,J], returns the set of (parent,jblock)'s that it pushes contributions into.
        // Implementation detail of finalize()
        std::vector<int> jblock_up_finalize(int jb, int blocksize) const; 

        // Given a block jb, returns the set of children[c] jblocks it pulls contributions from.
        // Implementation detail of finalize()
        std::vector<int> jblock_down_finalize(int jb, int c, int blocksize) const;

        // Given a sorted range of j-indices, finds the match for each on in this->j_pattern (or -1, if !found)
        // Implementation detail of finalize()
        std::vector<int> jmap(Range range) const;

        // Populates the caches above.
        void jfinalize(int blocksize);

        friend class SchurTree;
      };

    SchurTree();

    SchurTree(const AssemblyTree& A, const PatternRange& B, Options options = Options::create());

    SchurTree(const AssemblyTree& A, const intCRange& B, int J, Options options = Options::create());

    SchurTree(const SchurTree& that);

    SchurTree(InputStream& in);

    void write(OutputStream& out) const;

    void swap(SchurTree& that);

#ifdef MYRAMATH_ENABLE_CPP11
    SchurTree(SchurTree&& that);
#endif
    
    SchurTree& operator = (SchurTree that);

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

    int n_nodes() const;

    const Node& node(int n) const;
    
    std::vector<int> roots() const;

    std::vector<int> leaves() const;

    std::vector<int> postorder() const;
    std::vector<int> preorder() const;

    const std::vector<int>& s2a() const;

    std::map<int,int> a2s() const;

    WriteMap writemap() const;
    WriteMap2 writemap2() const;
    
    // High water mark should be approximately the sum of .first + .second
    std::pair<uint64_t,uint64_t> n_words() const;

    uint64_t n_work_solve() const;

    uint64_t n_work_rankk() const;

    static uint64_t n_work_gemm(const SchurTree& X, const SchurTree& Y);

    static std::vector<int> intersect(const SchurTree& X, const SchurTree& Y);

    Array2<int> rankk_j_intersects() const;
    Array2<int> rankk_j_intersects_old() const; // TODO remove

    static Array2<int> gemm_j_intersects(const SchurTree& X, const SchurTree& Y);
    static Array2<int> gemm_j_intersects_old(const SchurTree& X, const SchurTree& Y); // TODO remove

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

   ~SchurTree();
  
  private:

    // Implementation details of traversals.
    void postorder_detail(int n, std::vector<int>& answer) const;
    void preorder_detail(int n, std::vector<int>& answer) const;

    // Size of B.
    int I;
    int J;

    // Container of Nodes.
    std::vector<Node*> nodes;

    // Mapping from SchurTree::Node's to AssemblyTree::Node's
    std::vector<int> schur2assembly;

    // Blocksize, induced by underlying AssemblyTree.
    int blocksize;
    
  };

std::ostream& operator<< (std::ostream& out, const SchurTree::WriteMapItem& item);
std::ostream& operator<< (std::ostream& out, const SchurTree::WriteMap& map);
std::ostream& operator<< (std::ostream& out, const SchurTree::WriteMap2& map);
std::ostream& operator<< (std::ostream& out, const SchurTree::Node& stree);
std::ostream& operator<< (std::ostream& out, const SchurTree& stree);

} // namespace myra

#endif