libMesh
sparse_matrix_test.h
Go to the documentation of this file.
1 #ifndef SPARSE_MATRIX_TEST_H
2 #define SPARSE_MATRIX_TEST_H
3 
4 // test includes
5 #include "test_comm.h"
6 
7 // libMesh includes
8 #include <libmesh/dense_matrix.h>
9 #include <libmesh/sparse_matrix.h>
10 #include <libmesh/fuzzy_equals.h>
11 
12 #include "libmesh_cppunit.h"
13 
14 // C++ includes
15 #include <vector>
16 #ifdef LIBMESH_HAVE_CXX11_THREAD
17 #include <thread>
18 #include <algorithm>
19 #endif
20 
21 
22 #define SPARSEMATRIXTEST \
23  CPPUNIT_TEST(testGetAndSet); \
24  CPPUNIT_TEST(testReadMatlab1); \
25  CPPUNIT_TEST(testReadMatlab2); \
26  CPPUNIT_TEST(testReadMatlab4); \
27  CPPUNIT_TEST(testReadHDF5); \
28  CPPUNIT_TEST(testWriteAndRead); \
29  CPPUNIT_TEST(testClone);
30 
31 
32 
33 template <class DerivedClass>
34 class SparseMatrixTest : public CppUnit::TestCase
35 {
36 public:
37 
38  void setUp()
39  {
40  // By default we'll use the whole communicator in parallel;
41  // Serial-only NumericVector subclasses will need to override
42  // this and set something else first.
43  if (!my_comm)
45 
46  matrix = std::make_unique<DerivedClass>(*my_comm);
47 
48  // Use the same even partitioning that we'll auto-deduce in matrix
49  // read, to make it easier to test parallel reads
50 
51  global_m = global_n = my_comm->size() * 5.75 - 1;
52  if (nonsquare)
53  global_n *= 1.3;
54 
56  row_start = my_comm->rank() * global_m / my_comm->size(),
57  row_stop = (my_comm->rank()+1) * global_m / my_comm->size();
59  col_start = my_comm->rank() * global_n / my_comm->size(),
60  col_stop = (my_comm->rank()+1) * global_n / my_comm->size();
61 
62  local_m = row_stop - row_start;
63  local_n = col_stop - col_start;
64 
65  // Let's just play around with locally dense blocks for now
66  matrix->init(global_m,
67  global_n,
68  local_m,
69  local_n,
70  /*nnz=*/local_n,
71  /*noz=*/5);
72  }
73 
74 
75  void tearDown() {}
76 
77 
78  void setValues()
79  {
80  std::vector<libMesh::numeric_index_type> rows(local_m);
81  std::iota(rows.begin(), rows.end(), matrix->row_start());
82  std::vector<libMesh::numeric_index_type> cols(local_n);
83  std::iota(cols.begin(), cols.end(), matrix->col_start());
84 
86 
87  for (auto i : libMesh::make_range(local_m))
88  for (auto j : libMesh::make_range(local_n))
89  local(i, j) = (i + 1) * (j + 1) * (my_comm->rank() + 1);
90 
91  matrix->zero();
92 
93  matrix->add_matrix(local, rows, cols);
94  matrix->close();
95  }
96 
97 
98  void testValues()
99  {
100  auto functor = [this]()
101  {
102  std::vector<libMesh::numeric_index_type> cols_to_get;
103  std::vector<libMesh::Number> values;
104  const libMesh::numeric_index_type col_start =
105  matrix->col_start();
107  libMesh::make_range(matrix->row_start(),
108  matrix->row_stop()))
109  {
110  matrix->get_row(i, cols_to_get, values);
111  for (libMesh::numeric_index_type col_j :
112  libMesh::index_range(cols_to_get))
113  {
114  CPPUNIT_ASSERT_EQUAL(cols_to_get[col_j], col_start + col_j);
115  LIBMESH_ASSERT_NUMBERS_EQUAL
116  ((i - matrix->row_start() + 1) * (col_j + 1) * (my_comm->rank() + 1),
117  values[col_j], _tolerance);
118  }
119  }
120  };
121 
122 #ifdef LIBMESH_HAVE_CXX11_THREAD
123  auto num_threads = std::min(unsigned(2),
124  std::max(
125  std::thread::hardware_concurrency(),
126  unsigned(1)));
127  std::vector<std::thread> threads(num_threads);
128  for (unsigned int thread = 0; thread < num_threads; ++thread)
129  threads[thread] = std::thread(functor);
130  std::for_each(threads.begin(), threads.end(),
131  [](std::thread & x){x.join();});
132 #else
133  functor();
134 #endif
135  }
136 
137 
139  {
140  LOG_UNIT_TEST;
141 
142  setValues();
143 
144  testValues();
145  }
146 
147  void testReadMatlab(const std::string & filename)
148  {
149  // Laspack doesn't handle non-square matrices)
150  if (matrix->solver_package() == libMesh::LASPACK_SOLVERS)
151  return;
152 
153  matrix->clear();
154 
155  auto matrix2 = std::make_unique<DerivedClass>(*my_comm);
156 
157  matrix->read(filename);
158 
159 #ifndef LIBMESH_HAVE_GZSTREAM
160  return;
161 #endif
162 
163  matrix2->read(std::string(filename)+".gz");
164 
165  // We need some more SparseMatrix operators, but not today
166  CPPUNIT_ASSERT(matrix->l1_norm() == matrix2->l1_norm());
167  CPPUNIT_ASSERT(matrix->linfty_norm() == matrix2->linfty_norm());
168  }
169 
170 
172  {
173  LOG_UNIT_TEST;
174  testReadMatlab("matrices/geom_1_extraction_op.m");
175  }
176 
177 
179  {
180  LOG_UNIT_TEST;
181  testReadMatlab("matrices/geom_2_extraction_op.m");
182  }
183 
184 
186  {
187  LOG_UNIT_TEST;
188  testReadMatlab("matrices/geom_4_extraction_op.m");
189  }
190 
191 
193  {
194 #ifdef LIBMESH_HAVE_HDF5
195  LOG_UNIT_TEST;
196 
197  matrix->clear();
198  auto matrix2 = std::make_unique<DerivedClass>(*my_comm);
199  matrix->read("matrices/geom_1_extraction_op.m");
200  matrix2->read("matrices/geom_1_extraction_op.h5");
201 
202  // We need some more SparseMatrix operators, but not today
203  CPPUNIT_ASSERT(matrix->l1_norm() == matrix2->l1_norm());
204  CPPUNIT_ASSERT(matrix->linfty_norm() == matrix2->linfty_norm());
205 #endif // LIBMESH_HAVE_HDF5
206  }
207 
208 
210  {
211  LOG_UNIT_TEST;
212 
213  setValues();
214 
215  // If we're working with serial matrices then just print one of
216  // them so they don't step on the others' toes.
217  //
218  // Use a very short filename, because we had a bug with that and
219  // we want to test it.
220  if (matrix->n_processors() > 1 ||
221  TestCommWorld->rank() == 0)
222  matrix->print_matlab("M.m");
223 
224  matrix->clear();
225 
226  // Let's make sure we don't have any race conditions; we have
227  // multiple SparseMatrix subclasses that might be trying to read
228  // and write the same file.
229 
231 
232 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
233  // We're not supporting complex reads quite yet
234  return;
235 #endif
236 
237  matrix->read("M.m");
238 
240 
241  testValues();
242  }
243 
244  void testClone()
245  {
246  LOG_UNIT_TEST;
247 
248  setValues();
249 
250  // Matrix must be closed before it can be cloned.
251  matrix->close();
252 
253  {
254  // Create copy, test that it can go out of scope
255  auto copy = matrix->clone();
256 
257  // Check that matrices have the same local/global sizes
258  CPPUNIT_ASSERT_EQUAL(copy->m(), matrix->m());
259  CPPUNIT_ASSERT_EQUAL(copy->n(), matrix->n());
260  CPPUNIT_ASSERT_EQUAL(copy->local_m(), matrix->local_m());
261  CPPUNIT_ASSERT_EQUAL(copy->row_start(), matrix->row_start());
262  CPPUNIT_ASSERT_EQUAL(copy->row_stop(), matrix->row_stop());
263 
264  // Check that copy has same values as original
265  LIBMESH_ASSERT_FP_EQUAL(copy->l1_norm(), matrix->l1_norm(), _tolerance);
266  CPPUNIT_ASSERT(relative_fuzzy_equals(*matrix, *copy));
267  copy->scale(2);
268  CPPUNIT_ASSERT(!relative_fuzzy_equals(*matrix, *copy));
269  }
270 
271  {
272  // Create zero copy
273  auto zero_copy = matrix->zero_clone();
274 
275  // Check that matrices have the same local/global sizes
276  CPPUNIT_ASSERT_EQUAL(zero_copy->m(), matrix->m());
277  CPPUNIT_ASSERT_EQUAL(zero_copy->n(), matrix->n());
278  CPPUNIT_ASSERT_EQUAL(zero_copy->local_m(), matrix->local_m());
279  CPPUNIT_ASSERT_EQUAL(zero_copy->row_start(), matrix->row_start());
280  CPPUNIT_ASSERT_EQUAL(zero_copy->row_stop(), matrix->row_stop());
281 
282  // Check that zero_copy has same values as original
283  LIBMESH_ASSERT_FP_EQUAL(0.0, zero_copy->l1_norm(), _tolerance);
284  }
285  }
286 
287 protected:
288 
289  std::string libmesh_suite_name;
290 
292 
293  std::unique_ptr<DerivedClass> matrix;
294 
296  local_m, local_n,
298 
300 };
301 
302 #endif // SPARSE_MATRIX_TEST_H
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
static constexpr Real TOLERANCE
void barrier() const
libMesh::Parallel::Communicator * my_comm
processor_id_type rank() const
libMesh::numeric_index_type local_n
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
processor_id_type size() const
dof_id_type numeric_index_type
Definition: id_types.h:99
void testReadMatlab(const std::string &filename)
std::unique_ptr< DerivedClass > matrix
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
libMesh::numeric_index_type global_m
std::string libmesh_suite_name
const libMesh::Real _tolerance
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
libMesh::numeric_index_type nonsquare
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
Definition: fuzzy_equals.h:78
libMesh::numeric_index_type global_n
Defines a dense matrix for use in Finite Element-type computations.
Definition: dof_map.h:75
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
libMesh::numeric_index_type local_m
uint8_t dof_id_type
Definition: id_types.h:67