From a987300bd514595f827d514a50cc12c6c45c9e66 Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Wed, 31 Oct 2018 11:33:59 -0700 Subject: [PATCH 01/10] Clean up CMake files, add examples to tests --- examples/dense/CMakeLists.txt | 4 ++++ examples/dense/dense_lu.cpp | 2 +- examples/dense/dense_qr.cpp | 4 ++-- examples/sparse/CMakeLists.txt | 2 ++ examples/vector/CMakeLists.txt | 2 ++ modules/parlinalgcpp/CMakeLists.txt | 6 ------ modules/parlinalgcpp/tests/CMakeLists.txt | 6 ------ modules/partition/tests/CMakeLists.txt | 7 ------- modules/sparsesolver/CMakeLists.txt | 9 --------- modules/sparsesolver/tests/CMakeLists.txt | 8 -------- tests/CMakeLists.txt | 6 ------ 11 files changed, 11 insertions(+), 45 deletions(-) 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/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/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/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/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/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/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) From d1d79a9842fd7739f8837e807472c4c7bcc47e6d Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Fri, 2 Nov 2018 20:20:51 -0700 Subject: [PATCH 02/10] Initial setup for block diag comp smoother --- include/blockoperator.hpp | 3 + include/utilities.hpp | 63 +++++++++ modules/parlinalgcpp/CMakeLists.txt | 1 + modules/parlinalgcpp/include/parlinalgcpp.hpp | 1 + modules/parlinalgcpp/include/parmatrix.hpp | 27 ++++ modules/parlinalgcpp/include/parprecond.hpp | 40 ++++++ modules/parlinalgcpp/src/parprecond.cpp | 120 ++++++++++++++++++ modules/parlinalgcpp/tests/test_parlinalg.cpp | 38 +++++- src/blockoperator.cpp | 7 + src/utilities.cpp | 7 + 10 files changed, 303 insertions(+), 4 deletions(-) create mode 100644 modules/parlinalgcpp/include/parprecond.hpp create mode 100644 modules/parlinalgcpp/src/parprecond.cpp create mode 100644 src/utilities.cpp diff --git a/include/blockoperator.hpp b/include/blockoperator.hpp index c5fbf1a..e8f72d4 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; diff --git a/include/utilities.hpp b/include/utilities.hpp index cde6cb1..9418eb0 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 { @@ -28,6 +41,56 @@ void linalgcpp_verify(bool expression, const std::string& message = "linalgcpp v } } +/** @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/modules/parlinalgcpp/CMakeLists.txt b/modules/parlinalgcpp/CMakeLists.txt index 1037790..e752f50 100644 --- a/modules/parlinalgcpp/CMakeLists.txt +++ b/modules/parlinalgcpp/CMakeLists.txt @@ -38,6 +38,7 @@ add_library(parlinalgcpp src/parcommpkg.cpp src/parmatrix.cpp src/paroperator.cpp + src/parprecond.cpp src/parsmoother.cpp src/parsolvers.cpp src/parutilities.cpp diff --git a/modules/parlinalgcpp/include/parlinalgcpp.hpp b/modules/parlinalgcpp/include/parlinalgcpp.hpp index 59feb1a..3e27cdf 100644 --- a/modules/parlinalgcpp/include/parlinalgcpp.hpp +++ b/modules/parlinalgcpp/include/parlinalgcpp.hpp @@ -5,5 +5,6 @@ #include "parsolvers.hpp" #include "parutilities.hpp" #include "paroperator.hpp" +#include "parprecond.hpp" #include "parsmoother.hpp" #include "parcommpkg.hpp" diff --git a/modules/parlinalgcpp/include/parmatrix.hpp b/modules/parlinalgcpp/include/parmatrix.hpp index 193aa7f..12573b7 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 @@ -307,6 +316,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/parprecond.hpp b/modules/parlinalgcpp/include/parprecond.hpp new file mode 100644 index 0000000..cff78ee --- /dev/null +++ b/modules/parlinalgcpp/include/parprecond.hpp @@ -0,0 +1,40 @@ +/*! @file */ + +#ifndef PARPRECOND_HPP__ +#define PARPRECOND_HPP__ + +#include "parmatrix.hpp" + +namespace linalgcpp +{ + +class ParBlockDiagComp : public ParOperator +{ + public: + ParBlockDiagComp() = default; + ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vertex); + + ~ParBlockDiagComp() = default; + + /*! @brief Apply the action of the solver + @param input input vector x + @param output output vector y + */ + void Mult(const linalgcpp::VectorView& input, + linalgcpp::VectorView output) const; + + private: + ParMatrix MakeRedistributor(const ParMatrix& agg_vertex); + + ParMatrix redist_; + + BlockOperator block_op_; + std::vector> solvers_; + + mutable Vector sol_r_; + mutable Vector rhs_r_; +}; + +} //namespace linalgcpp + +#endif // PARPRECOND_HPP diff --git a/modules/parlinalgcpp/src/parprecond.cpp b/modules/parlinalgcpp/src/parprecond.cpp new file mode 100644 index 0000000..6301319 --- /dev/null +++ b/modules/parlinalgcpp/src/parprecond.cpp @@ -0,0 +1,120 @@ +/*! @file */ + +#include "parprecond.hpp" + +namespace linalgcpp +{ + +ParBlockDiagComp::ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vertex) + : ParOperator(A), redist_(MakeRedistributor(agg_vertex)), + sol_r_(redist_.Rows()), rhs_r_(redist_.Rows()) +{ + ParMatrix redist_T = redist_.Transpose(); + ParMatrix agg_vertex_r = agg_vertex.Mult(redist_T); + ParMatrix vertex_agg_r = agg_vertex_r.Transpose(); + ParMatrix A_r = redist_.Mult(A).Mult(redist_T); + + linalgcpp_verify(agg_vertex_r.GetOffd().nnz() == 0); + + block_op_ = BlockOperator(agg_vertex_r.GetDiag().GetIndptr()); + solvers_.resize(agg_vertex_r.Rows()); + + std::vector marker(redist_.Rows(), -1); + int num_aggs = agg_vertex.Rows(); + int num_vertices = agg_vertex.Cols(); + + // Diagonal compensation + std::vector A_diag(A_r.Rows()); + { + auto& diag_indptr = A_r.GetDiag().GetIndptr(); + auto& diag_indices = A_r.GetDiag().GetIndices(); + auto& diag_data = A_r.GetDiag().GetData(); + auto& offd_indptr = A_r.GetOffd().GetIndptr(); + auto& offd_indices = A_r.GetOffd().GetIndices(); + auto& offd_data = A_r.GetOffd().GetData(); + auto& col_map = A_r.GetColMap(); + auto& part = vertex_agg_r.GetDiag().GetIndices(); + + for (int row = 0; row < num_vertices; ++row) + { + for (int j = diag_indptr[row]; j < diag_indptr[row + 1]; ++j) + { + int col = diag_indices[j]; + + if (col > row) + { + int agg_i = part[row]; + int agg_j = part[col]; + + if (agg_i != agg_j) + { + A_diag[row] += std::fabs(diag_data[j]); + A_diag[col] += std::fabs(diag_data[j]); + + diag_data[j] = 0.0; // Not necessary + } + } + } + + for (int j = offd_indptr[row]; j < offd_indptr[row + 1]; ++j) + { + int agg_i = part[row]; + A_diag[row] += std::fabs(offd_data[j]); + offd_data[j] = 0.0; // Not necessary + } + } + } + + A_r.AddDiag(A_diag); + + if (A_r.GetMyId() == 0) + { + A_r.GetDiag().ToDense().Print("A diag:"); + A_r.GetOffd().ToDense().Print("A offd:"); + } + + for (int i = 0; i < num_aggs; ++i) + { + auto indices = agg_vertex_r.GetDiag().GetIndices(i); + auto A_agg = A_r.GetDiag().GetSubMatrix(indices, indices, marker); + + DenseMatrix A_dense(A_agg.Rows(), A_agg.Cols()); + A_agg.ToDense(A_dense); + + A_dense.Print("A sub:"); + A_dense.Invert(); + A_dense.Print("A sub inv:"); + + solvers_[i] = make_unique(std::move(A_dense)); + block_op_.SetBlock(i, i, *solvers_[i]); + } +} + +ParMatrix ParBlockDiagComp::MakeRedistributor(const ParMatrix& agg_vertex) +{ + int diag_rows = agg_vertex.GetDiag().nnz(); + int diag_cols = agg_vertex.GetDiag().Cols(); + + int offd_rows = agg_vertex.GetOffd().nnz(); + int offd_cols = agg_vertex.GetOffd().Cols(); + + int total_rows = diag_rows + offd_rows; + + SparseMatrix diag = SparseIdentity(total_rows, diag_cols); + SparseMatrix offd = SparseIdentity(total_rows, offd_cols, diag_rows); + + auto row_starts = GenerateOffsets(agg_vertex.GetComm(), total_rows); + + return ParMatrix(agg_vertex.GetComm(), row_starts, agg_vertex.GetColStarts(), + std::move(diag), std::move(offd), agg_vertex.GetColMap()); +} + +void ParBlockDiagComp::Mult(const linalgcpp::VectorView& input, + linalgcpp::VectorView output) const +{ + redist_.Mult(input, rhs_r_); + block_op_.Mult(rhs_r_, sol_r_); + redist_.MultAT(sol_r_, output); +} + +} // namespace linalgcpp diff --git a/modules/parlinalgcpp/tests/test_parlinalg.cpp b/modules/parlinalgcpp/tests/test_parlinalg.cpp index 5773e16..1a37b0d 100644 --- a/modules/parlinalgcpp/tests/test_parlinalg.cpp +++ b/modules/parlinalgcpp/tests/test_parlinalg.cpp @@ -8,6 +8,7 @@ using namespace linalgcpp; void PrintVector(MPI_Comm comm, const std::vector>& vects, const std::vector& evals); +SparseMatrix make_agg_vertex(const ParMatrix& A); int main(int argc, char** argv) { @@ -88,10 +89,14 @@ int main(int argc, char** argv) ParCG pcg_direct(pmat, smoother_direct, cg_max_iter); ParCG pcg_direct_blas(pmat, smoother_direct_blas, cg_max_iter); + SparseMatrix agg_vertex = make_agg_vertex(pmat); + ParMatrix par_agg_vertex(pmat.GetComm(), std::move(agg_vertex)); + ParBlockDiagComp block_comp(pmat, par_agg_vertex); + test_pcg_copy = cg; // rofl - std::vector ops { + std::vector ops { &cg, &pcg_boomer, &pcg_parasails, @@ -126,6 +131,7 @@ int main(int argc, char** argv) &pcg_direct_blas, &test_smooth_copy, &test_pcg_copy, + &block_comp, }; //ops = {&smoother_direct_blas}; @@ -136,7 +142,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 +173,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 +193,22 @@ void PrintVector(MPI_Comm comm, const std::vector>& ev } } + +SparseMatrix make_agg_vertex(const ParMatrix& A) +{ + int num_aggs = 1; + 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; + + 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/src/blockoperator.cpp b/src/blockoperator.cpp index a7cc0a4..9de24e9 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)); 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 From d99703f0616dfc6b2dd0c2c6891ffe0c847333c2 Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Sun, 4 Nov 2018 13:00:38 -0800 Subject: [PATCH 03/10] Add BlockDiagOperator and utilities --- include/blockoperator.hpp | 104 +++++++++++++- include/linalgcpp.hpp | 1 + include/solvers.hpp | 7 +- include/vectorview.hpp | 48 +++++++ modules/parlinalgcpp/include/parprecond.hpp | 2 +- modules/parlinalgcpp/include/parutilities.hpp | 29 ++++ modules/parlinalgcpp/src/parprecond.cpp | 9 +- modules/parlinalgcpp/tests/test_parlinalg.cpp | 9 +- src/blockoperator.cpp | 134 ++++++++++++++++++ src/solvers.cpp | 25 +++- tests/test_linalg.cpp | 19 ++- 11 files changed, 373 insertions(+), 14 deletions(-) diff --git a/include/blockoperator.hpp b/include/blockoperator.hpp index e8f72d4..e837649 100644 --- a/include/blockoperator.hpp +++ b/include/blockoperator.hpp @@ -65,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); @@ -100,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/vectorview.hpp b/include/vectorview.hpp index 0a7e721..033f586 100644 --- a/include/vectorview.hpp +++ b/include/vectorview.hpp @@ -87,6 +87,26 @@ class VectorView */ virtual const T* end() const; + /*! @brief STL like back. Points to the last element + @retval pointer to the last element + */ + virtual T& back(); + + /*! @brief STL like back. Points to the last element + @retval pointer to the last element + */ + virtual const T& back() const; + + /*! @brief STL like front. Points to the first element + @retval pointer to the first element + */ + virtual T& front(); + + /*! @brief STL like front. Points to the first element + @retval pointer to the first element + */ + virtual const T& front() 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::back() +{ + assert(size_ >= 1); + return data_[size_ - 1]; +} + +template +const T& VectorView::back() const +{ + assert(size_ >= 1); + return data_[size_ - 1]; +} + +template +T& VectorView::front() +{ + assert(size_ >= 1); + return data_[0]; +} + +template +const T& VectorView::front() const +{ + assert(size_ >= 1); + return data_[0]; +} + template T& VectorView::operator[](int i) { diff --git a/modules/parlinalgcpp/include/parprecond.hpp b/modules/parlinalgcpp/include/parprecond.hpp index cff78ee..747e755 100644 --- a/modules/parlinalgcpp/include/parprecond.hpp +++ b/modules/parlinalgcpp/include/parprecond.hpp @@ -28,7 +28,7 @@ class ParBlockDiagComp : public ParOperator ParMatrix redist_; - BlockOperator block_op_; + BlockDiagOperator block_op_; std::vector> solvers_; mutable Vector sol_r_; diff --git a/modules/parlinalgcpp/include/parutilities.hpp b/modules/parlinalgcpp/include/parutilities.hpp index 3f02740..fbe69dc 100644 --- a/modules/parlinalgcpp/include/parutilities.hpp +++ b/modules/parlinalgcpp/include/parutilities.hpp @@ -12,6 +12,10 @@ namespace linalgcpp { +/// Call output only on processor 0 +#define ParPrint(myid, output) if (myid == 0) output + + class ParMatrix; /// Split a square matrix between processes @@ -69,5 +73,30 @@ 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); +/** @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/parprecond.cpp b/modules/parlinalgcpp/src/parprecond.cpp index 6301319..13fe275 100644 --- a/modules/parlinalgcpp/src/parprecond.cpp +++ b/modules/parlinalgcpp/src/parprecond.cpp @@ -16,7 +16,7 @@ ParBlockDiagComp::ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vert linalgcpp_verify(agg_vertex_r.GetOffd().nnz() == 0); - block_op_ = BlockOperator(agg_vertex_r.GetDiag().GetIndptr()); + block_op_ = BlockDiagOperator(agg_vertex_r.GetDiag().GetIndptr()); solvers_.resize(agg_vertex_r.Rows()); std::vector marker(redist_.Rows(), -1); @@ -35,6 +35,8 @@ ParBlockDiagComp::ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vert auto& col_map = A_r.GetColMap(); auto& part = vertex_agg_r.GetDiag().GetIndices(); + const auto& A_ii = A_r.GetDiag().GetDiag(); + for (int row = 0; row < num_vertices; ++row) { for (int j = diag_indptr[row]; j < diag_indptr[row + 1]; ++j) @@ -48,8 +50,9 @@ ParBlockDiagComp::ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vert if (agg_i != agg_j) { - A_diag[row] += std::fabs(diag_data[j]); - A_diag[col] += std::fabs(diag_data[j]); + double tau = std::sqrt(A_ii[row] / A_ii[col]); + A_diag[row] += std::fabs(diag_data[j]) * tau; + A_diag[col] += std::fabs(diag_data[j]) / tau; diag_data[j] = 0.0; // Not necessary } diff --git a/modules/parlinalgcpp/tests/test_parlinalg.cpp b/modules/parlinalgcpp/tests/test_parlinalg.cpp index 1a37b0d..9a0ab1a 100644 --- a/modules/parlinalgcpp/tests/test_parlinalg.cpp +++ b/modules/parlinalgcpp/tests/test_parlinalg.cpp @@ -196,15 +196,20 @@ void PrintVector(MPI_Comm comm, const std::vector>& ev SparseMatrix make_agg_vertex(const ParMatrix& A) { - int num_aggs = 1; + 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; + //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); diff --git a/src/blockoperator.cpp b/src/blockoperator.cpp index 9de24e9..96d80a0 100644 --- a/src/blockoperator.cpp +++ b/src/blockoperator.cpp @@ -156,4 +156,138 @@ 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::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..21d1a31 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,7 +205,7 @@ 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-15, "PCG is not positive definite!"); x.Add(alpha, p_); @@ -206,7 +226,8 @@ void PCGSolver::Mult(const VectorView& b, VectorView x) const if (verbose_) { - printf("PCG %d: %.2e\n", num_iter_, numer / r0); + printf("PCG %d: reduction: %.2e numer: %.2e, tol2: %.2e\n", + num_iter_, numer / r0, numer, tol_tol); } if (numer < tol_tol) diff --git a/tests/test_linalg.cpp b/tests/test_linalg.cpp index 839cc30..7e2eadd 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() From 3e446d9ac370f0a4e1698cbae7a73dcacd5c83db Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Sun, 4 Nov 2018 13:13:24 -0800 Subject: [PATCH 04/10] Fix PCGSolver, add utilities --- include/vectorview.hpp | 48 +++++++++++++++++++ modules/parlinalgcpp/include/parutilities.hpp | 28 +++++++++++ src/solvers.cpp | 2 +- 3 files changed, 77 insertions(+), 1 deletion(-) 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/include/parutilities.hpp b/modules/parlinalgcpp/include/parutilities.hpp index 3f02740..72791fd 100644 --- a/modules/parlinalgcpp/include/parutilities.hpp +++ b/modules/parlinalgcpp/include/parutilities.hpp @@ -12,6 +12,9 @@ namespace linalgcpp { +/// Call output only on processor 0 +#define ParPrint(myid, output) if (myid == 0) { output; } + class ParMatrix; /// Split a square matrix between processes @@ -69,5 +72,30 @@ 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); +/** @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/src/solvers.cpp b/src/solvers.cpp index ba2fe0b..1f617c3 100644 --- a/src/solvers.cpp +++ b/src/solvers.cpp @@ -185,7 +185,7 @@ 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_); From 886db9a1af740fe04e519add7971a4c186f0651c Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Sun, 4 Nov 2018 13:50:48 -0800 Subject: [PATCH 05/10] Add assert/verify tests --- examples/include/ex_utilities.hpp | 9 ++---- include/linalgcpp.hpp | 1 + include/utilities.hpp | 2 +- tests/test_linalg.cpp | 50 +++++++++++++++++++++++++++++++ 4 files changed, 55 insertions(+), 7 deletions(-) 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/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/utilities.hpp b/include/utilities.hpp index cde6cb1..90a016d 100644 --- a/include/utilities.hpp +++ b/include/utilities.hpp @@ -10,7 +10,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); diff --git a/tests/test_linalg.cpp b/tests/test_linalg.cpp index 839cc30..25cdfba 100644 --- a/tests/test_linalg.cpp +++ b/tests/test_linalg.cpp @@ -1521,6 +1521,55 @@ 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 + { + 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"; + } + + 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 +1586,7 @@ int main(int argc, char** argv) test_timer(); test_eigensolve(); test_argparser(argc, argv); + test_assert(); return EXIT_SUCCESS; } From 793b37ae0e046a34b9d82319aff9a61d5da642d0 Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Wed, 7 Nov 2018 10:35:14 -0800 Subject: [PATCH 06/10] Add utilites to ParMatrix --- modules/parlinalgcpp/include/parmatrix.hpp | 39 ++++++++++++++++++++++ modules/parlinalgcpp/src/parmatrix.cpp | 13 ++++++++ 2 files changed, 52 insertions(+) diff --git a/modules/parlinalgcpp/include/parmatrix.hpp b/modules/parlinalgcpp/include/parmatrix.hpp index 193aa7f..bc89040 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,17 @@ 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 + */ + void EliminateZeros(double tol = 0.0); + private: void Init(); void HypreCreate(); @@ -307,6 +327,25 @@ 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/src/parmatrix.cpp b/modules/parlinalgcpp/src/parmatrix.cpp index 9ec0582..f02aaf4 100644 --- a/modules/parlinalgcpp/src/parmatrix.cpp +++ b/modules/parlinalgcpp/src/parmatrix.cpp @@ -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) +{ + diag_.EliminateZeros(tol); + offd_.EliminateZeros(tol); + + HypreCreate(); +} + } // namespace linalgcpp From b89c03f7f88e3aa9e97ae096213bacfcad3f5c80 Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Mon, 10 Dec 2018 10:21:45 -0800 Subject: [PATCH 07/10] Add parallel utilities, remove parprecond for now --- include/sparsematrix.hpp | 10 +- include/utilities.hpp | 22 ++++ modules/parlinalgcpp/CMakeLists.txt | 1 - modules/parlinalgcpp/include/parcommpkg.hpp | 1 + modules/parlinalgcpp/include/parmatrix.hpp | 3 +- modules/parlinalgcpp/include/paroperator.hpp | 6 + modules/parlinalgcpp/include/parprecond.hpp | 40 ------ modules/parlinalgcpp/include/parutilities.hpp | 55 ++++++++ modules/parlinalgcpp/src/parcommpkg.cpp | 4 +- modules/parlinalgcpp/src/parmatrix.cpp | 6 +- modules/parlinalgcpp/src/paroperator.cpp | 9 ++ modules/parlinalgcpp/src/parprecond.cpp | 123 ------------------ modules/parlinalgcpp/src/parutilities.cpp | 32 ++++- modules/parlinalgcpp/tests/test_parlinalg.cpp | 6 - modules/partition/include/partition.hpp | 9 +- modules/sparsesolver/include/sparsesolve.hpp | 5 + modules/sparsesolver/src/sparsesolve.cpp | 4 + src/blockoperator.cpp | 7 +- tests/test_linalg.cpp | 26 ++++ 19 files changed, 183 insertions(+), 186 deletions(-) delete mode 100644 modules/parlinalgcpp/include/parprecond.hpp delete mode 100644 modules/parlinalgcpp/src/parprecond.cpp 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 630ec58..8cccebb 100644 --- a/include/utilities.hpp +++ b/include/utilities.hpp @@ -41,6 +41,28 @@ 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 diff --git a/modules/parlinalgcpp/CMakeLists.txt b/modules/parlinalgcpp/CMakeLists.txt index e752f50..1037790 100644 --- a/modules/parlinalgcpp/CMakeLists.txt +++ b/modules/parlinalgcpp/CMakeLists.txt @@ -38,7 +38,6 @@ add_library(parlinalgcpp src/parcommpkg.cpp src/parmatrix.cpp src/paroperator.cpp - src/parprecond.cpp src/parsmoother.cpp src/parsolvers.cpp src/parutilities.cpp 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 55ca90f..4dedb46 100644 --- a/modules/parlinalgcpp/include/parmatrix.hpp +++ b/modules/parlinalgcpp/include/parmatrix.hpp @@ -293,8 +293,9 @@ class ParMatrix: public ParOperator /*! @brief Remove entries less than tolerance @param tol tolerance + @param keep_diag keep diagonal element, regardless of value */ - void EliminateZeros(double tol = 0.0); + void EliminateZeros(double tol = 0.0, bool keep_diag = false); private: void Init(); 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/parprecond.hpp b/modules/parlinalgcpp/include/parprecond.hpp deleted file mode 100644 index 747e755..0000000 --- a/modules/parlinalgcpp/include/parprecond.hpp +++ /dev/null @@ -1,40 +0,0 @@ -/*! @file */ - -#ifndef PARPRECOND_HPP__ -#define PARPRECOND_HPP__ - -#include "parmatrix.hpp" - -namespace linalgcpp -{ - -class ParBlockDiagComp : public ParOperator -{ - public: - ParBlockDiagComp() = default; - ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vertex); - - ~ParBlockDiagComp() = default; - - /*! @brief Apply the action of the solver - @param input input vector x - @param output output vector y - */ - void Mult(const linalgcpp::VectorView& input, - linalgcpp::VectorView output) const; - - private: - ParMatrix MakeRedistributor(const ParMatrix& agg_vertex); - - ParMatrix redist_; - - BlockDiagOperator block_op_; - std::vector> solvers_; - - mutable Vector sol_r_; - mutable Vector rhs_r_; -}; - -} //namespace linalgcpp - -#endif // PARPRECOND_HPP diff --git a/modules/parlinalgcpp/include/parutilities.hpp b/modules/parlinalgcpp/include/parutilities.hpp index 72791fd..a6b39a0 100644 --- a/modules/parlinalgcpp/include/parutilities.hpp +++ b/modules/parlinalgcpp/include/parutilities.hpp @@ -15,11 +15,60 @@ 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); @@ -72,6 +121,12 @@ 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 { 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 f02aaf4..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_); } @@ -402,9 +402,9 @@ int ParMatrix::RowSize(int row) return diag_.RowSize(row) + offd_.RowSize(row); } -void ParMatrix::EliminateZeros(double tol) +void ParMatrix::EliminateZeros(double tol, bool keep_diag) { - diag_.EliminateZeros(tol); + diag_.EliminateZeros(tol, keep_diag); offd_.EliminateZeros(tol); HypreCreate(); 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/parprecond.cpp b/modules/parlinalgcpp/src/parprecond.cpp deleted file mode 100644 index 13fe275..0000000 --- a/modules/parlinalgcpp/src/parprecond.cpp +++ /dev/null @@ -1,123 +0,0 @@ -/*! @file */ - -#include "parprecond.hpp" - -namespace linalgcpp -{ - -ParBlockDiagComp::ParBlockDiagComp(const ParMatrix& A, const ParMatrix& agg_vertex) - : ParOperator(A), redist_(MakeRedistributor(agg_vertex)), - sol_r_(redist_.Rows()), rhs_r_(redist_.Rows()) -{ - ParMatrix redist_T = redist_.Transpose(); - ParMatrix agg_vertex_r = agg_vertex.Mult(redist_T); - ParMatrix vertex_agg_r = agg_vertex_r.Transpose(); - ParMatrix A_r = redist_.Mult(A).Mult(redist_T); - - linalgcpp_verify(agg_vertex_r.GetOffd().nnz() == 0); - - block_op_ = BlockDiagOperator(agg_vertex_r.GetDiag().GetIndptr()); - solvers_.resize(agg_vertex_r.Rows()); - - std::vector marker(redist_.Rows(), -1); - int num_aggs = agg_vertex.Rows(); - int num_vertices = agg_vertex.Cols(); - - // Diagonal compensation - std::vector A_diag(A_r.Rows()); - { - auto& diag_indptr = A_r.GetDiag().GetIndptr(); - auto& diag_indices = A_r.GetDiag().GetIndices(); - auto& diag_data = A_r.GetDiag().GetData(); - auto& offd_indptr = A_r.GetOffd().GetIndptr(); - auto& offd_indices = A_r.GetOffd().GetIndices(); - auto& offd_data = A_r.GetOffd().GetData(); - auto& col_map = A_r.GetColMap(); - auto& part = vertex_agg_r.GetDiag().GetIndices(); - - const auto& A_ii = A_r.GetDiag().GetDiag(); - - for (int row = 0; row < num_vertices; ++row) - { - for (int j = diag_indptr[row]; j < diag_indptr[row + 1]; ++j) - { - int col = diag_indices[j]; - - if (col > row) - { - int agg_i = part[row]; - int agg_j = part[col]; - - if (agg_i != agg_j) - { - double tau = std::sqrt(A_ii[row] / A_ii[col]); - A_diag[row] += std::fabs(diag_data[j]) * tau; - A_diag[col] += std::fabs(diag_data[j]) / tau; - - diag_data[j] = 0.0; // Not necessary - } - } - } - - for (int j = offd_indptr[row]; j < offd_indptr[row + 1]; ++j) - { - int agg_i = part[row]; - A_diag[row] += std::fabs(offd_data[j]); - offd_data[j] = 0.0; // Not necessary - } - } - } - - A_r.AddDiag(A_diag); - - if (A_r.GetMyId() == 0) - { - A_r.GetDiag().ToDense().Print("A diag:"); - A_r.GetOffd().ToDense().Print("A offd:"); - } - - for (int i = 0; i < num_aggs; ++i) - { - auto indices = agg_vertex_r.GetDiag().GetIndices(i); - auto A_agg = A_r.GetDiag().GetSubMatrix(indices, indices, marker); - - DenseMatrix A_dense(A_agg.Rows(), A_agg.Cols()); - A_agg.ToDense(A_dense); - - A_dense.Print("A sub:"); - A_dense.Invert(); - A_dense.Print("A sub inv:"); - - solvers_[i] = make_unique(std::move(A_dense)); - block_op_.SetBlock(i, i, *solvers_[i]); - } -} - -ParMatrix ParBlockDiagComp::MakeRedistributor(const ParMatrix& agg_vertex) -{ - int diag_rows = agg_vertex.GetDiag().nnz(); - int diag_cols = agg_vertex.GetDiag().Cols(); - - int offd_rows = agg_vertex.GetOffd().nnz(); - int offd_cols = agg_vertex.GetOffd().Cols(); - - int total_rows = diag_rows + offd_rows; - - SparseMatrix diag = SparseIdentity(total_rows, diag_cols); - SparseMatrix offd = SparseIdentity(total_rows, offd_cols, diag_rows); - - auto row_starts = GenerateOffsets(agg_vertex.GetComm(), total_rows); - - return ParMatrix(agg_vertex.GetComm(), row_starts, agg_vertex.GetColStarts(), - std::move(diag), std::move(offd), agg_vertex.GetColMap()); -} - -void ParBlockDiagComp::Mult(const linalgcpp::VectorView& input, - linalgcpp::VectorView output) const -{ - redist_.Mult(input, rhs_r_); - block_op_.Mult(rhs_r_, sol_r_); - redist_.MultAT(sol_r_, output); -} - -} // namespace linalgcpp 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/test_parlinalg.cpp b/modules/parlinalgcpp/tests/test_parlinalg.cpp index 9a0ab1a..4f185c4 100644 --- a/modules/parlinalgcpp/tests/test_parlinalg.cpp +++ b/modules/parlinalgcpp/tests/test_parlinalg.cpp @@ -8,7 +8,6 @@ using namespace linalgcpp; void PrintVector(MPI_Comm comm, const std::vector>& vects, const std::vector& evals); -SparseMatrix make_agg_vertex(const ParMatrix& A); int main(int argc, char** argv) { @@ -89,10 +88,6 @@ int main(int argc, char** argv) ParCG pcg_direct(pmat, smoother_direct, cg_max_iter); ParCG pcg_direct_blas(pmat, smoother_direct_blas, cg_max_iter); - SparseMatrix agg_vertex = make_agg_vertex(pmat); - ParMatrix par_agg_vertex(pmat.GetComm(), std::move(agg_vertex)); - ParBlockDiagComp block_comp(pmat, par_agg_vertex); - test_pcg_copy = cg; // rofl @@ -131,7 +126,6 @@ int main(int argc, char** argv) &pcg_direct_blas, &test_smooth_copy, &test_pcg_copy, - &block_comp, }; //ops = {&smoother_direct_blas}; 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/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/src/blockoperator.cpp b/src/blockoperator.cpp index 96d80a0..971631e 100644 --- a/src/blockoperator.cpp +++ b/src/blockoperator.cpp @@ -181,7 +181,7 @@ BlockDiagOperator::BlockDiagOperator(std::vector row_offsets, std::vector& input, VectorView output) const { assert(input.size() == cols_); diff --git a/tests/test_linalg.cpp b/tests/test_linalg.cpp index 1b30d70..f3c55e8 100644 --- a/tests/test_linalg.cpp +++ b/tests/test_linalg.cpp @@ -1551,6 +1551,17 @@ void test_assert() 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!"); @@ -1571,6 +1582,21 @@ void test_assert() 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!"); From 70cdee8f84bd318b7640b339b056d005a913b13a Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Thu, 13 Dec 2018 08:04:59 -0800 Subject: [PATCH 08/10] Remove a dot product in PCG --- src/solvers.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/solvers.cpp b/src/solvers.cpp index 1f617c3..3352204 100644 --- a/src/solvers.cpp +++ b/src/solvers.cpp @@ -189,7 +189,7 @@ void PCGSolver::Mult(const VectorView& b, VectorView x) const x.Add(alpha, p_); - double denom = (*Dot_)(z_, r_); + double denom = r_z; r_.Sub(alpha, Ap_); @@ -202,19 +202,19 @@ 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: %.2e\n", num_iter_, r_z_next / r0); } - 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_; From 8112b023b40829c2938f5692708b4fdc5bd00d9d Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Thu, 13 Dec 2018 08:10:12 -0800 Subject: [PATCH 09/10] Fix print statement in PCG --- modules/parlinalgcpp/include/parlinalgcpp.hpp | 1 - src/solvers.cpp | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/modules/parlinalgcpp/include/parlinalgcpp.hpp b/modules/parlinalgcpp/include/parlinalgcpp.hpp index 3e27cdf..59feb1a 100644 --- a/modules/parlinalgcpp/include/parlinalgcpp.hpp +++ b/modules/parlinalgcpp/include/parlinalgcpp.hpp @@ -5,6 +5,5 @@ #include "parsolvers.hpp" #include "parutilities.hpp" #include "paroperator.hpp" -#include "parprecond.hpp" #include "parsmoother.hpp" #include "parcommpkg.hpp" diff --git a/src/solvers.cpp b/src/solvers.cpp index 3f7c5e2..140a09a 100644 --- a/src/solvers.cpp +++ b/src/solvers.cpp @@ -226,8 +226,8 @@ void PCGSolver::Mult(const VectorView& b, VectorView x) const if (verbose_) { - printf("PCG %d: reduction: %.2e numer: %.2e, tol2: %.2e\n", - num_iter_, numer / r0, numer, tol_tol); + printf("PCG %d: reduction: %.2e r_z_next: %.2e, tol2: %.2e\n", + num_iter_, r_z_next / r0, r_z_next, tol_tol); } if (r_z_next < tol_tol) From 4df81717c19147e1a2db136d9bf0fcdabcae60fb Mon Sep 17 00:00:00 2001 From: Stephan Gelever Date: Thu, 13 Dec 2018 08:31:46 -0800 Subject: [PATCH 10/10] Fix compile error w/ templates --- include/utilities.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/utilities.hpp b/include/utilities.hpp index 8cccebb..4c23f4b 100644 --- a/include/utilities.hpp +++ b/include/utilities.hpp @@ -67,7 +67,7 @@ void linalgcpp_verify(F&& lambda, const std::string& message = "linalgcpp verifi @param size square size of identity @return identity matrix */ -template +template SparseMatrix SparseIdentity(int size); /** @brief Construct an rectangular identity matrix (as a SparseMatrix) @@ -76,7 +76,7 @@ SparseMatrix SparseIdentity(int size); @param row_offset offset row where diagonal identity starts @param col_offset offset column where diagonal identity starts */ -template +template SparseMatrix SparseIdentity(int rows, int cols, int row_offset = 0, int col_offset = 0); template