diff --git a/src/fvm.hpp b/src/fvm.hpp
index e72bd5e82cf0098c22fc5d92d20c61d1a76e0ba9..bd1339c01c09539e52539c418728fe9952863486 100644
--- a/src/fvm.hpp
+++ b/src/fvm.hpp
@@ -2,14 +2,14 @@
#include <algorithm>
-#include "cell.hpp"
-#include "segment.hpp"
+#include <cell.hpp>
+#include <segment.hpp>
+#include <math.hpp>
+#include <matrix.hpp>
+#include <util.hpp>
+#include <algorithms.hpp>
-#include "math.hpp"
-#include "util.hpp"
-#include "algorithms.hpp"
-
-#include "../vector/include/Vector.hpp"
+#include <vector/include/Vector.hpp>
namespace nest {
namespace mc {
@@ -33,7 +33,7 @@ class fvm_cell {
private:
- index_type parent_index_;
+ matrix<value_type, size_type> matrix_;
};
////////////////////////////////////////////////////////////
@@ -42,8 +42,8 @@ class fvm_cell {
template <typename T, typename I>
fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
+: matrix_(cell.parent_index())
{
- parent_index_ = cell.parent_index();
auto const& segment_index = cell.segment_index();
auto num_fv = segment_index.back();
@@ -69,7 +69,7 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
for(auto c : cable->compartments()) {
auto rhs = segment_index[seg_idx] + c.index;
- auto lhs = parent_index_[rhs];
+ auto lhs = matrix_.p()[rhs];
auto rad_C = math::mean(left(c.radius), right(c.radius));
auto len = c.length/2;
@@ -89,8 +89,10 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
++seg_idx;
}
- //std::cout << "areas " << algorithms::sum(fv_volumes) << ", " << cell.volume() << "\n";
- //std::cout << "volumes " << algorithms::sum(fv_areas) << ", " << cell.area() << "\n";
+ for(auto i=0; i<num_fv; ++i) {
+ //auto area = fv_areas[i];
+ }
+
}
} // namespace fvm
diff --git a/src/matrix.hpp b/src/matrix.hpp
index 1b1daad2fc17af2d8a3ad87ced5fc03435deca96..cd8fc30a9c18ebb90e846c0510cca7ac81855560 100644
--- a/src/matrix.hpp
+++ b/src/matrix.hpp
@@ -1,6 +1,10 @@
#pragma once
-#include "../vector/include/Vector.hpp"
+#include <type_traits>
+
+#include "util.hpp"
+
+#include <vector/include/Vector.hpp>
namespace nest {
namespace mc {
@@ -20,67 +24,144 @@ class matrix {
using index_type = memory::HostVector<size_type>;
using index_view = typename index_type::view_type;
+ /// construct matrix for one or more cells, with combined parent index
+ /// and a cell index
+ template <
+ typename LHS, typename RHS,
+ typename = typename
+ std::enable_if<
+ util::is_container<LHS>::value &&
+ util::is_container<RHS>::value
+ >
+ >
+ matrix(LHS&& pi, RHS&& ci)
+ : parent_index_(std::forward<LHS>(pi))
+ , cell_index_(std::forward<RHS>(ci))
+ {
+ setup();
+ }
+
+ /// construct matrix for a single cell described by a parent index
+ template <
+ typename IDX,
+ typename = typename
+ std::enable_if< util::is_container<IDX>::value >
+ >
+ matrix(IDX&& pi)
+ : parent_index_(std::forward<IDX>(pi))
+ , cell_index_(2)
+ {
+ cell_index_[0] = 0;
+ cell_index_[1] = size();
+ setup();
+ }
/// the dimension of the matrix (i.e. the number of rows or colums)
size_type size() const
{
- return parent_indices_.size();
+ return parent_index_.size();
}
/// the total memory used to store the matrix
std::size_t memory() const
{
- auto s = sizeof(value_type) * data_.size();
- s += sizeof(size_type) * parent_indices_.size();
- return s + sizeof(matrix);
+ auto s = 5 * (sizeof(value_type) * size() + sizeof(vector_type));
+ s += sizeof(size_type) * (parent_index_.size() + cell_index_.size())
+ + 2*sizeof(index_type);
+ s += sizeof(matrix);
+ return s;
}
/// the number of cell matrices that have been packed together
size_type num_cells() const
{
- return cell_index_.size();
+ return cell_index_.size() - 1;
}
- private:
+ /// the vector holding the lower part of the matrix
+ view_type l()
+ {
+ return l_;
+ }
+
+ /// the vector holding the diagonal of the matrix
+ view_type d()
+ {
+ return d_;
+ }
+
+ /// the vector holding the upper part of the matrix
+ view_type u()
+ {
+ return u_;
+ }
+
+ /// the vector holding the right hand side of the linear equation system
+ view_type rhs()
+ {
+ return rhs_;
+ }
+
+ /// the vector holding the parent index
+ index_view p()
+ {
+ return parent_index_;
+ }
- /// calculate the number of values per sub-array plus padding
- size_type reservation() const
+ /// solve the linear system
+ /// upon completion the solution is stored in the RHS storage
+ /// and can be accessed via rhs()
+ void solve()
{
- constexpr auto alignment = vector_type::coordinator_type::alignment();
- auto const padding = memory::impl::get_padding<value_type>(alignment, size());
- return size() + padding;
+ index_view const& p = parent_index_;
+ auto const ncells = num_cells();
+
+ // loop over submatrices
+ for(auto m=0; m<ncells; ++m) {
+ auto first = cell_index_[m];
+ auto last = cell_index_[m+1];
+
+ // backward sweep
+ for(auto i=last-1; i>first; --i) {
+ auto factor = l_[i] / d_[i];
+ d_[p[i]] -= factor * u_[i];
+ rhs_[p[i]] -= factor * rhs_[i];
+ }
+ rhs_[first] /= d_[first];
+
+ // forward sweep
+ for(auto i=first+1; i<last; ++i) {
+ rhs_[i] -= u_[i] * rhs_[p[i]];
+ rhs_[i] /= d_[i];
+ }
+ }
}
- /// subdivide the memory in data_ into the sub-arrays that store the matrix
- /// diagonals and rhs/solution vectors
+ private:
+
+ /// allocate memory for storing matrix and right hand side vector
void setup()
{
const auto n = size();
- const auto n_alloc = reservation();
-
- // point sub-vectors into storage
- a_ = data_( 0, n);
- d_ = data_( n_alloc, n_alloc + n);
- b_ = data_(2*n_alloc, 2*n_alloc + n);
- rhs_ = data_(3*n_alloc, 3*n_alloc + n);
- v_ = data_(4*n_alloc, 4*n_alloc + n);
- }
- // the parent indices are stored as an index array
- index_type parent_indices_;
- index_type cell_index_;
+ l_ = vector_type(n);
+ d_ = vector_type(n);
+ u_ = vector_type(n);
+ rhs_ = vector_type(n);
+ }
- size_type num_cells_;
+ /// the parent indice that describe matrix structure
+ index_type parent_index_;
- // all the data fields point into a single block of storage
- vector_type data_;
+ /// indexes that point to the start of each cell in the index
+ index_type cell_index_;
- // the individual data files that point into block storage
- view_type a_;
- view_type d_;
- view_type b_;
- view_type rhs_;
- view_type v_;
+ /// storage for lower, diagonal, upper and rhs
+ vector_type l_;
+ vector_type d_;
+ vector_type u_;
+ /// after calling solve, the solution is stored in rhs_
+ vector_type rhs_;
};
} // namespace nest
diff --git a/src/util.hpp b/src/util.hpp
index cb68c2d1f7aa05c5ad4da7a37a5a1abae11c19dc..d92392799cf98da27ee1ca62788686a308c02c1d 100644
--- a/src/util.hpp
+++ b/src/util.hpp
@@ -65,6 +65,21 @@ namespace util {
return p.second;
}
+ template <typename T>
+ struct is_std_vector : std::false_type {};
+
+ template <typename T, typename C>
+ struct is_std_vector<std::vector<T,C>> : std::true_type {};
+
+ template <typename T>
+ struct is_container :
+ std::conditional<
+ is_std_vector<typename std::decay<T>::type>::value ||
+ memory::is_array<T>::value,
+ std::true_type,
+ std::false_type
+ >::type
+ {};
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 15d3ec160e9f9a701564c802ee188ebd8845b48e..b555fc465f3e3292ae2fa48f11dfad17bbe0e6e5 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -13,10 +13,11 @@ set(TEST_SOURCES
gtest-all.cpp
# unit tests
+ test_algorithms.cpp
test_cell.cpp
test_compartments.cpp
+ test_matrix.cpp
test_point.cpp
- test_algorithms.cpp
test_run.cpp
test_segment.cpp
test_stimulus.cpp
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index decf1316702dc5e04133131e838caa09e2539373..9328b315f992a9b6b12dd9e88b0e13ef6ace5d22 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -2,8 +2,8 @@
#include <vector>
-#include "../src/algorithms.hpp"
-#include "../src/util.hpp"
+#include <algorithms.hpp>
+#include <util.hpp>
TEST(algorithms, sum)
{
diff --git a/tests/test_matrix.cpp b/tests/test_matrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ba6794f1bd4ab095cc542bee340bb8a4e5ff15cc
--- /dev/null
+++ b/tests/test_matrix.cpp
@@ -0,0 +1,123 @@
+#include <numeric>
+
+#include "gtest.h"
+
+#include <matrix.hpp>
+#include <math.hpp>
+#include <vector/include/Vector.hpp>
+
+TEST(matrix, construct_from_parent_only)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // pass parent index as a std::vector by reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(p.size(), 3);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a std::vector by rvalue reference
+ // the input should not be moved from
+ {
+ std::vector<int> p = {0,0,1};
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(p.size(), 3);
+ EXPECT_EQ(m.size(), 3);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+
+ // pass parent index as a HostVector by reference
+ // expect that the input is not moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(p.size(), 3);
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+ // pass parent index as a HostVector by rvalue reference
+ // expect that the input is moved from
+ {
+ memory::HostVector<int> p(3, 0);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{std::move(p)};
+ EXPECT_EQ(m.num_cells(), 1);
+ EXPECT_EQ(m.size(), 3);
+ EXPECT_EQ(p.size(), 0); // 0 implies moved from
+
+ auto mp = m.p();
+ EXPECT_EQ(mp[0], 0);
+ EXPECT_EQ(mp[1], 0);
+ EXPECT_EQ(mp[2], 1);
+ }
+}
+
+TEST(matrix, solve)
+{
+ using matrix_type = nest::mc::matrix<double, int>;
+
+ // trivial case : 1x1 matrix
+ {
+ matrix_type m{std::vector<int>{0}};
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ EXPECT_EQ(m.rhs()[0], 0.5);
+ }
+ // matrices in the range of 2x2 to 1000x1000
+ {
+ using namespace nest::mc::math;
+ for(auto n : memory::Range(2,1001)) {
+ auto p = std::vector<int>(n);
+ std::iota(p.begin()+1, p.end(), 0);
+ matrix_type m{p};
+
+ EXPECT_EQ(m.size(), n);
+ EXPECT_EQ(m.num_cells(), 1);
+
+ m.d()(memory::all) = 2;
+ m.l()(memory::all) = -1;
+ m.u()(memory::all) = -1;
+ m.rhs()(memory::all) = 1;
+
+ m.solve();
+
+ auto x = m.rhs();
+ auto err = square(std::fabs(2.*x[0] - x[1] - 1.));
+ for(auto i : memory::Range(1,n-1)) {
+ err += square(std::fabs(2.*x[i] - x[i-1] - x[i+1] - 1.));
+ }
+ err += square(std::fabs(2.*x[n-1] - x[n-2] - 1.));
+
+ EXPECT_NEAR(0., std::sqrt(err), 1e-8);
+ }
+
+ }
+}
+
diff --git a/vector b/vector
index c2fc6318f6d65ae9d7aad60309f3b53d590a00ec..e95bcbf1de6f833b5113a82dcc6d47bb7f41c397 160000
--- a/vector
+++ b/vector
@@ -1 +1 @@
-Subproject commit c2fc6318f6d65ae9d7aad60309f3b53d590a00ec
+Subproject commit e95bcbf1de6f833b5113a82dcc6d47bb7f41c397