diff --git a/examples/dense/CMakeLists.txt b/examples/dense/CMakeLists.txt index cf71d13..18eefc7 100644 --- a/examples/dense/CMakeLists.txt +++ b/examples/dense/CMakeLists.txt @@ -7,6 +7,10 @@ target_link_libraries(dense_lu linalgcpp) add_executable(dense_qr dense_qr.cpp) target_link_libraries(dense_qr linalgcpp) +add_test(example_dense_basics dense_basics) +add_test(example_dense_lu dense_lu) +add_test(example_dense_qr dense_qr) + function(copy_file filename) configure_file( "${PROJECT_SOURCE_DIR}/examples/dense/${filename}" diff --git a/examples/dense/dense_lu.cpp b/examples/dense/dense_lu.cpp index 3a98126..4eb0a82 100644 --- a/examples/dense/dense_lu.cpp +++ b/examples/dense/dense_lu.cpp @@ -125,7 +125,7 @@ void LU() void time_LU() { - for (int i = 10; i < 2000; i *= 1.5) + for (int i = 10; i < 800; i *= 1.5) { std::cout << "LU " << i << "\t"; diff --git a/examples/dense/dense_qr.cpp b/examples/dense/dense_qr.cpp index fa0ee30..9d4fdf4 100644 --- a/examples/dense/dense_qr.cpp +++ b/examples/dense/dense_qr.cpp @@ -86,7 +86,7 @@ double test_QR(int m, int n) void time_QR() { - for (int i = 10; i < 2000; i *= 1.5) + for (int i = 10; i < 800; i *= 1.5) { int m = i; int n = i * 1.5; @@ -152,7 +152,7 @@ void QR_eig() } // Run QR - for (int i = 0; i < 10000; ++i) + for (int i = 0; i < 1000; ++i) { auto qr = qr_decomp(A); auto& Q = qr.first; diff --git a/examples/include/ex_utilities.hpp b/examples/include/ex_utilities.hpp index c8c49c1..7425901 100644 --- a/examples/include/ex_utilities.hpp +++ b/examples/include/ex_utilities.hpp @@ -3,12 +3,7 @@ #include -#include "parser.hpp" -#include "vector.hpp" -#include "densematrix.hpp" -#include "sparsematrix.hpp" -#include "eigensolver.hpp" -#include "timer.hpp" +#include "linalgcpp.hpp" using Vector = linalgcpp::Vector; using VectorView = linalgcpp::VectorView; @@ -17,6 +12,8 @@ using SparseMatrix = linalgcpp::SparseMatrix; using CooMatrix = linalgcpp::CooMatrix; using Timer = linalgcpp::Timer; using Operator = linalgcpp::Operator; +using linalgcpp::linalgcpp_assert; +using linalgcpp::linalgcpp_verify; inline diff --git a/examples/sparse/CMakeLists.txt b/examples/sparse/CMakeLists.txt index d3fcf8f..2cfc2ef 100644 --- a/examples/sparse/CMakeLists.txt +++ b/examples/sparse/CMakeLists.txt @@ -1,6 +1,8 @@ add_executable(sparse_basics sparse_basics.cpp) target_link_libraries(sparse_basics linalgcpp) +add_test(example_sparse_basics sparse_basics) + function(copy_file filename) configure_file( "${PROJECT_SOURCE_DIR}/examples/sparse/${filename}" diff --git a/examples/vector/CMakeLists.txt b/examples/vector/CMakeLists.txt index af6be49..2c2b084 100644 --- a/examples/vector/CMakeLists.txt +++ b/examples/vector/CMakeLists.txt @@ -1,6 +1,8 @@ add_executable(vector_ops vector_ops.cpp) target_link_libraries(vector_ops linalgcpp) +add_test(example_vector_basics vector_ops) + function(copy_file filename) configure_file( "${PROJECT_SOURCE_DIR}/examples/vector/${filename}" diff --git a/include/blockoperator.hpp b/include/blockoperator.hpp index c5fbf1a..e837649 100644 --- a/include/blockoperator.hpp +++ b/include/blockoperator.hpp @@ -43,6 +43,9 @@ class BlockOperator : public Operator /*! @brief Move deconstructor */ BlockOperator(BlockOperator&& other) noexcept; + /*! @brief Assignment operator */ + BlockOperator& operator=(BlockOperator&& other) noexcept; + /*! @brief Default deconstructor */ ~BlockOperator() noexcept = default; @@ -62,14 +65,14 @@ class BlockOperator : public Operator /*! @brief Get a block @param i row index @param j column index - @retval SparseMatrix at index (i, j) + @retval op Operator at index (i, j) */ const Operator* GetBlock(int i, int j) const; /*! @brief Set a block @param i row index @param j column index - @param SparseMatrix at index (i, j) + @param op Operator at index (i, j) */ void SetBlock(int i, int j, const Operator& op); @@ -97,6 +100,106 @@ class BlockOperator : public Operator mutable Vector tmp_; }; +/*! @brief Diagonal ense block matrix where each block is a sparse matrix + Blocks of zero size are assumed to be zero matrix + Offsets are to be of size (blocks + 1), where the last + entry is the total size. + + @Note This is a view. BlockDiagOperator does not take + possession of any given input! +*/ +class BlockDiagOperator : public Operator +{ + public: + /*! @brief Default Constructor of zero size */ + BlockDiagOperator(); + + /*! @brief Square Constructor with given symmetric offsets*/ + explicit BlockDiagOperator(std::vector offsets); + + /*! @brief Rectangle Constructor with given offsets*/ + BlockDiagOperator(std::vector row_offsets, std::vector col_offsets); + + /*! @brief Copy deconstructor */ + BlockDiagOperator(const BlockDiagOperator& other) noexcept; + + /*! @brief Move deconstructor */ + BlockDiagOperator(BlockDiagOperator&& other) noexcept; + + /*! @brief Assignment operator */ + BlockDiagOperator& operator=(BlockDiagOperator&& other) noexcept; + + /*! @brief Default deconstructor */ + ~BlockDiagOperator() noexcept = default; + + /*! @brief Swap two block operators */ + friend void swap(BlockDiagOperator& lhs, BlockDiagOperator& rhs) noexcept; + + /*! @brief Get the row offsets + @retval the row offsets + */ + const std::vector& GetRowOffsets() const; + + /*! @brief Get the col offsets + @retval the col offsets + */ + const std::vector& GetColOffsets() const; + + /*! @brief Get a block + @param i row, column index + @retval op Operator at index (i, i) + + */ + const Operator* GetBlock(int i) const; + + /*! @brief Get a block + @param i row index + @param j column index + @retval op Operator at index (i, j) + + @note i must equal j + */ + const Operator* GetBlock(int i, int j) const; + + /*! @brief Set a block + @param i diagonal row, col index + @param op Operator at index (i, j) + */ + void SetBlock(int i, const Operator& op); + + /*! @brief Set a block + @param i row index + @param j column index + @param op Operator at index (i, j) + + @note i must equal j + */ + void SetBlock(int i, int j, const Operator& op); + + /*! @brief Get the transpose matrix + */ + BlockDiagOperator Transpose() const; + + /// Operator Requirement + void Mult(const VectorView& input, VectorView output) const override; + + /// Operator Requirement + void MultAT(const VectorView& input, VectorView output) const override; + + using Operator::Mult; + + private: + std::vector row_offsets_; + std::vector col_offsets_; + + std::vector A_; + + mutable BlockVector x_; + mutable BlockVector y_; + + mutable Vector tmp_; +}; + } //namespace linalgcpp #endif // BLOCKOPERATOR_HPP__ diff --git a/include/linalgcpp.hpp b/include/linalgcpp.hpp index c7ba71e..647b743 100644 --- a/include/linalgcpp.hpp +++ b/include/linalgcpp.hpp @@ -17,3 +17,4 @@ #include "eigensolver.hpp" #include "timer.hpp" #include "argparser.hpp" +#include "utilities.hpp" diff --git a/include/solvers.hpp b/include/solvers.hpp index 63ee2d1..a95d7e5 100644 --- a/include/solvers.hpp +++ b/include/solvers.hpp @@ -18,10 +18,13 @@ class Solver : public Operator Solver(); /*! @brief Copy Constructor */ - Solver(const Solver& other) noexcept = default; + Solver(const Solver& other) noexcept; /*! @brief Move Constructor */ - Solver(Solver&& other) noexcept = default; + Solver(Solver&& other) noexcept; + + /*! @brief Assignement Operator */ + Solver& operator=(Solver&& other) noexcept; /*! @brief Destructor */ virtual ~Solver() noexcept = default; diff --git a/include/sparsematrix.hpp b/include/sparsematrix.hpp index 8a79419..732f878 100644 --- a/include/sparsematrix.hpp +++ b/include/sparsematrix.hpp @@ -424,8 +424,9 @@ class SparseMatrix : public Operator /*! @brief Remove entries less than a tolerance @param tol tolerance to remove + @param keep_diag keep diagonal element, regardless of value */ - void EliminateZeros(T tolerance = 0); + void EliminateZeros(T tolerance = 0, bool keep_diag = false); /// Operator Requirement, calls the templated Mult void Mult(const VectorView& input, VectorView output) const override; @@ -1470,7 +1471,7 @@ void SparseMatrix::EliminateCol(const U& marker, const VectorView& x, } template -void SparseMatrix::EliminateZeros(T tolerance) +void SparseMatrix::EliminateZeros(T tolerance, bool keep_diag) { std::vector elim_indptr(1, 0); std::vector elim_indices; @@ -1484,7 +1485,10 @@ void SparseMatrix::EliminateZeros(T tolerance) { for (int j = indptr_[i]; j < indptr_[i + 1]; ++j) { - if (std::fabs(data_[j]) > tolerance) + int col = indices_[j]; + + if (std::fabs(data_[j]) > tolerance || + (keep_diag && col == i)) { elim_indices.push_back(indices_[j]); elim_data.push_back(data_[j]); diff --git a/include/utilities.hpp b/include/utilities.hpp index cde6cb1..4c23f4b 100644 --- a/include/utilities.hpp +++ b/include/utilities.hpp @@ -3,6 +3,19 @@ #ifndef UTILITIES_HPP__ #define UTILITIES_HPP__ +#include +#include "sparsematrix.hpp" + +#if __cplusplus > 201103L +using std::make_unique; +#else +template +std::unique_ptr make_unique(Ts&& ... params) +{ + return std::unique_ptr(new T(std::forward(params)...)); +} +#endif + namespace linalgcpp { @@ -10,7 +23,7 @@ namespace linalgcpp inline void linalgcpp_assert(bool expression, const std::string& message = "linalgcpp assertion failed") { -#ifdef NDEBUG +#ifndef NDEBUG if (!expression) { throw std::runtime_error(message); @@ -28,6 +41,78 @@ void linalgcpp_verify(bool expression, const std::string& message = "linalgcpp v } } +/*! @brief Throw if false in debug mode only */ +template +void linalgcpp_assert(F&& lambda, const std::string& message = "linalgcpp assertion failed") +{ +#ifndef NDEBUG + if (!lambda()) + { + throw std::runtime_error(message); + } +#endif +} + +/*! @brief Throw if false unconditionally */ +template +void linalgcpp_verify(F&& lambda, const std::string& message = "linalgcpp verification failed") +{ + if (!lambda()) + { + throw std::runtime_error(message); + } +} + +/** @brief Sparse identity of given size + @param size square size of identity + @return identity matrix +*/ +template +SparseMatrix 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 +SparseMatrix SparseIdentity(int rows, int cols, int row_offset = 0, int col_offset = 0); + +template +SparseMatrix SparseIdentity(int size) +{ + assert(size >= 0); + + return SparseMatrix(std::vector(size, (T)1.0)); +} + +template +SparseMatrix SparseIdentity(int rows, int cols, int row_offset, int col_offset) +{ + assert(rows >= 0); + assert(cols >= 0); + assert(row_offset <= rows); + assert(row_offset >= 0); + assert(col_offset <= cols); + assert(col_offset >= 0); + + const int diag_size = std::min(rows - row_offset, cols - col_offset); + + std::vector indptr(rows + 1); + + std::fill(std::begin(indptr), std::begin(indptr) + row_offset, 0); + std::iota(std::begin(indptr) + row_offset, std::begin(indptr) + row_offset + diag_size, 0); + std::fill(std::begin(indptr) + row_offset + diag_size, std::begin(indptr) + rows + 1, diag_size); + + std::vector indices(diag_size); + std::iota(std::begin(indices), std::begin(indices) + diag_size, col_offset); + + std::vector data(diag_size, 1.0); + + return SparseMatrix(std::move(indptr), std::move(indices), std::move(data), rows, cols); +} + } //namespace linalgcpp #endif // UTILITIES_HPP__ diff --git a/include/vectorview.hpp b/include/vectorview.hpp index 0a7e721..389264f 100644 --- a/include/vectorview.hpp +++ b/include/vectorview.hpp @@ -87,6 +87,26 @@ class VectorView */ virtual const T* end() const; + /*! @brief STL like front. Points to first element of data + @retval reference to the first element of data + */ + virtual T& front(); + + /*! @brief STL like front. Points to first element of data + @retval reference to the first element of data + */ + virtual const T& front() const; + + /*! @brief STL like back. Points to last element of data + @retval reference to the last element of data + */ + virtual T& back(); + + /*! @brief STL like back. Points to last element of data + @retval reference to the last element of data + */ + virtual const T& back() const; + /*! @brief Index operator @param i index into vector @retval reference to value at index i @@ -316,6 +336,34 @@ const T* VectorView::end() const return data_ + size_; } +template +T& VectorView::front() +{ + assert(size_ > 0); + return data_[0]; +} + +template +const T& VectorView::front() const +{ + assert(size_ > 0); + return data_[0]; +} + +template +T& VectorView::back() +{ + assert(size_ > 0); + return data_[size_ - 1]; +} + +template +const T& VectorView::back() const +{ + assert(size_ > 0); + return data_[size_ - 1]; +} + template T& VectorView::operator[](int i) { diff --git a/modules/parlinalgcpp/CMakeLists.txt b/modules/parlinalgcpp/CMakeLists.txt index 9f92c94..1037790 100644 --- a/modules/parlinalgcpp/CMakeLists.txt +++ b/modules/parlinalgcpp/CMakeLists.txt @@ -46,12 +46,6 @@ add_library(parlinalgcpp add_dependencies(parlinalgcpp linalgcpp) -set_target_properties(parlinalgcpp PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - target_compile_options(parlinalgcpp PRIVATE $<$:-Wall -O2>) target_compile_options(parlinalgcpp PRIVATE $<$:-Wall -O2>) diff --git a/modules/parlinalgcpp/include/parcommpkg.hpp b/modules/parlinalgcpp/include/parcommpkg.hpp index 02d50f1..9c4bbaf 100644 --- a/modules/parlinalgcpp/include/parcommpkg.hpp +++ b/modules/parlinalgcpp/include/parcommpkg.hpp @@ -7,6 +7,7 @@ #include #include +#include "utilities.hpp" #include "_hypre_parcsr_mv.h" namespace linalgcpp diff --git a/modules/parlinalgcpp/include/parmatrix.hpp b/modules/parlinalgcpp/include/parmatrix.hpp index 193aa7f..4dedb46 100644 --- a/modules/parlinalgcpp/include/parmatrix.hpp +++ b/modules/parlinalgcpp/include/parmatrix.hpp @@ -205,6 +205,15 @@ class ParMatrix: public ParOperator /*! @brief Access the column map */ const std::vector& GetColMap() const; + /*! @brief Access the diagonal block */ + linalgcpp::SparseMatrix& GetDiag(); + + /*! @brief Access the off-diagonal block */ + linalgcpp::SparseMatrix& GetOffd(); + + /*! @brief Access the column map */ + std::vector& GetColMap(); + /*! @brief Print entries in both blocks and the column map @param label Label to show prior to printing @param out Output stream to print to @@ -276,6 +285,18 @@ class ParMatrix: public ParOperator */ void EliminateRow(int index); + /*! @brief Return the row number of nonzero elements combining both the + diagonal and off diagonal blocks. + @param row row to to get size of + */ + int RowSize(int row); + + /*! @brief Remove entries less than tolerance + @param tol tolerance + @param keep_diag keep diagonal element, regardless of value + */ + void EliminateZeros(double tol = 0.0, bool keep_diag = false); + private: void Init(); void HypreCreate(); @@ -307,6 +328,24 @@ const linalgcpp::SparseMatrix& ParMatrix::GetOffd() const return offd_; } +inline +std::vector& ParMatrix::GetColMap() +{ + return col_map_; +} + +inline +linalgcpp::SparseMatrix& ParMatrix::GetDiag() +{ + return diag_; +} + +inline +linalgcpp::SparseMatrix& ParMatrix::GetOffd() +{ + return offd_; +} + inline ParMatrix ParMatrix::operator*(const ParMatrix& other) const { diff --git a/modules/parlinalgcpp/include/paroperator.hpp b/modules/parlinalgcpp/include/paroperator.hpp index de6eb6e..c069285 100644 --- a/modules/parlinalgcpp/include/paroperator.hpp +++ b/modules/parlinalgcpp/include/paroperator.hpp @@ -66,6 +66,12 @@ class ParOperator: public linalgcpp::Operator virtual void Mult(const linalgcpp::VectorView& input, linalgcpp::VectorView output) const = 0; + /*! @brief Apply the action of this operator Ax = y + @param input vector x + @param output vector y + */ + virtual Vector Mult(const linalgcpp::VectorView& input) const; + /* @brief Global number of rows */ virtual int GlobalRows() const; diff --git a/modules/parlinalgcpp/include/parutilities.hpp b/modules/parlinalgcpp/include/parutilities.hpp index 3f02740..a6b39a0 100644 --- a/modules/parlinalgcpp/include/parutilities.hpp +++ b/modules/parlinalgcpp/include/parutilities.hpp @@ -12,11 +12,63 @@ namespace linalgcpp { +/// Call output only on processor 0 +#define ParPrint(myid, output) if (myid == 0) { output; } + +/*! @brief Throw if false in debug mode only */ +template +void linalgcpp_parassert(MPI_Comm comm, F&& lambda, const std::string& message = "linalgcpp assertion failed") +{ +#ifndef NDEBUG + int local_exp = !lambda(); + int global_exp; + MPI_Allreduce(&local_exp, &global_exp, 1, MPI_INT, MPI_SUM, comm); + + if (global_exp) + { + throw std::runtime_error(message); + } +#endif // NDEBUG +} + +/*! @brief Throw if false unconditionally */ +template +void linalgcpp_parverify(MPI_Comm comm, F&& lambda, const std::string& message = "linalgcpp verification failed") +{ + int local_exp = !lambda(); + int global_exp; + MPI_Allreduce(&local_exp, &global_exp, 1, MPI_INT, MPI_SUM, comm); + + if (global_exp) + { + throw std::runtime_error(message); + } +} + +/*! @brief Throw if false unconditionally */ +inline +void linalgcpp_parverify(MPI_Comm comm, bool expression, const std::string& message = "linalgcpp verification failed") +{ + int local_exp = !expression; + int global_exp; + MPI_Allreduce(&local_exp, &global_exp, 1, MPI_INT, MPI_SUM, comm); + + if (global_exp) + { + throw std::runtime_error(message); + } +} + + + class ParMatrix; /// Split a square matrix between processes ParMatrix ParSplit(MPI_Comm comm, const linalgcpp::SparseMatrix& A_global, const std::vector& proc_part); +/// Split a square matrix between processes, including local to global map +ParMatrix ParSplit(MPI_Comm comm, const linalgcpp::SparseMatrix& A_global, const std::vector& proc_part, std::vector& local_part); + /// Generate offsets given local sizes std::vector GenerateOffsets(MPI_Comm comm, int local_size); @@ -69,5 +121,36 @@ ParMatrix ParSub(double alpha, const ParMatrix& A, double beta, const ParMatrix& /** @note Row starts must match between A and B */ ParMatrix ParSub(double alpha, const ParMatrix& A, double beta, const ParMatrix& B, std::vector& marker); +/// Othogonalize against constant vector +void ParSubAvg(MPI_Comm comm, VectorView x); + +/// Othogonalize against constant vector, with known global size +void ParSubAvg(MPI_Comm comm, VectorView x, int global_size); + +/** @brief Handles mpi initialization and finalization */ +struct MpiSession +{ + /** @brief Constructor + + @param argc argc from command line + @param argv argv from command line + @param comm MPI Communicator to use + */ + MpiSession(int argc, char** argv, MPI_Comm comm_in = MPI_COMM_WORLD) + : comm(comm_in) + { + MPI_Init(&argc, &argv); + MPI_Comm_size(comm, &num_procs); + MPI_Comm_rank(comm, &myid); + } + + /** @brief Destructor */ + ~MpiSession() { MPI_Finalize(); } + + MPI_Comm comm; + int num_procs; + int myid; +}; + } // namespace linalgcpp #endif // PARUTILITIES_HPP__ diff --git a/modules/parlinalgcpp/src/parcommpkg.cpp b/modules/parlinalgcpp/src/parcommpkg.cpp index cf3e511..eab3a78 100644 --- a/modules/parlinalgcpp/src/parcommpkg.cpp +++ b/modules/parlinalgcpp/src/parcommpkg.cpp @@ -6,11 +6,11 @@ namespace linalgcpp ParCommPkg::ParCommPkg(const hypre_ParCSRMatrix* A) : comm_(A->comm) { - assert(A); + linalgcpp_assert(A != NULL); const hypre_ParCSRCommPkg* comm_pkg = A->comm_pkg; - assert(comm_pkg); + linalgcpp_assert(comm_pkg != NULL); num_sends_ = comm_pkg->num_sends; diff --git a/modules/parlinalgcpp/src/parmatrix.cpp b/modules/parlinalgcpp/src/parmatrix.cpp index 9ec0582..c8ab11d 100644 --- a/modules/parlinalgcpp/src/parmatrix.cpp +++ b/modules/parlinalgcpp/src/parmatrix.cpp @@ -346,7 +346,7 @@ void ParMatrix::AddDiag(const std::vector& diag_vals) ParCommPkg ParMatrix::MakeCommPkg() const { - assert(A_); + linalgcpp_assert(A_ != NULL, "ParMatrix::A_ not valid!"); return ParCommPkg(A_); } @@ -397,4 +397,17 @@ void ParMatrix::EliminateRow(int index) offd_.EliminateRow(index); } +int ParMatrix::RowSize(int row) +{ + return diag_.RowSize(row) + offd_.RowSize(row); +} + +void ParMatrix::EliminateZeros(double tol, bool keep_diag) +{ + diag_.EliminateZeros(tol, keep_diag); + offd_.EliminateZeros(tol); + + HypreCreate(); +} + } // namespace linalgcpp diff --git a/modules/parlinalgcpp/src/paroperator.cpp b/modules/parlinalgcpp/src/paroperator.cpp index 5bcb953..1cc8195 100644 --- a/modules/parlinalgcpp/src/paroperator.cpp +++ b/modules/parlinalgcpp/src/paroperator.cpp @@ -98,6 +98,15 @@ void swap(ParOperator& lhs, ParOperator& rhs) noexcept std::swap(lhs.b_, rhs.b_); } +Vector ParOperator::Mult(const linalgcpp::VectorView& input) const +{ + Vector output(Cols(), 0.0); + + Mult(input, output); + + return output; +} + int ParOperator::GlobalRows() const { int size = 0.0; diff --git a/modules/parlinalgcpp/src/parutilities.cpp b/modules/parlinalgcpp/src/parutilities.cpp index cdebca8..852e0d9 100644 --- a/modules/parlinalgcpp/src/parutilities.cpp +++ b/modules/parlinalgcpp/src/parutilities.cpp @@ -5,8 +5,15 @@ namespace linalgcpp { +ParMatrix ParSplit(MPI_Comm comm, const linalgcpp::SparseMatrix& A_global, const std::vector& proc_part) +{ + std::vector local_part; + return ParSplit(comm, A_global, proc_part, local_part); +} + + ParMatrix ParSplit(MPI_Comm comm, const linalgcpp::SparseMatrix& A_global, - const std::vector& proc_part) + const std::vector& proc_part, std::vector& local_part) { assert(A_global.Rows() == A_global.Cols()); @@ -19,7 +26,7 @@ ParMatrix ParSplit(MPI_Comm comm, const linalgcpp::SparseMatrix& A_globa linalgcpp::CooMatrix local_offd; std::vector col_map; - std::vector local_part; + local_part.clear(); const int proc_size = proc_part.size(); @@ -319,4 +326,25 @@ ParMatrix ParSub(double alpha, const ParMatrix& A, double beta, const ParMatrix& return ParAdd(alpha, A, -beta, B, marker); } +void ParSubAvg(MPI_Comm comm, VectorView x) +{ + int local_size = x.size(); + int global_size; + + MPI_Allreduce(&local_size, &global_size, 1, MPI_INT, MPI_SUM, comm); + + ParSubAvg(comm, x, global_size); +} + +void ParSubAvg(MPI_Comm comm, VectorView x, int global_size) +{ + double local_sum = Sum(x); + double global_sum; + + MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, comm); + global_sum /= global_size; + + x -= global_sum; +} + } //namsespace linalgcpp diff --git a/modules/parlinalgcpp/tests/CMakeLists.txt b/modules/parlinalgcpp/tests/CMakeLists.txt index 0bdcdfc..3ed90e3 100644 --- a/modules/parlinalgcpp/tests/CMakeLists.txt +++ b/modules/parlinalgcpp/tests/CMakeLists.txt @@ -1,12 +1,6 @@ add_executable(test_parlinalgcpp test_parlinalg.cpp) target_link_libraries(test_parlinalgcpp parlinalgcpp) -set_target_properties(test_parlinalgcpp PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - configure_file(${CMAKE_CURRENT_LIST_DIR}/test_data/test_mat.bin ${CMAKE_CURRENT_BINARY_DIR}/test_data/test_mat.bin COPYONLY) diff --git a/modules/parlinalgcpp/tests/test_parlinalg.cpp b/modules/parlinalgcpp/tests/test_parlinalg.cpp index 5773e16..4f185c4 100644 --- a/modules/parlinalgcpp/tests/test_parlinalg.cpp +++ b/modules/parlinalgcpp/tests/test_parlinalg.cpp @@ -91,7 +91,7 @@ int main(int argc, char** argv) test_pcg_copy = cg; // rofl - std::vector ops { + std::vector ops { &cg, &pcg_boomer, &pcg_parasails, @@ -136,7 +136,7 @@ int main(int argc, char** argv) for (const auto& op : ops) { if (myid == 0) - printf("Op: %d\n", count++); + printf("Op: %d\t", count++); for (auto& x_i : evects) { @@ -167,9 +167,14 @@ void PrintVector(MPI_Comm comm, const std::vector>& ev if (myid == 0) { std::cout.precision(10); - std::cout << "Evals: " << evals; + std::cout << "Evals: "; + for (auto&& ev : evals) + { + std::cout << ev << " "; + } + std::cout << "\n"; - int subset_size = std::min(10, evects[0].size()); + int subset_size = std::min(0, evects[0].size()); for (int j = 0; j < subset_size; ++j) { for(size_t i = 0; i < evects.size(); ++i) @@ -182,3 +187,27 @@ void PrintVector(MPI_Comm comm, const std::vector>& ev } } + +SparseMatrix make_agg_vertex(const ParMatrix& A) +{ + int num_aggs = 3; + int num_rows = A.Rows(); + + std::vector indptr(num_aggs + 1); + indptr[0] = 0; + //indptr[1] = num_rows; + + //indptr[1] = num_rows / 2; + //indptr[2] = num_rows; + + indptr[1] = num_rows / 3; + indptr[2] = 2 * num_rows / 3; + indptr[3] = num_rows; + + std::vector indices(num_rows); + std::iota(std::begin(indices), std::end(indices), 0); + std::vector data(num_rows, 1.0); + + return SparseMatrix(std::move(indptr), std::move(indices), std::move(data), + num_aggs, num_rows); +} diff --git a/modules/partition/include/partition.hpp b/modules/partition/include/partition.hpp index 04c895e..3155030 100644 --- a/modules/partition/include/partition.hpp +++ b/modules/partition/include/partition.hpp @@ -107,9 +107,9 @@ template std::vector Partition(const linalgcpp::SparseMatrix& mat, int num_parts = 2, double unbalance_factor = 1.0, bool contig = true, bool weighted = false) { - assert(num_parts > 0); - assert(num_parts <= mat.Rows()); - assert(mat.Cols() == mat.Rows()); + linalgcpp_assert(num_parts > 0, "Non-Positive Partition Parts"); + linalgcpp_assert(num_parts <= mat.Rows(), "Too Many Partition Parts"); + linalgcpp_assert(mat.Cols() == mat.Rows(), "Non-Square Partition Relationship"); const int size = mat.Rows(); @@ -141,8 +141,9 @@ std::vector Partition(const linalgcpp::SparseMatrix& mat, int num_parts idx_t objval; - METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, + int metis_status = METIS_PartGraphKway(&nvtxs, &ncon, xadj, adjncy, vwgt, vsize, adjwgt, &nparts, tpwgts, &ubvec, options, &objval, part.data()); + linalgcpp_verify(metis_status == METIS_OK, "METIS Failed to Partition!"); RemoveEmpty(part); diff --git a/modules/partition/tests/CMakeLists.txt b/modules/partition/tests/CMakeLists.txt index 1d416db..65c431e 100644 --- a/modules/partition/tests/CMakeLists.txt +++ b/modules/partition/tests/CMakeLists.txt @@ -4,11 +4,4 @@ target_link_libraries(test_partition partition) target_compile_options(test_partition PRIVATE $<$:-Wall -Wextra -O2>) target_compile_options(test_partition PRIVATE $<$:-Wall -Wextra -O2>) -set_target_properties(test_partition PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - - add_test(test_partition test_partition) diff --git a/modules/sparsesolver/CMakeLists.txt b/modules/sparsesolver/CMakeLists.txt index c49ec17..19bafcc 100644 --- a/modules/sparsesolver/CMakeLists.txt +++ b/modules/sparsesolver/CMakeLists.txt @@ -28,24 +28,15 @@ ############################################################### cmake_minimum_required(VERSION 3.5) -#project(sparsesolve VERSION 1.0.0 LANGUAGES CXX) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/modules") find_package(SuiteSparse REQUIRED) -#find_package(linalgcpp REQUIRED) add_library(sparsesolve src/sparsesolve.cpp) add_dependencies(sparsesolve linalgcpp) -set_target_properties(sparsesolve PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - target_link_libraries(sparsesolve - #PUBLIC SuiteSparse::SuiteSparse linalgcpp::linalgcpp PUBLIC SuiteSparse::SuiteSparse linalgcpp ) diff --git a/modules/sparsesolver/include/sparsesolve.hpp b/modules/sparsesolver/include/sparsesolve.hpp index 6ee7c64..7da8338 100644 --- a/modules/sparsesolver/include/sparsesolve.hpp +++ b/modules/sparsesolver/include/sparsesolve.hpp @@ -6,6 +6,9 @@ #include "umfpack.h" #include "linalgcpp.hpp" +namespace linalgcpp +{ + class SparseSolver : public linalgcpp::Operator { public: @@ -39,4 +42,6 @@ class SparseSolver : public linalgcpp::Operator mutable std::vector info_; }; +} // namespace linalgcpp + #endif // SPARSESOLVE_HPP__ diff --git a/modules/sparsesolver/src/sparsesolve.cpp b/modules/sparsesolver/src/sparsesolve.cpp index b354fd5..64f4aa2 100644 --- a/modules/sparsesolver/src/sparsesolve.cpp +++ b/modules/sparsesolver/src/sparsesolve.cpp @@ -2,6 +2,9 @@ #include "sparsesolve.hpp" +namespace linalgcpp +{ + SparseSolver::SparseSolver(linalgcpp::SparseMatrix A) : Operator(A), A_(std::move(A)), numeric_(nullptr), control_(UMFPACK_CONTROL, 0), info_(UMFPACK_INFO, 0) @@ -138,3 +141,4 @@ void SparseSolver::MultAT(const linalgcpp::VectorView& input, assert(status == 0); } +} // namespace linalgcpp diff --git a/modules/sparsesolver/tests/CMakeLists.txt b/modules/sparsesolver/tests/CMakeLists.txt index 5ba03fd..37214eb 100644 --- a/modules/sparsesolver/tests/CMakeLists.txt +++ b/modules/sparsesolver/tests/CMakeLists.txt @@ -1,12 +1,4 @@ add_executable(test_sparsesolver test_sparsesolver.cpp) target_link_libraries(test_sparsesolver sparsesolve) -set_target_properties(test_sparsesolver PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - add_test(test_sparsesolver test_sparsesolver) - - diff --git a/src/blockoperator.cpp b/src/blockoperator.cpp index a7cc0a4..971631e 100644 --- a/src/blockoperator.cpp +++ b/src/blockoperator.cpp @@ -42,6 +42,13 @@ BlockOperator::BlockOperator(BlockOperator&& other) noexcept swap(*this, other); } +BlockOperator& BlockOperator::operator=(BlockOperator&& other) noexcept +{ + swap(*this, other); + + return *this; +} + void swap(BlockOperator& lhs, BlockOperator& rhs) noexcept { swap(static_cast(lhs), static_cast(rhs)); @@ -149,4 +156,143 @@ void BlockOperator::MultAT(const VectorView& input, VectorView o output = x_; } +BlockDiagOperator::BlockDiagOperator() + : row_offsets_(1, 0), col_offsets_(1, 0) +{ + +} + +BlockDiagOperator::BlockDiagOperator(std::vector offsets) : + Operator(offsets.back()), + row_offsets_(offsets), col_offsets_(offsets), + A_(row_offsets_.size() - 1, nullptr), + x_(col_offsets_), y_(row_offsets_) +{ + +} + +BlockDiagOperator::BlockDiagOperator(std::vector row_offsets, std::vector col_offsets) + : Operator(row_offsets.back(), col_offsets.back()), + row_offsets_(std::move(row_offsets)), col_offsets_(std::move(col_offsets)), + A_(row_offsets_.size() - 1, nullptr), + x_(col_offsets_), y_(row_offsets_) +{ + +} + +BlockDiagOperator::BlockDiagOperator(const BlockDiagOperator& other) noexcept + : Operator(other), + row_offsets_(other.row_offsets_), col_offsets_(other.col_offsets_), + A_(other.A_), x_(other.x_), y_(other.y_) +{ + +} + +BlockDiagOperator::BlockDiagOperator(BlockDiagOperator&& other) noexcept +{ + swap(*this, other); +} + +BlockDiagOperator& BlockDiagOperator::operator=(BlockDiagOperator&& other) noexcept +{ + swap(*this, other); + + return *this; +} + +void swap(BlockDiagOperator& lhs, BlockDiagOperator& rhs) noexcept +{ + swap(static_cast(lhs), static_cast(rhs)); + + swap(lhs.row_offsets_, rhs.row_offsets_); + swap(lhs.col_offsets_, rhs.col_offsets_); + swap(lhs.A_, rhs.A_); + swap(lhs.x_, rhs.x_); + swap(lhs.y_, rhs.y_); +} + + +const std::vector& BlockDiagOperator::GetRowOffsets() const +{ + return row_offsets_; +} + +const std::vector& BlockDiagOperator::GetColOffsets() const +{ + return col_offsets_; +} + +const Operator* BlockDiagOperator::GetBlock(int i, int j) const +{ + assert(i < static_cast(row_offsets_.size()) - 1); + assert(j < static_cast(col_offsets_.size()) - 1); + assert(i == j); + + return A_[i]; +} + +void BlockDiagOperator::SetBlock(int i, int j, const Operator& op) +{ + assert(i < static_cast(row_offsets_.size()) - 1); + assert(j < static_cast(col_offsets_.size()) - 1); + + assert(op.Rows() == (row_offsets_[i + 1] - row_offsets_[i])); + assert(op.Cols() == (col_offsets_[j + 1] - col_offsets_[j])); + + assert(i == j); + + A_[i] = &op; +} + +void BlockDiagOperator::SetBlock(int i, const Operator& op) +{ + SetBlock(i, i, op); +} + +void BlockDiagOperator::Mult(const VectorView& input, VectorView output) const +{ + assert(input.size() == cols_); + assert(output.size() == rows_); + + x_ = input; + y_ = 0.0; + + const int num_blocks = row_offsets_.size() - 1; + + for (int i = 0; i < num_blocks; ++i) + { + const Operator* op = A_[i]; + + if (op) + { + op->Mult(x_.GetBlock(i), y_.GetBlock(i)); + } + } + + output = y_; +} + +void BlockDiagOperator::MultAT(const VectorView& input, VectorView output) const +{ + assert(input.size() == rows_); + assert(output.size() == cols_); + + y_ = input; + x_ = 0.0; + + const int num_blocks = row_offsets_.size() - 1; + + for (int i = 0; i < num_blocks; ++i) + { + const Operator* op = A_[i]; + + if (op) + { + op->MultAT(y_.GetBlock(i), x_.GetBlock(i)); + } + } + + output = x_; +} + } // namespace linalgcpp diff --git a/src/solvers.cpp b/src/solvers.cpp index ba2fe0b..140a09a 100644 --- a/src/solvers.cpp +++ b/src/solvers.cpp @@ -23,6 +23,26 @@ Solver::Solver(const Operator& A, int max_iter, double rel_tol, assert(A_->Rows() == A_->Cols()); } +Solver::Solver(const Solver& other) noexcept + : Operator(other), A_(other.A_), max_iter_(other.max_iter_), + verbose_(other.verbose_), rel_tol_(other.rel_tol_), + abs_tol_(other.abs_tol_), Dot_(other.Dot_), + num_iter_(other.num_iter_) +{ +} + +Solver::Solver(Solver&& other) noexcept +{ + swap(*this, other); +} + +Solver& Solver::operator=(Solver&& other) noexcept +{ + swap(*this, other); + + return *this; +} + void swap(Solver& lhs, Solver& rhs) noexcept { swap(static_cast(lhs), static_cast(rhs)); @@ -185,11 +205,11 @@ void PCGSolver::Mult(const VectorView& b, VectorView x) const double r_z = (*Dot_)(r_ , z_); double alpha = r_z / pAp; - linalgcpp_verify(pAp > -1e12, "PCG is not positive definite!"); + linalgcpp_verify(pAp > -1e-12, "PCG is not positive definite!"); x.Add(alpha, p_); - double denom = (*Dot_)(z_, r_); + double denom = r_z; r_.Sub(alpha, Ap_); @@ -202,19 +222,20 @@ void PCGSolver::Mult(const VectorView& b, VectorView x) const z_ = r_; } - double numer = (*Dot_)(z_, r_); + double r_z_next = (*Dot_)(z_, r_); if (verbose_) { - printf("PCG %d: %.2e\n", num_iter_, numer / r0); + printf("PCG %d: reduction: %.2e r_z_next: %.2e, tol2: %.2e\n", + num_iter_, r_z_next / r0, r_z_next, tol_tol); } - if (numer < tol_tol) + if (r_z_next < tol_tol) { break; } - double beta = numer / denom; + double beta = r_z_next / denom; p_ *= beta; p_ += z_; diff --git a/src/utilities.cpp b/src/utilities.cpp new file mode 100644 index 0000000..4aa4232 --- /dev/null +++ b/src/utilities.cpp @@ -0,0 +1,7 @@ +#include "utilities.hpp" + +namespace linalgcpp +{ + + +} // namespace linalgcpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 35bc427..738f114 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,10 +1,4 @@ add_executable(test_linalg test_linalg.cpp) target_link_libraries(test_linalg linalgcpp) -set_target_properties(test_linalg PROPERTIES - CXX_STANDARD 11 - CXX_STANDARD_REQUIRED YES - CXX_EXTENSIONS NO -) - add_test(test_linalg test_linalg) diff --git a/tests/test_linalg.cpp b/tests/test_linalg.cpp index 839cc30..f3c55e8 100644 --- a/tests/test_linalg.cpp +++ b/tests/test_linalg.cpp @@ -1243,8 +1243,10 @@ void test_blockvector() void test_blockoperator() { CooMatrix coo(2, 2); - coo.Add(0, 0, 1.0); - coo.Add(1, 1, 1.0); + coo.Add(0, 0, 3.0); + coo.Add(0, 1, -1.0); + coo.Add(1, 0, -1.0); + coo.Add(1, 1, 4.0); auto sparse = coo.ToSparse(); auto dense = coo.ToDense(); @@ -1296,6 +1298,19 @@ void test_blockoperator() y.Print("y T:"); } + // Diag Block + { + BlockDiagOperator b(offsets, offsets); + b.SetBlock(0, 0, coo); + b.SetBlock(1, 1, dense); + + b.Mult(x, y); + y.Print("y:"); + + b.MultAT(x, y); + y.Print("y T:"); + } + } void test_timer() @@ -1521,6 +1536,81 @@ void test_argparser(int argc, char** argv) } } +void test_assert() +{ + bool verify_works = true; + bool assert_works = true; + + try + { + linalgcpp_verify(false, "Catch this!"); + verify_works = false; + } + catch (const std::runtime_error& e) + { + std::cout << "Verification caught: " << e.what() << "\n"; + } + + try + { + bool val = false; + linalgcpp_verify([&val]() { return val; } , "Catch this lambda!"); + verify_works = false; + } + catch (const std::runtime_error& e) + { + std::cout << "Verification caught: " << e.what() << "\n"; + } + + try + { + linalgcpp_verify(true, "Don't catch this!"); + } + catch (const std::runtime_error& e) + { + std::cout << "Verification caught: " << e.what() << "\n"; + verify_works = false; + } + + try + { + linalgcpp_assert(false, "Catch this!"); + assert_works = false; + } + catch (const std::runtime_error& e) + { + std::cout << "Assertion caught: " << e.what() << "\n"; + } + + bool eval = false; + try + { + bool val = false; + // In Debug mode, lambda is evaluated + // In Release mode, lambda is NOT evaluated + linalgcpp_assert([&]() { eval = true; return val; } , "Catch this lambda!"); + assert_works = false; + } + catch (const std::runtime_error& e) + { + std::cout << "Assertion caught: " << e.what() << "\n"; + } + std::cout << "Lambda Evaluated: " << std::boolalpha << eval << "\n"; + + try + { + linalgcpp_assert(true, "Don't catch this!"); + } + catch (const std::runtime_error& e) + { + assert_works = false; + std::cout << "Assertion caught: " << e.what() << "\n"; + } + + std::cout << "Verfication Works: " << std::boolalpha << verify_works << "\n"; + std::cout << "Assertation Works: " << std::boolalpha << assert_works << "\n"; +} + int main(int argc, char** argv) { test_coo(); @@ -1537,6 +1627,7 @@ int main(int argc, char** argv) test_timer(); test_eigensolve(); test_argparser(argc, argv); + test_assert(); return EXIT_SUCCESS; }