1 #include <libmesh/petsc_matrix.h>
2 
3 #ifdef LIBMESH_HAVE_PETSC
4 
5 #include <libmesh/parallel.h>
6 #include <libmesh/dense_matrix.h>
7 
8 #include "libmesh_cppunit.h"
9 #include "test_comm.h"
10 
11 #include <vector>
12 #ifdef LIBMESH_HAVE_CXX11_THREAD
13 #include <thread>
14 #include <algorithm>
15 #endif
16 
17 using namespace libMesh;
18 
19 class PetscMatrixTest : public CppUnit::TestCase
20 {
21 public:
22   CPPUNIT_TEST_SUITE(PetscMatrixTest);
23 
24   CPPUNIT_TEST(testGetAndSet);
25   CPPUNIT_TEST(testClone);
26 
27   CPPUNIT_TEST_SUITE_END();
28 
29 public:
setUp()30   void setUp()
31   {
32     _comm = TestCommWorld;
33     _matrix = libmesh_make_unique<PetscMatrix<Number>>(*_comm);
34 
35     numeric_index_type root_block_size = 2;
36     _local_size = root_block_size + static_cast<numeric_index_type>(_comm->rank());
37     _global_size = 0;
38     _i.push_back(0);
39 
40     for (processor_id_type p = 0; p < _comm->size(); ++p)
41     {
42       numeric_index_type block_size = root_block_size + static_cast<numeric_index_type>(p);
43       _global_size += block_size;
44       _i.push_back(_global_size);
45     }
46 
47     _matrix->init(_global_size,
48                   _global_size,
49                   _local_size,
50                   _local_size,
51                   /*nnz=*/_local_size,
52                   /*noz=*/0);
53   }
54 
tearDown()55   void tearDown() {}
56 
testGetAndSet()57   void testGetAndSet()
58   {
59     std::vector<numeric_index_type> rows(_local_size);
60     std::vector<numeric_index_type> cols(_local_size);
61     DenseMatrix<Number> local(_local_size, _local_size);
62 
63     numeric_index_type index = _i[_comm->rank()], count = 0;
64     for (; count < _local_size; ++count, ++index)
65     {
66       rows[count] = index;
67       cols[count] = index;
68       for (numeric_index_type j = 0; j < _local_size; ++j)
69         local(count, j) = (count + 1) * (j + 1) * (_comm->rank() + 1);
70     }
71 
72     _matrix->add_matrix(local, rows, cols);
73     _matrix->close();
74 
75     index = _i[_comm->rank()], count = 0;
76 
77     auto functor = [this]()
78     {
79       std::vector<numeric_index_type> cols_to_get;
80       std::vector<Number> values;
81       numeric_index_type local_index = _i[_comm->rank()], local_count = 0;
82       for (; local_count < _local_size; ++local_count, ++local_index)
83       {
84         _matrix->get_row(local_index, cols_to_get, values);
85         for (numeric_index_type j = 0; j < _local_size; ++j)
86         {
87           LIBMESH_ASSERT_FP_EQUAL((local_count + 1) * (j + 1) * (_comm->rank() + 1),
88                                   libMesh::libmesh_real(values[j]),
89                                   _tolerance);
90           CPPUNIT_ASSERT_EQUAL(cols_to_get[local_count], local_index);
91         }
92       }
93     };
94 
95 #ifdef LIBMESH_HAVE_CXX11_THREAD
96     auto num_threads = std::min(unsigned(2),
97                                 std::max(
98                                   std::thread::hardware_concurrency(),
99                                   unsigned(1)));
100     std::vector<std::thread> threads(num_threads);
101     for (unsigned int thread = 0; thread < num_threads; ++thread)
102       threads[thread] = std::thread(functor);
103     std::for_each(threads.begin(), threads.end(),
104                   [](std::thread & x){x.join();});
105 #else
106     functor();
107 #endif
108   }
109 
testClone()110   void testClone()
111   {
112     // Matrix must be closed before it can be cloned.
113     _matrix->close();
114 
115     {
116       // Create copy, test that it can go out of scope
117       auto copy = _matrix->clone();
118 
119       // Check that matrices have the same local/global sizes
120       CPPUNIT_ASSERT_EQUAL(copy->m(), _matrix->m());
121       CPPUNIT_ASSERT_EQUAL(copy->n(), _matrix->n());
122       CPPUNIT_ASSERT_EQUAL(copy->local_m(), _matrix->local_m());
123       CPPUNIT_ASSERT_EQUAL(copy->row_start(), _matrix->row_start());
124       CPPUNIT_ASSERT_EQUAL(copy->row_stop(), _matrix->row_stop());
125 
126       // Check that copy has same values as original
127       LIBMESH_ASSERT_FP_EQUAL(copy->l1_norm(), _matrix->l1_norm(), _tolerance);
128     }
129 
130     {
131       // Create zero copy
132       auto zero_copy = _matrix->zero_clone();
133 
134       // Check that matrices have the same local/global sizes
135       CPPUNIT_ASSERT_EQUAL(zero_copy->m(), _matrix->m());
136       CPPUNIT_ASSERT_EQUAL(zero_copy->n(), _matrix->n());
137       CPPUNIT_ASSERT_EQUAL(zero_copy->local_m(), _matrix->local_m());
138       CPPUNIT_ASSERT_EQUAL(zero_copy->row_start(), _matrix->row_start());
139       CPPUNIT_ASSERT_EQUAL(zero_copy->row_stop(), _matrix->row_stop());
140 
141       // Check that zero_copy has same values as original
142       LIBMESH_ASSERT_FP_EQUAL(0.0, zero_copy->l1_norm(), _tolerance);
143     }
144   }
145 
146 private:
147 
148   Parallel::Communicator * _comm;
149   std::unique_ptr<PetscMatrix<Number>> _matrix;
150   numeric_index_type _local_size, _global_size;
151   std::vector<numeric_index_type> _i;
152   const Real _tolerance = TOLERANCE * TOLERANCE;
153 
154 };
155 
156 CPPUNIT_TEST_SUITE_REGISTRATION(PetscMatrixTest);
157 
158 #endif
159