Graph tools

CMakeLists.txt

@@ -44,6 +44,7 @@ add_library(linalgcpp
     src/solvers.cpp
     src/sparsematrix.cpp
     src/timer.cpp
+    src/utilities.cpp
     src/vector.cpp
     src/vectorview.cpp
 )
@@ -58,10 +59,6 @@ target_include_directories(linalgcpp
         ${CMAKE_CURRENT_SOURCE_DIR}/src
 )
 
-add_subdirectory(tests)
-add_subdirectory(modules)
-add_subdirectory(examples)
-
 target_link_libraries(linalgcpp
     PRIVATE ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES}
 )
@@ -109,5 +106,8 @@ export(EXPORT linalgcpp-targets
        FILE ${CMAKE_CURRENT_BINARY_DIR}/linalgcppTargets.cmake
        NAMESPACE linalgcpp::
        )
-
 export(PACKAGE linalgcpp)
+
+add_subdirectory(tests)
+add_subdirectory(modules)
+add_subdirectory(examples)

config/extern_build.sh

@@ -42,6 +42,7 @@ CC=mpicc CXX=mpic++ cmake \
     -DLINALGCPP_ENABLE_MPI=Yes \
     -DLINALGCPP_ENABLE_METIS=Yes \
     -DLINALGCPP_ENABLE_SUITESPARSE=Yes \
+    -DLINALGCPP_ENABLE_GRAPHTOOLS=Yes \
     -DHypre_INC_DIR=$HYPRE_DIR/include \
     -DHypre_LIB_DIR=$HYPRE_DIR/lib \
     -DMETIS_DIR=$METIS_DIR \

config/gelever.sh

@@ -11,6 +11,7 @@ CXX=mpic++ CC=mpicc cmake .. \
     -DLINALGCPP_ENABLE_MPI=Yes \
     -DLINALGCPP_ENABLE_METIS=Yes \
     -DLINALGCPP_ENABLE_SUITESPARSE=Yes \
+    -DLINALGCPP_ENABLE_GRAPHTOOLS=Yes \
     -DHypre_DIR=${HYPRE_DIR} \
     -DSUITESPARSE_INCLUDE_DIR_HINTS=${SUITESPARSE_DIR}/include \
     -DSUITESPARSE_LIBRARY_DIR_HINTS=${SUITESPARSE_DIR}/lib

include/coomatrix.hpp

@@ -246,8 +246,11 @@ class CooMatrix : public Operator
         */
         void EliminateZeros(double tolerance = 0);
 
-    private:
+        /*! @brief Obtain the current size from maximum elements
+            @returns (rows, cols) current size
+        */
         std::tuple<int, int> FindSize() const;
+    private:
 
         bool size_set_;
 
@@ -707,8 +710,8 @@ std::tuple<int, int> CooMatrix<T>::FindSize() const
         return std::tuple<int, int> {rows_, cols_};
     }
 
-    int rows = 0;
-    int cols = 0;
+    int rows = -1;
+    int cols = -1;
 
     for (const auto& entry : entries_)
     {

include/densematrix.hpp

@@ -498,6 +498,19 @@ class DenseMatrix : public Operator
             @param dense dense matrix that holds the submatrix
         */
         void AddSubMatrix(const std::vector<int>& rows, std::vector<int>& cols, const DenseMatrix& dense);
+        /*! @brief Get a non-contigious submatrix, given the rows and cols
+            @param rows rows to get
+            @param cols cols to get
+            @param dense dense matrix that holds the submatrix
+        */
+        void GetSubMatrix(const std::vector<int>& rows, std::vector<int>& cols, DenseMatrix& dense) const;
+
+        /*! @brief Get a non-contigious submatrix, given the rows and cols
+            @param rows rows to get
+            @param cols cols to get
+            @returns dense dense matrix that holds the submatrix
+        */
+        DenseMatrix GetSubMatrix(const std::vector<int>& rows, std::vector<int>& cols) const;
 
         /*! @brief Compute singular values and vectors A = U * S * VT
                    Where S is returned and A is replaced with VT

include/eigensolver.hpp

@@ -77,6 +77,13 @@ class EigenSolver
         */
         void Solve(const DenseMatrix& mat, double rel_tol, int max_evect,
                    EigenPair& eigen_pair);
+
+        /*! @brief Compute full spectrum of eigen pairs
+
+            @param mat Dense matrix to solve for
+            @param EigenPair pair of eigenvalues and eigenvectors
+        */
+        void Solve(const DenseMatrix& mat, EigenPair& eigen_pair);
 
 
     private:

include/utilities.hpp

@@ -5,6 +5,7 @@
 
 #include <memory>
 #include "sparsematrix.hpp"
+#include "coomatrix.hpp"
 
 #if __cplusplus > 201103L
 using std::make_unique;
@@ -19,6 +20,85 @@ std::unique_ptr<T> make_unique(Ts&& ... params)
 namespace linalgcpp
 {
 
+/** @brief Set a global marker with given local indices
+    Such that marker[global_index] = local_index and
+    all other entries are -1
+
+    @param marker global marker to set
+    @param indices local indices
+*/
+void SetMarker(std::vector<int>& marker, const std::vector<int>& indices);
+
+/** @brief Clear a global marker with given local indices
+    Such that marker[global_index] = -1
+
+    @param marker global marker to clear
+    @param indices local indices
+*/
+void ClearMarker(std::vector<int>& marker, const std::vector<int>& indices);
+
+
+/** @brief Sparse identity of given size
+    @param size square size of identity
+    @return identity matrix
+*/
+template <typename T = double>
+SparseMatrix<T> SparseIdentity(int size);
+
+/** @brief Construct an rectangular identity matrix (as a SparseMatrix)
+    @param rows number of row
+    @param cols number of columns
+    @param row_offset offset row where diagonal identity starts
+    @param col_offset offset column where diagonal identity starts
+*/
+template <typename T = double>
+SparseMatrix<T> SparseIdentity(int rows, int cols, int row_offset = 0, int col_offset = 0);
+
+/** @brief Adds two sparse matrices C = alpha * A + beta * B
+           Nonzero entries do not have to match between A and B
+
+    @param alpha scale for A
+    @param A A matrix
+    @param beta scale for B
+    @param B B matrix
+    @returns C such that C = alpha * A + beta * B
+*/
+template <typename T = double>
+SparseMatrix<T> Add(double alpha, const SparseMatrix<T>& A,
+                    double beta, const SparseMatrix<T>& B);
+
+/** @brief Adds two sparse matrices C = A + B
+           Nonzero entries do not have to match between A and B
+
+    @param A A matrix
+    @param B B matrix
+    @returns C such that C = A + B
+*/
+template <typename T = double>
+SparseMatrix<T> Add(const SparseMatrix<T>& A,
+                    const SparseMatrix<T>& B);
+
+/** @brief Change the type of sparse matrix through copy
+
+    @param input Input matrix
+    @param output Output matrix of desired output type
+    @returns C such that C = A + B
+*/
+template <typename T, typename U>
+SparseMatrix<U> Duplicate(const SparseMatrix<T>& input);
+
+/** @brief Check if each entry is finite
+
+    @param x collection to check
+    @throws if any element is not finite
+*/
+template <typename T>
+void VerifyFinite(const T& x);
+
+///
+/// Inline implementations:
+///
+
 /*! @brief Throw if false in debug mode only */
 inline
 void linalgcpp_assert(bool expression, const std::string& message = "linalgcpp assertion failed")
@@ -63,31 +143,15 @@ void linalgcpp_verify(F&& lambda, const std::string& message = "linalgcpp verifi
     }
 }
 
-/** @brief Sparse identity of given size
-    @param size square size of identity
-    @return identity matrix
-*/
 template <typename T>
-SparseMatrix<T> SparseIdentity(int size);
-
-/** @brief Construct an rectangular identity matrix (as a SparseMatrix)
-    @param rows number of row
-    @param cols number of columns
-    @param row_offset offset row where diagonal identity starts
-    @param col_offset offset column where diagonal identity starts
-*/
-template <typename T>
-SparseMatrix<T> SparseIdentity(int rows, int cols, int row_offset = 0, int col_offset = 0);
-
-template <typename T = double>
 SparseMatrix<T> SparseIdentity(int size)
 {
     assert(size >= 0);
 
     return SparseMatrix<T>(std::vector<T>(size, (T)1.0));
 }
 
-template <typename T = double>
+template <typename T>
 SparseMatrix<T> SparseIdentity(int rows, int cols, int row_offset, int col_offset)
 {
     assert(rows >= 0);
@@ -113,6 +177,68 @@ SparseMatrix<T> SparseIdentity(int rows, int cols, int row_offset, int col_offse
     return SparseMatrix<T>(std::move(indptr), std::move(indices), std::move(data), rows, cols);
 }
 
+template <typename T>
+SparseMatrix<T> Add(const SparseMatrix<T>& A,
+                    const SparseMatrix<T>& B)
+{
+    return Add(1.0, A, 1.0, B);
+}
+
+template <typename T>
+SparseMatrix<T> Add(double alpha, const SparseMatrix<T>& A,
+                    double beta, const SparseMatrix<T>& B)
+{
+    assert(A.Rows() == B.Rows());
+    assert(A.Cols() == B.Cols());
+
+    CooMatrix<T> coo(A.Rows(), A.Cols());
+
+    auto add_mat = [&coo](double scale, const SparseMatrix<T> & mat)
+    {
+        const auto& indptr = mat.GetIndptr();
+        const auto& indices = mat.GetIndices();
+        const auto& data = mat.GetData();
+
+        int rows = mat.Rows();
+
+        for (int i = 0; i < rows; ++i)
+        {
+            for (int j = indptr[i]; j < indptr[i + 1]; ++j)
+            {
+                coo.Add(i, indices[j], scale * data[j]);
+            }
+        }
+    };
+
+    add_mat(alpha, A);
+    add_mat(beta, B);
+
+    return coo.ToSparse();
+}
+
+template <typename T, typename U>
+SparseMatrix<U> Duplicate(const SparseMatrix<T>& input)
+{
+    auto indptr = input.GetIndptr();
+    auto indices = input.GetIndices();
+    std::vector<U> data(input.GetData().size());
+
+    const auto& input_data = input.GetData();
+    std::copy(std::begin(input_data), std::end(input_data), std::begin(data));
+
+    return SparseMatrix<U>(std::move(indptr), std::move(indices), std::move(data),
+                           input.Rows(), input.Cols());
+}
+
+template <typename T>
+void VerifyFinite(const T& x)
+{
+    for (auto&& x_i : x)
+    {
+        linalgcpp_verify(std::isfinite(x_i), "X is not finite!");
+    }
+}
+
 } //namespace linalgcpp
 
 #endif // UTILITIES_HPP__

modules/CMakeLists.txt

@@ -10,4 +10,6 @@ if (LINALGCPP_ENABLE_METIS)
     add_subdirectory(partition)
 endif(LINALGCPP_ENABLE_METIS)
 
-
+if (LINALGCPP_ENABLE_GRAPHTOOLS)
+    add_subdirectory(graphtools)
+endif(LINALGCPP_ENABLE_GRAPHTOOLS)