Develop

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}"

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";
 

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;

examples/include/ex_utilities.hpp

@@ -3,12 +3,7 @@
 
 #include <iostream>
 
-#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<double>;
 using VectorView = linalgcpp::VectorView<double>;
@@ -17,6 +12,8 @@ using SparseMatrix = linalgcpp::SparseMatrix<double>;
 using CooMatrix = linalgcpp::CooMatrix<double>;
 using Timer = linalgcpp::Timer;
 using Operator = linalgcpp::Operator;
+using linalgcpp::linalgcpp_assert;
+using linalgcpp::linalgcpp_verify;
 
 
 inline

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}"

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}"

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<double> 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<int> offsets);
+
+        /*! @brief Rectangle Constructor with given offsets*/
+        BlockDiagOperator(std::vector<int> row_offsets, std::vector<int> 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<int>& GetRowOffsets() const;
+
+        /*! @brief Get the col offsets
+            @retval the col offsets
+        */
+        const std::vector<int>& 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<double>& input, VectorView<double> output) const override;
+
+        /// Operator Requirement
+        void MultAT(const VectorView<double>& input, VectorView<double> output) const override;
+
+        using Operator::Mult;
+
+    private:
+        std::vector<int> row_offsets_;
+        std::vector<int> col_offsets_;
+
+        std::vector<const Operator*> A_;
+
+        mutable BlockVector<double> x_;
+        mutable BlockVector<double> y_;
+
+        mutable Vector<double> tmp_;
+};
+
 } //namespace linalgcpp
 
 #endif // BLOCKOPERATOR_HPP__

include/linalgcpp.hpp

@@ -17,3 +17,4 @@
 #include "eigensolver.hpp"
 #include "timer.hpp"
 #include "argparser.hpp"
+#include "utilities.hpp"

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;

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<double>& input, VectorView<double> output) const override;
@@ -1470,7 +1471,7 @@ void SparseMatrix<T>::EliminateCol(const U& marker, const VectorView<double>& x,
 }
 
 template <typename T>
-void SparseMatrix<T>::EliminateZeros(T tolerance)
+void SparseMatrix<T>::EliminateZeros(T tolerance, bool keep_diag)
 {
     std::vector<int> elim_indptr(1, 0);
     std::vector<int> elim_indices;
@@ -1484,7 +1485,10 @@ void SparseMatrix<T>::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]);

include/utilities.hpp

@@ -3,14 +3,27 @@
 #ifndef UTILITIES_HPP__
 #define UTILITIES_HPP__
 
+#include <memory>
+#include "sparsematrix.hpp"
+
+#if __cplusplus > 201103L
+using std::make_unique;
+#else
+template<typename T, typename... Ts>
+std::unique_ptr<T> make_unique(Ts&& ... params)
+{
+    return std::unique_ptr<T>(new T(std::forward<Ts>(params)...));
+}
+#endif
+
 namespace linalgcpp
 {
 
 /*! @brief Throw if false in debug mode only */
 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 <typename F>
+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 <typename F>
+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 <typename T>
+SparseMatrix<T> SparseIdentity(int size);
+
+/** @brief Construct an rectangular identity matrix (as a SparseMatrix)
+    @param rows number of row
+    @param cols number of columns
+    @param row_offset offset row where diagonal identity starts
+    @param col_offset offset column where diagonal identity starts
+*/
+template <typename T>
+SparseMatrix<T> SparseIdentity(int rows, int cols, int row_offset = 0, int col_offset = 0);
+
+template <typename T = double>
+SparseMatrix<T> SparseIdentity(int size)
+{
+    assert(size >= 0);
+
+    return SparseMatrix<T>(std::vector<T>(size, (T)1.0));
+}
+
+template <typename T = double>
+SparseMatrix<T> 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<int> 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<int> indices(diag_size);
+    std::iota(std::begin(indices), std::begin(indices) + diag_size, col_offset);
+
+    std::vector<T> data(diag_size, 1.0);
+
+    return SparseMatrix<T>(std::move(indptr), std::move(indices), std::move(data), rows, cols);
+}
+
 } //namespace linalgcpp
 
 #endif // UTILITIES_HPP__