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/decl.hpp

@@ -15,6 +15,45 @@
 
 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
@@ -134,6 +173,10 @@ class DistGraph
     // 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;
 
@@ -170,4 +213,4 @@ class DistGraph
 
 } // namespace El
 
-#endif // ifndef EL_CORE_DISTGRAPH_DECL_HPP
+#endif // ifndef EL_CORE_DISTGRAPH_HPP

include/El/core/DistSparseMatrix/decl.hpp

@@ -11,44 +11,8 @@
 #ifndef EL_CORE_DISTSPARSEMATRIX_DECL_HPP
 #define EL_CORE_DISTSPARSEMATRIX_DECL_HPP
 
-namespace El {
+namespace El  {
 
-struct DistSparseMultMeta
-{
-    bool ready;
-    // NOTE: The 'send' and 'recv' roles reverse for adjoint multiplication
-    Int numRecvInds;
-    vector<int> sendSizes, sendOffs,
-                recvSizes, recvOffs;
-    vector<Int> sendInds, colOffs;
-
-    DistSparseMultMeta() : ready(false), numRecvInds(0) { }
-
-    void Clear()
-    {
-        ready = false;
-        numRecvInds = 0; 
-        SwapClear( sendSizes );
-        SwapClear( recvSizes );
-        SwapClear( sendOffs );
-        SwapClear( recvOffs );
-        SwapClear( sendInds );
-        SwapClear( colOffs );
-    }
-
-    const DistSparseMultMeta& operator=( const DistSparseMultMeta& 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;
-    }
-};
 
 // Use a simple 1d distribution where each process owns a fixed number of rows,
 //     if last process,  height - (commSize-1)*floor(height/commSize)
@@ -194,8 +158,7 @@ class DistSparseMatrix
     // total number of nonzeros divided by the number of processes
     double Imbalance() const EL_NO_RELEASE_EXCEPT;
 
-    mutable DistSparseMultMeta multMeta;
-    DistSparseMultMeta InitializeMultMeta() const;
+    DistGraphMultMeta InitializeMultMeta() const;
 
     void MappedSources
     ( const DistMap& reordering, vector<Int>& mappedSources ) const;

include/El/core/DistSparseMatrix/impl.hpp

@@ -19,6 +19,8 @@
    which can be found in the LICENSE file in the root directory, or at 
    http://opensource.org/licenses/BSD-2-Clause
 */
+#ifndef EL_CORE_DISTSPARSEMATRIX_IMPL_HPP
+#define EL_CORE_DISTSPARSEMATRIX_IMPL_HPP
 
 #include <El/blas_like/level1/Axpy.hpp>
 #include <El/blas_like/level1/Scale.hpp>
@@ -70,7 +72,7 @@ void DistSparseMatrix<T>::Empty( bool freeMemory )
         SwapClear( vals_ );
     else
         vals_.resize( 0 );
-    multMeta.Clear();
+    distGraph_.multMeta.Clear();
 
     SwapClear( remoteVals_ );
 }
@@ -199,7 +201,7 @@ EL_NO_RELEASE_EXCEPT
     {
         distGraph_.QueueLocalConnection( localRow, col );
         vals_.push_back( value );
-        multMeta.ready = false;
+        distGraph_.multMeta.ready = false;
     }
 }
 
@@ -234,7 +236,7 @@ EL_NO_RELEASE_EXCEPT
     else
     {
         distGraph_.QueueLocalDisconnection( localRow, col );
-        multMeta.ready = false;
+        distGraph_.multMeta.ready = false;
     }
 }
 
@@ -404,7 +406,6 @@ DistSparseMatrix<T>::operator=( const DistSparseMatrix<T>& A )
     distGraph_ = A.distGraph_;
     vals_ = A.vals_;
     remoteVals_ = A.remoteVals_;
-    multMeta = A.multMeta;
     return *this;
 }
 
@@ -716,97 +717,8 @@ void DistSparseMatrix<T>::AssertLocallyConsistent() const
     if( !LocallyConsistent() )
         LogicError("Distributed sparse matrix must be consistent");
 }
-
 template<typename T>
-DistSparseMultMeta DistSparseMatrix<T>::InitializeMultMeta() const
-{
-    DEBUG_CSE
-    if( multMeta.ready )
-        return multMeta;
-    mpi::Comm comm = Comm();
-    const int commSize = distGraph_.commSize_;
-    auto& meta = multMeta;
-
-    // Compute the set of row indices that we need from X in a normal
-    // multiply or update of Y in the adjoint case
-    const Int* colBuffer = LockedTargetBuffer();
-    const Int numLocalEntries = NumLocalEntries();
-    vector<ValueInt<Int>> uniqueCols(numLocalEntries);
-    for( Int e=0; e<numLocalEntries; ++e )
-        uniqueCols[e] = ValueInt<Int>{colBuffer[e],e};
-    std::sort( uniqueCols.begin(), uniqueCols.end(), ValueInt<Int>::Lesser );
-    meta.colOffs.resize(numLocalEntries);
-    {
-        Int uniqueOff=-1, lastUnique=-1;
-        for( Int e=0; e<numLocalEntries; ++e )    
-        {
-            if( lastUnique != uniqueCols[e].value )
-            {
-                ++uniqueOff;
-                lastUnique = uniqueCols[e].value;
-                uniqueCols[uniqueOff] = uniqueCols[e];
-            }
-            meta.colOffs[uniqueCols[e].index] = uniqueOff;
-        }
-        uniqueCols.resize( uniqueOff+1 );
-    }
-    const Int numRecvInds = uniqueCols.size();
-    meta.numRecvInds = numRecvInds;
-    vector<Int> recvInds( numRecvInds );
-    meta.recvSizes.clear();
-    meta.recvSizes.resize( commSize, 0 );
-    meta.recvOffs.resize( commSize );
-    Int vecBlocksize = Width() / commSize;
-    if( vecBlocksize*commSize < Width() || Width() == 0 ) 
-        ++vecBlocksize;
-
-    {
-        Int off=0, lastOff=0, qPrev=0;
-        for( ; off<numRecvInds; ++off )
-        {
-            const Int j = uniqueCols[off].value;
-            const int q = j / vecBlocksize;
-            while( qPrev != q )
-            {
-                meta.recvSizes[qPrev] = off - lastOff;
-                meta.recvOffs[qPrev+1] = off;
-
-                lastOff = off;
-                ++qPrev;
-            }
-            recvInds[off] = j;
-        }
-        while( qPrev != commSize-1 )
-        {
-            meta.recvSizes[qPrev] = off - lastOff;
-            meta.recvOffs[qPrev+1] = off;
-            lastOff = off;
-            ++qPrev;
-        }
-        meta.recvSizes[commSize-1] = off - lastOff;
-    }
-
-    // Coordinate
-    meta.sendSizes.resize( commSize );
-    mpi::AllToAll( meta.recvSizes.data(), 1, meta.sendSizes.data(), 1, comm );
-    Int numSendInds=0;
-    meta.sendOffs.resize( commSize );
-    for( int q=0; q<commSize; ++q )
-    {
-        meta.sendOffs[q] = numSendInds;
-        numSendInds += meta.sendSizes[q];
-    }
-    meta.sendInds.resize( numSendInds );
-    mpi::AllToAll
-    ( recvInds.data(),      meta.recvSizes.data(), meta.recvOffs.data(),
-      meta.sendInds.data(), meta.sendSizes.data(), meta.sendOffs.data(),
-      comm );
-
-    meta.numRecvInds = numRecvInds;
-    meta.ready = true;
-
-    return meta;
-}
+DistGraphMultMeta DistSparseMatrix<T>::InitializeMultMeta() const {  return distGraph_.InitializeMultMeta();  }
 
 template<typename T>
 void DistSparseMatrix<T>::MappedSources
@@ -908,3 +820,4 @@ bool DistSparseMatrix<T>::CompareEntries( const Entry<T>& a, const Entry<T>& b )
 #undef EL_EXTERN
 
 } // namespace El
+#endif //EL_CORE_DISTSPARSEMATRIX_IMPL_HPP

src/blas_like/level1/ColumnMinAbs.cpp

@@ -187,7 +187,7 @@ void ColumnMinAbs
     mins.Resize( A.Width(), 1 );
     Fill( mins, limits::Max<Real>() );
     A.InitializeMultMeta();
-    const auto& meta = A.multMeta;
+    const auto& meta = A.LockedDistGraph().multMeta;
 
     // Pack the send values 
     // --------------------
@@ -243,7 +243,7 @@ void ColumnMinAbsNonzero
     // =========================================================
     mins = upperBounds;
     A.InitializeMultMeta();
-    const auto& meta = A.multMeta;
+    const auto& meta = A.LockedDistGraph().multMeta;
 
     // Pack the send values 
     // --------------------

src/blas_like/level1/ColumnNorms.cpp

@@ -191,8 +191,7 @@ void ColumnTwoNorms
     norms.Resize( A.Width(), 1 );
     Zero( norms );
     A.InitializeMultMeta();
-    const auto& meta = A.multMeta;
-
+    const auto& meta = A.LockedDistGraph().multMeta;
     // Pack the send values 
     // --------------------
     vector<Real> sendScales( meta.numRecvInds, 0 ),
@@ -263,7 +262,7 @@ void ColumnMaxNorms
     norms.Resize( A.Width(), 1 );
     Zero( norms );
     A.InitializeMultMeta();
-    const auto& meta = A.multMeta;
+    const auto& meta = A.LockedDistGraph().multMeta;
 
     // Pack the send values 
     // --------------------

src/blas_like/level1/Copy.cpp

@@ -71,6 +71,7 @@ void Copy( const DistGraph& A, DistGraph& B )
     // Directly assign instead of queueing up the individual edges
     B.sources_ = A.sources_;
     B.targets_ = A.targets_;
+    B.multMeta = A.multMeta;
     B.locallyConsistent_ = A.locallyConsistent_;
     B.localSourceOffsets_ = A.localSourceOffsets_;
     B.ProcessLocalQueues();