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