Multiply sparse matrices

include/El/blas_like/level1/DiagonalScale.hpp

@@ -231,8 +231,7 @@ void DiagonalScale
               LogicError("The size of d must match the width of A");
         )
         A.InitializeMultMeta();
-        const auto& meta = A.multMeta;
-
+        const auto& meta = A.LockedDistGraph().multMeta;
         // Pack the send values 
         const Int numSendInds = meta.sendInds.size();
         vector<T> sendVals;

include/El/blas_like/level1/DiagonalSolve.hpp

@@ -257,8 +257,7 @@ void DiagonalSolve
     else
     {
         A.InitializeMultMeta();
-        const auto& meta = A.multMeta;
-
+        const auto& meta = A.LockedDistGraph().multMeta;
         // Pack the send values
         const Int numSendInds = meta.sendInds.size();
         vector<F> sendVals;
@@ -308,8 +307,7 @@ void SymmetricDiagonalSolve
     const Int firstLocalRow = d.FirstLocalRow();
 
     A.InitializeMultMeta();
-    const auto& meta = A.multMeta;
-
+    const auto& meta = A.LockedDistGraph().multMeta;
     // Pack the send values
     const Int numSendInds = meta.sendInds.size();
     vector<Real> sendVals( numSendInds );

include/El/blas_like/level1/EntrywiseMap.hpp

@@ -160,7 +160,6 @@ void EntrywiseMap
     B.remoteVals_.resize( numRemoteEntries );
     for( Int k=0; k<numRemoteEntries; ++k )
         B.remoteVals_[k] = func(A.remoteVals_[k]);
-    B.multMeta = A.multMeta;
     B.ProcessQueues();
 }
 

include/El/blas_like/level3.hpp

@@ -159,6 +159,12 @@ void Multiply
   T alpha, const SparseMatrix<T>& A, const Matrix<T>& X,
   T beta,                                  Matrix<T>& Y );
 
+template<typename T>
+void Multiply
+( Orientation orientation,
+  T alpha, const Graph& A, const Matrix<T>& X,
+  T beta,                        Matrix<T>& Y );
+
 template<typename T>
 void Multiply
 ( Orientation orientation,

include/El/core.hpp

@@ -235,7 +235,7 @@ class DistMatrix;
 #include <El/core/Matrix/decl.hpp>
 #include <El/core/Graph/decl.hpp>
 #include <El/core/DistMap/decl.hpp>
-#include <El/core/DistGraph/decl.hpp>
+#include <El/core/DistGraph.hpp>
 #include <El/core/SparseMatrix/decl.hpp>
 #include <El/core/DistSparseMatrix/decl.hpp>
 #include <El/core/DistMultiVec/decl.hpp>

include/El/core/DistGraph.hpp

@@ -11,6 +11,206 @@
 #ifndef EL_CORE_DISTGRAPH_HPP
 #define EL_CORE_DISTGRAPH_HPP
 
-#include <El/core/DistGraph/decl.hpp>
+#include <set>
+
+namespace El {
+
+struct DistGraphMultMeta
+{
+    bool ready;
+    // NOTE: The 'send' and 'recv' roles reverse for adjoint multiplication
+    Int numRecvInds;
+    vector<int> sendSizes, sendOffs,
+                recvSizes, recvOffs;
+    vector<Int> sendInds, colOffs;
+
+    DistGraphMultMeta() : ready(false), numRecvInds(0) { }
+
+    void Clear()
+    {
+        ready = false;
+        numRecvInds = 0;
+        SwapClear( sendSizes );
+        SwapClear( recvSizes );
+        SwapClear( sendOffs );
+        SwapClear( recvOffs );
+        SwapClear( sendInds );
+        SwapClear( colOffs );
+    }
+
+    const DistGraphMultMeta& operator=( const DistGraphMultMeta& meta )
+    {
+        ready = meta.ready;
+        numRecvInds = meta.numRecvInds;
+        sendSizes = meta.sendSizes;
+        sendOffs = meta.sendOffs;
+        recvSizes = meta.recvSizes;
+        recvOffs = meta.recvOffs;
+        sendInds = meta.sendInds;
+        colOffs = meta.colOffs;
+        return *this;
+    }
+};
+
+
+
+using std::set;
+
+// Forward declare ldl::DistFront
+namespace ldl { template<typename F> struct DistFront; }
+
+// Use a simple 1d distribution where each process owns a fixed number of 
+// sources:
+//     if last process,  numSources - (commSize-1)*floor(numSources/commSize)
+//     otherwise,        floor(numSources/commSize)
+class DistGraph
+{
+public:
+    // Constructors and destructors
+    // ============================
+    DistGraph( mpi::Comm comm=mpi::COMM_WORLD );
+    DistGraph( Int numSources, mpi::Comm comm=mpi::COMM_WORLD );
+    DistGraph( Int numSources, Int numTargets, mpi::Comm comm=mpi::COMM_WORLD );
+    DistGraph( const Graph& graph );
+    // TODO: Move constructor
+    DistGraph( const DistGraph& graph );
+    ~DistGraph();
+
+    // Assignment and reconfiguration
+    // ==============================
+
+    // Making a copy
+    // -------------
+    const DistGraph& operator=( const Graph& graph );
+    const DistGraph& operator=( const DistGraph& graph );
+    // TODO: Move assignment
+
+    // Make a copy of a subgraph
+    // -------------------------
+    DistGraph operator()( Range<Int> I, Range<Int> J ) const;
+    DistGraph operator()( Range<Int> I, const vector<Int>& J ) const;
+    DistGraph operator()( const vector<Int>& I, Range<Int> J ) const;
+    DistGraph operator()( const vector<Int>& I, const vector<Int>& J ) const;
+
+    // Changing the graph size
+    // -----------------------
+    void Empty( bool freeMemory=true );
+    void Resize( Int numVertices );
+    void Resize( Int numSources, Int numTargets );
+
+    // Changing the distribution
+    // -------------------------
+    void SetComm( mpi::Comm comm );
+
+    // Assembly
+    // --------
+    void Reserve( Int numLocalEdges, Int numRemoteEdges=0 );
+
+    // Safe edge insertion/removal procedure
+    void Connect( Int source, Int target );
+    void ConnectLocal( Int localSource, Int target );
+    void Disconnect( Int source, Int target ); 
+    void DisconnectLocal( Int localSource, Int target );
+
+    void FreezeSparsity() EL_NO_EXCEPT;
+    void UnfreezeSparsity() EL_NO_EXCEPT;
+    bool FrozenSparsity() const EL_NO_EXCEPT;
+
+    // For inserting/removing a sequence of edges and then forcing consistency
+    void QueueConnection( Int source, Int target, bool passive=false ) 
+    EL_NO_RELEASE_EXCEPT;
+    void QueueLocalConnection( Int localSource, Int target )
+    EL_NO_RELEASE_EXCEPT; 
+    void QueueDisconnection( Int source, Int target, bool passive=false )
+    EL_NO_RELEASE_EXCEPT;
+    void QueueLocalDisconnection( Int localSource, Int target )
+    EL_NO_RELEASE_EXCEPT;
+    void ProcessQueues();
+    void ProcessLocalQueues();
+
+    // For manually modifying/accessing buffers
+    void ForceNumLocalEdges( Int numLocalEdges );
+    void ForceConsistency( bool consistent=true ) EL_NO_EXCEPT;
+    Int* SourceBuffer() EL_NO_EXCEPT;
+    Int* TargetBuffer() EL_NO_EXCEPT;
+    Int* OffsetBuffer() EL_NO_EXCEPT;
+    const Int* LockedSourceBuffer() const EL_NO_EXCEPT;
+    const Int* LockedTargetBuffer() const EL_NO_EXCEPT;
+    const Int* LockedOffsetBuffer() const EL_NO_EXCEPT;
+    void ComputeSourceOffsets();
+
+    // Queries
+    // =======
+
+    // High-level data
+    // ---------------
+    Int NumSources() const EL_NO_EXCEPT;
+    Int NumTargets() const EL_NO_EXCEPT;
+    Int NumEdges() const EL_NO_EXCEPT;
+    Int FirstLocalSource() const EL_NO_EXCEPT;
+    Int NumLocalSources() const EL_NO_EXCEPT;
+    Int NumLocalEdges() const EL_NO_EXCEPT;
+    Int Capacity() const EL_NO_EXCEPT;
+    bool LocallyConsistent() const EL_NO_EXCEPT;
+
+    // Distribution information
+    // ------------------------
+    mpi::Comm Comm() const EL_NO_EXCEPT;
+    Int Blocksize() const EL_NO_EXCEPT;
+    int SourceOwner( Int s ) const EL_NO_RELEASE_EXCEPT;
+    Int GlobalSource( Int sLoc ) const EL_NO_RELEASE_EXCEPT;
+    Int LocalSource( Int s ) const EL_NO_RELEASE_EXCEPT;
+
+    // Detailed local information
+    // --------------------------
+    Int Source( Int localEdge ) const EL_NO_RELEASE_EXCEPT;
+    Int Target( Int localEdge ) const EL_NO_RELEASE_EXCEPT;
+    Int SourceOffset( Int localSource ) const EL_NO_RELEASE_EXCEPT;
+    Int Offset( Int localSource, Int target ) const EL_NO_RELEASE_EXCEPT;
+    Int NumConnections( Int localSource ) const EL_NO_RELEASE_EXCEPT;
+
+    // Return the ratio of the maximum number of local edges to the 
+    // total number of edges divided by the number of processes
+    double Imbalance() const EL_NO_RELEASE_EXCEPT;
+
+    mutable DistGraphMultMeta multMeta;
+    DistGraphMultMeta InitializeMultMeta() const;
+
+
+    void AssertConsistent() const;
+    void AssertLocallyConsistent() const;
+
+private:
+    Int numSources_, numTargets_;
+    mpi::Comm comm_;
+    // Apparently calling MPI_Comm_size in an inner loop is a very bad idea...
+    int commSize_;
+    int commRank_;
+
+    Int blocksize_;
+    Int numLocalSources_;
+
+    bool frozenSparsity_ = false;
+    vector<Int> sources_, targets_;
+    set<pair<Int,Int>> markedForRemoval_;
+
+    vector<Int> remoteSources_, remoteTargets_;
+    vector<pair<Int,Int>> remoteRemovals_;
+
+    void InitializeLocalData();
+
+    // Helpers for local indexing
+    bool locallyConsistent_ = true;
+    vector<Int> localSourceOffsets_;
+
+    friend class Graph;
+    friend void Copy( const Graph& A, DistGraph& B );
+    friend void Copy( const DistGraph& A, Graph& B );
+    friend void Copy( const DistGraph& A, DistGraph& B );
+
+    template<typename F> friend class DistSparseMatrix;
+};
+
+} // namespace El
 
 #endif // ifndef EL_CORE_DISTGRAPH_HPP