From 2c48dc9a9a583795913e8d0cbf8b33d6665cd0d0 Mon Sep 17 00:00:00 2001 From: root Date: Wed, 26 Jul 2023 03:16:07 +0000 Subject: [PATCH 1/5] refactor sparse_matrix class to use map --- source/module_io/sparse_matrix.cpp | 58 +++++---- source/module_io/sparse_matrix.h | 55 ++++++-- source/module_io/test/sparse_matrix_test.cpp | 130 +++++++++++-------- 3 files changed, 148 insertions(+), 95 deletions(-) diff --git a/source/module_io/sparse_matrix.cpp b/source/module_io/sparse_matrix.cpp index 1781d09bcf3..b7585bd478a 100644 --- a/source/module_io/sparse_matrix.cpp +++ b/source/module_io/sparse_matrix.cpp @@ -1,5 +1,6 @@ #include "sparse_matrix.h" +#include #include #include "module_base/tool_quit.h" @@ -11,51 +12,39 @@ namespace ModuleIO * @brief Add value to the matrix with row and column indices */ template -void SparseMatrix::addValue(int row, int col, T value) +void SparseMatrix::insert(int row, int col, T value) { if (row < 0 || row >= _rows || col < 0 || col >= _cols) { ModuleBase::WARNING_QUIT("SparseMatrix::addValue", "row or col index out of range"); } - data.push_back(std::make_tuple(row, col, value)); + if (std::abs(value) > _sparse_threshold) + { + elements[std::make_pair(row, col)] = value; + } } /** * @brief Print to CSR format */ template -void SparseMatrix::printToCSR(std::ostream& ofs, double threshold, int precision) +void SparseMatrix::printToCSR(std::ostream& ofs, int precision) { - // Filter elements greater than the threshold - auto it = std::remove_if(data.begin(), data.end(), [threshold](const std::tuple& elem) { - return std::abs(std::get<2>(elem)) <= threshold; - }); - data.erase(it, data.end()); - - // Sort data by row, then by column - std::sort(data.begin(), data.end(), [](const std::tuple& a, const std::tuple& b) { - if (std::get<0>(a) == std::get<0>(b)) - { - return std::get<1>(a) < std::get<1>(b); - } - return std::get<0>(a) < std::get<0>(b); - }); - // Initialize the CSR arrays std::vector csr_row_ptr; csr_row_ptr.assign(_rows + 1, 0); // print the CSR values - for (int i = 0; i < data.size(); i++) + for (const auto &element : elements) { - ofs << " " << std::fixed << std::scientific << std::setprecision(precision) << std::get<2>(data[i]); + ofs << " " << std::fixed << std::scientific << std::setprecision(precision) << element.second; } ofs << std::endl; // print the CSR column indices - for (int i = 0; i < data.size(); i++) + for (const auto &element : elements) { - ofs << " " << std::get<1>(data[i]); - int row = std::get<0>(data[i]); + ofs << " " << element.first.second; + int row = element.first.first; csr_row_ptr[row + 1]++; } ofs << std::endl; @@ -91,16 +80,35 @@ void SparseMatrix::readCSR(const std::vector& values, ModuleBase::WARNING_QUIT("SparseMatrix::readCSR", "Column indices and values size mismatch"); } - data.clear(); + elements.clear(); for (int row = 0; row < _rows; row++) { for (int idx = row_ptr[row]; idx < row_ptr[row + 1]; idx++) { - data.push_back(std::make_tuple(row, col_ind[idx], values[idx])); + elements[std::make_pair(row, col_ind[idx])] = values[idx]; } } } +// define the operator to index a matrix element +template +T SparseMatrix::operator()(int row, int col) +{ + if (row < 0 || row >= _rows || col < 0 || col >= _cols) + { + ModuleBase::WARNING_QUIT("SparseMatrix::operator()", "row or col index out of range"); + } + auto it = elements.find(std::make_pair(row, col)); + if (it != elements.end()) + { + return it->second; + } + else + { + return static_cast(0); + } +} + // Explicit instantiation of template classes template class SparseMatrix; template class SparseMatrix>; diff --git a/source/module_io/sparse_matrix.h b/source/module_io/sparse_matrix.h index 40ef07bdeca..63b7f6343a9 100644 --- a/source/module_io/sparse_matrix.h +++ b/source/module_io/sparse_matrix.h @@ -1,18 +1,23 @@ #ifndef SPARSE_MATRIX_H #define SPARSE_MATRIX_H -#include #include -#include +#include +#include #include namespace ModuleIO { /** - * @brief Sparse matrix class - * - * @tparam T + * @brief Sparse matrix class designed mainly for csr format input and output. + * @details + * The sparse matrix is stored in a map. + * The map key is a pair of row and column indices. + * The map value is the matrix element. + * The matrix element is stored only if its absolute value is greater than the threshold. + * The threshold is set to 1.0e-10 by default. + * @tparam T data type, it can be double or std::complex */ template class SparseMatrix @@ -28,16 +33,10 @@ class SparseMatrix } // add value to the matrix with row and column indices - void addValue(int row, int col, T value); - - // get data vector - const std::vector>& getData() const - { - return data; - } + void insert(int row, int col, T value); // print data in CSR (Compressed Sparse Row) format - void printToCSR(std::ostream& ofs, double threshold, int precision = 8); + void printToCSR(std::ostream& ofs, int precision = 8); // read CSR data from arrays void readCSR(const std::vector& values, const std::vector& col_ind, const std::vector& row_ptr); @@ -66,10 +65,38 @@ class SparseMatrix return _cols; } + // define the operator to index a matrix element + T operator()(int row, int col); + + // set the threshold + void setSparseThreshold(double sparse_threshold) + { + _sparse_threshold = sparse_threshold; + } + + // get the threshold + double getSparseThreshold() const + { + return _sparse_threshold; + } + + // get the number of non-zero elements + int getNNZ() const + { + return elements.size(); + } + + // get elements + const std::map, T>& getElements() const + { + return elements; + } + private: int _rows; int _cols; - std::vector> data; + std::map, T> elements; + double _sparse_threshold = 1.0e-10; }; // class SparseMatrix } // namespace ModuleIO diff --git a/source/module_io/test/sparse_matrix_test.cpp b/source/module_io/test/sparse_matrix_test.cpp index 28a6d5ec41f..28699c0d9d9 100644 --- a/source/module_io/test/sparse_matrix_test.cpp +++ b/source/module_io/test/sparse_matrix_test.cpp @@ -11,14 +11,12 @@ /** * - Tested Functions: - * - addValue() + * - insert() * - Add value to the matrix with row and column indices - * - convertToCSR() - * - Convert to CSR format - * - printCSR() - * - Print data in CSR format + * - printToCSR() + * - Print data in CSR format * - readCSR() - * - Read CSR data from arrays + * - Read CSR data from arrays */ class SparseMatrixTest : public testing::Test @@ -29,87 +27,89 @@ class SparseMatrixTest : public testing::Test std::string output; }; -TEST_F(SparseMatrixTest, AddValueDouble) +TEST_F(SparseMatrixTest, InsertDouble) { // Add a value to the matrix with row and column indices - smd.addValue(2, 2, 3.0); - smd.addValue(3, 3, 4.0); + smd.insert(2, 2, 3.0); + smd.insert(3, 3, 4.0); - EXPECT_EQ(smd.getData().size(), 2); + EXPECT_EQ(smd.getNNZ(), 2); EXPECT_EQ(smd.getCols(), 4); EXPECT_EQ(smd.getRows(), 4); // Check if the value was added correctly - EXPECT_EQ(std::make_tuple(2, 2, 3.0), smd.getData()[0]); - EXPECT_EQ(std::make_tuple(3, 3, 4.0), smd.getData()[1]); + EXPECT_DOUBLE_EQ(smd(2, 2), 3.0); + EXPECT_DOUBLE_EQ(smd(3, 3), 4.0); } -TEST_F(SparseMatrixTest, AddValueDoubleWarning) +TEST_F(SparseMatrixTest, InsertDoubleWarning) { // case 1 testing::internal::CaptureStdout(); - EXPECT_EXIT(smd.addValue(2, -1, 3.0), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smd.insert(2, -1, 3.0), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 2 testing::internal::CaptureStdout(); - EXPECT_EXIT(smd.addValue(-1, 0, 3.0), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smd.insert(-1, 0, 3.0), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 3 testing::internal::CaptureStdout(); - EXPECT_EXIT(smd.addValue(2, 4, 3.0), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smd.insert(2, 4, 3.0), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 4 testing::internal::CaptureStdout(); - EXPECT_EXIT(smd.addValue(4, 2, 3.0), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smd.insert(4, 2, 3.0), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); } -TEST_F(SparseMatrixTest, AddValueComplex) +TEST_F(SparseMatrixTest, InsertComplex) { // Add a value to the matrix with row and column indices - smc.addValue(2, 2, std::complex(1.0, 0.0)); - smc.addValue(3, 3, std::complex(0.0, 0.2)); + smc.insert(2, 2, std::complex(1.0, 0.0)); + smc.insert(3, 3, std::complex(0.0, 0.2)); - EXPECT_EQ(smc.getData().size(), 2); + EXPECT_EQ(smc.getNNZ(), 2); EXPECT_EQ(smc.getCols(), 4); EXPECT_EQ(smc.getRows(), 4); // Check if the value was added correctly - EXPECT_EQ(std::make_tuple(2, 2, std::complex(1.0, 0.0)), smc.getData()[0]); - EXPECT_EQ(std::make_tuple(3, 3, std::complex(0.0, 0.2)), smc.getData()[1]); + EXPECT_DOUBLE_EQ(smc(2, 2).real(), 1.0); + EXPECT_DOUBLE_EQ(smc(2, 2).imag(), 0.0); + EXPECT_DOUBLE_EQ(smc(3, 3).real(), 0.0); + EXPECT_DOUBLE_EQ(smc(3, 3).imag(), 0.2); } -TEST_F(SparseMatrixTest, AddValueComplexWarning) +TEST_F(SparseMatrixTest, InsertComplexWarning) { // case 1 testing::internal::CaptureStdout(); - EXPECT_EXIT(smc.addValue(2, -1, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smc.insert(2, -1, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 2 testing::internal::CaptureStdout(); - EXPECT_EXIT(smc.addValue(-1, 0, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smc.insert(-1, 0, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 3 testing::internal::CaptureStdout(); - EXPECT_EXIT(smc.addValue(2, 4, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smc.insert(2, 4, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); // case 4 testing::internal::CaptureStdout(); - EXPECT_EXIT(smc.addValue(4, 2, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); + EXPECT_EXIT(smc.insert(4, 2, std::complex(1.0, 0.0)), ::testing::ExitedWithCode(0), ""); output = testing::internal::GetCapturedStdout(); // test output on screen EXPECT_THAT(output, testing::HasSubstr("row or col index out of range")); @@ -118,7 +118,7 @@ TEST_F(SparseMatrixTest, AddValueComplexWarning) TEST_F(SparseMatrixTest, CSRDouble) { // test the printCSR function for the following matrix - // 1.0 2.0 0.0 0.0 0.0 0.0 + // 1.0 2.0 0.5 0.4 0.0 0.0 // 0.0 3.0 0.0 4.0 0.0 0.0 // 0.0 0.0 5.0 6.0 7.0 0.0 // 0.0 0.0 0.0 0.0 0.0 8.0 @@ -126,20 +126,25 @@ TEST_F(SparseMatrixTest, CSRDouble) // csr_values = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] // csr_col_ind = [0, 1, 1, 3, 2, 3, 4, 5] // csr_row_ptr = [0, 2, 4, 7, 8] + // set the sparse threshold to 0.6 and the precision to 2 ModuleIO::SparseMatrix smd0; ModuleIO::SparseMatrix smd1; smd0.setRows(4); smd0.setCols(6); - smd0.addValue(0, 0, 1.0); - smd0.addValue(0, 1, 2.0); - smd0.addValue(1, 1, 3.0); - smd0.addValue(1, 3, 4.0); - smd0.addValue(2, 2, 5.0); - smd0.addValue(2, 3, 6.0); - smd0.addValue(2, 4, 7.0); - smd0.addValue(3, 5, 8.0); + smd0.setSparseThreshold(0.6); + smd0.insert(0, 0, 1.0); + smd0.insert(0, 1, 2.0); + smd0.insert(0, 2, 0.5); + smd0.insert(0, 3, 0.4); + smd0.insert(1, 1, 3.0); + smd0.insert(1, 3, 4.0); + smd0.insert(2, 2, 5.0); + smd0.insert(2, 3, 6.0); + smd0.insert(2, 4, 7.0); + smd0.insert(3, 5, 8.0); + EXPECT_DOUBLE_EQ(smd0.getSparseThreshold(), 0.6); testing::internal::CaptureStdout(); - smd0.printToCSR(std::cout, 0.0, 2); + smd0.printToCSR(std::cout, 2); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("1.00e+00 2.00e+00 3.00e+00 4.00e+00 5.00e+00 6.00e+00 7.00e+00 8.00e+00")); EXPECT_THAT(output, testing::HasSubstr("0 1 1 3 2 3 4 5")); @@ -153,11 +158,13 @@ TEST_F(SparseMatrixTest, CSRDouble) smd1.setCols(6); smd1.readCSR(expected_csr_values, expected_csr_col_ind, expected_csr_row_ptr); - for (int i = 0; i < 8; i++) + EXPECT_EQ(smd0.getNNZ(), 8); + EXPECT_EQ(smd1.getNNZ(), 8); + + for (const auto& element : smd0.getElements()) { - EXPECT_EQ(std::get<0>(smd0.getData()[i]), std::get<0>(smd1.getData()[i])); - EXPECT_EQ(std::get<1>(smd0.getData()[i]), std::get<1>(smd1.getData()[i])); - EXPECT_DOUBLE_EQ(std::get<2>(smd0.getData()[i]), std::get<2>(smd1.getData()[i])); + auto it = smd1.getElements().find(element.first); + EXPECT_DOUBLE_EQ(it->second, element.second); } } @@ -172,20 +179,25 @@ TEST_F(SparseMatrixTest, CSRComplex) // csr_values = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0] // csr_col_ind = [0, 1, 1, 3, 2, 3, 4, 5] // csr_row_ptr = [0, 2, 4, 7, 8] + // set the sparse threshold to 0.6 and the precision to 1 ModuleIO::SparseMatrix> smc0; ModuleIO::SparseMatrix> smc1; smc0.setRows(4); smc0.setCols(6); - smc0.addValue(0, 0, std::complex(1.0, 0.0)); - smc0.addValue(0, 1, std::complex(2.0, 0.0)); - smc0.addValue(1, 1, std::complex(3.0, 0.0)); - smc0.addValue(1, 3, std::complex(4.0, 0.0)); - smc0.addValue(2, 2, std::complex(5.0, 0.0)); - smc0.addValue(2, 3, std::complex(6.0, 0.0)); - smc0.addValue(2, 4, std::complex(7.0, 0.0)); - smc0.addValue(3, 5, std::complex(8.0, 0.0)); + smc0.setSparseThreshold(0.6); + EXPECT_DOUBLE_EQ(smc0.getSparseThreshold(), 0.6); + smc0.insert(0, 0, std::complex(1.0, 0.0)); + smc0.insert(0, 1, std::complex(2.0, 0.0)); + smc0.insert(0, 2, std::complex(0.5, 0.0)); + smc0.insert(0, 3, std::complex(0.0, 0.4)); + smc0.insert(1, 1, std::complex(3.0, 0.0)); + smc0.insert(1, 3, std::complex(4.0, 0.0)); + smc0.insert(2, 2, std::complex(5.0, 0.0)); + smc0.insert(2, 3, std::complex(6.0, 0.0)); + smc0.insert(2, 4, std::complex(7.0, 0.0)); + smc0.insert(3, 5, std::complex(8.0, 0.0)); testing::internal::CaptureStdout(); - smc0.printToCSR(std::cout, 0.5, 1); + smc0.printToCSR(std::cout, 1); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("(1.0e+00,0.0e+00) (2.0e+00,0.0e+00) (3.0e+00,0.0e+00) (4.0e+00,0.0e+00) (5.0e+00,0.0e+00) (6.0e+00,0.0e+00) (7.0e+00,0.0e+00) (8.0e+00,0.0e+00)")); EXPECT_THAT(output, testing::HasSubstr("0 1 1 3 2 3 4 5")); @@ -206,11 +218,17 @@ TEST_F(SparseMatrixTest, CSRComplex) smc1.setCols(6); smc1.readCSR(expected_csr_values, expected_csr_col_ind, expected_csr_row_ptr); - for (int i = 0; i < 8; i++) + for (const auto& element : smc0.getElements()) { - EXPECT_EQ(std::get<0>(smc0.getData()[i]), std::get<0>(smc1.getData()[i])); - EXPECT_EQ(std::get<1>(smc0.getData()[i]), std::get<1>(smc1.getData()[i])); - EXPECT_DOUBLE_EQ((std::get<2>(smc0.getData()[i])).real(), (std::get<2>(smc1.getData()[i])).real()); - EXPECT_DOUBLE_EQ((std::get<2>(smc0.getData()[i])).imag(), (std::get<2>(smc1.getData()[i])).imag()); + auto it = smc1.getElements().find(element.first); + EXPECT_DOUBLE_EQ(it->second.real(), element.second.real()); + EXPECT_DOUBLE_EQ(it->second.imag(), element.second.imag()); } + // index matrix elements + EXPECT_DOUBLE_EQ(smc1(0, 1).imag(), 0.0); + EXPECT_DOUBLE_EQ(smc1(0, 1).real(), 2.0); + EXPECT_DOUBLE_EQ(smc1(1, 1).imag(), 0.0); + EXPECT_DOUBLE_EQ(smc1(1, 1).real(), 3.0); + EXPECT_DOUBLE_EQ(smc1(1, 2).imag(), 0.0); + EXPECT_DOUBLE_EQ(smc1(1, 2).real(), 0.0); } \ No newline at end of file From abf0b620b5119b43e4bf8345dd195dfa06fa34dd Mon Sep 17 00:00:00 2001 From: root Date: Wed, 26 Jul 2023 03:18:50 +0000 Subject: [PATCH 2/5] update output_hcontainer using new sparse matrix --- .../module_hcontainer/output_hcontainer.cpp | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp index d60421c8d54..9a2bbd570c3 100644 --- a/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/output_hcontainer.cpp @@ -66,7 +66,7 @@ void Output_HContainer::write_single_R(int rx, int ry, int rz) { this->_hcontainer->fix_R(rx, ry, rz); ModuleIO::SparseMatrix sparse_matrix = ModuleIO::SparseMatrix(this->_ParaV->nrow, this->_ParaV->ncol); - int nonzero = 0; + sparse_matrix.setSparseThreshold(this->_sparse_threshold); for (int iap = 0; iap < this->_hcontainer->size_atom_pairs(); ++iap) { auto atom_pair = this->_hcontainer->get_atom_pair(iap); @@ -86,22 +86,16 @@ void Output_HContainer::write_single_R(int rx, int ry, int rz) for (int iw2 = 0; iw2 < size2; ++iw2) { const int global_index2 = start2 + iw2; - T tmp_matrix_value = atom_pair.get_matrix_value(global_index1, global_index2); - - if (std::abs(tmp_matrix_value) > _sparse_threshold) - { - nonzero++; - sparse_matrix.addValue(global_index1, global_index2, tmp_matrix_value); - // to do: consider 2D block-cyclic distribution - } + sparse_matrix.insert(global_index1, global_index2, tmp_matrix_value); + // to do: consider 2D block-cyclic distribution } } } - if (nonzero != 0) + if (sparse_matrix.getNNZ() != 0) { - _ofs << rx << " " << ry << " " << rz << " " << nonzero << std::endl; - sparse_matrix.printToCSR(_ofs, _sparse_threshold, _precision); + _ofs << rx << " " << ry << " " << rz << " " << sparse_matrix.getNNZ() << std::endl; + sparse_matrix.printToCSR(_ofs, _precision); } this->_hcontainer->unfix_R(); } From 5a63d2c253f12948510ba08b4a78f9d66ab2faa4 Mon Sep 17 00:00:00 2001 From: root Date: Wed, 26 Jul 2023 04:18:08 +0000 Subject: [PATCH 3/5] add csr file reader --- source/module_io/CMakeLists.txt | 2 + source/module_io/csr_reader.cpp | 140 +++++++++++++++++++++ source/module_io/csr_reader.h | 71 +++++++++++ source/module_io/file_reader.cpp | 50 ++++++++ source/module_io/file_reader.h | 39 ++++++ source/module_io/test/CMakeLists.txt | 12 ++ source/module_io/test/csr_reader_test.cpp | 100 +++++++++++++++ source/module_io/test/file_reader_test.cpp | 73 +++++++++++ source/module_io/test/support/SR.csr | 11 ++ 9 files changed, 498 insertions(+) create mode 100644 source/module_io/csr_reader.cpp create mode 100644 source/module_io/csr_reader.h create mode 100644 source/module_io/file_reader.cpp create mode 100644 source/module_io/file_reader.h create mode 100644 source/module_io/test/csr_reader_test.cpp create mode 100644 source/module_io/test/file_reader_test.cpp create mode 100644 source/module_io/test/support/SR.csr diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 823f66512c9..da92f571c7d 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -51,6 +51,8 @@ if(ENABLE_LCAO) output_dm.cpp output_dm1.cpp sparse_matrix.cpp + file_reader.cpp + csr_reader.cpp ) list(APPEND objects_advanced unk_overlap_lcao.cpp diff --git a/source/module_io/csr_reader.cpp b/source/module_io/csr_reader.cpp new file mode 100644 index 00000000000..5d6ee372b52 --- /dev/null +++ b/source/module_io/csr_reader.cpp @@ -0,0 +1,140 @@ +#include "csr_reader.h" + +#include "module_base/tool_quit.h" + +namespace ModuleIO +{ + +// constructor +template +csrFileReader::csrFileReader(const std::string& filename) : FileReader(filename) +{ + parseFile(); +} + +// function to parse file +template +void csrFileReader::parseFile() +{ + // Check if file is open + if (!isOpen()) + { + ModuleBase::WARNING_QUIT("csrFileReader::parseFile", "File is not open"); + } + + std::string tmp_string; + + // Read the step + readLine(); + ss >> tmp_string >> step; + + // Read the matrix dimension + readLine(); + ss >> tmp_string >> tmp_string >> tmp_string >> tmp_string >> matrixDimension; + + // Read the number of R + readLine(); + ss >> tmp_string >> tmp_string >> tmp_string >> tmp_string >> numberOfR; + + // Read the matrices + for (int i = 0; i < numberOfR; i++) + { + std::vector RCoord(3); + int nonZero; + + readLine(); + ss >> RCoord[0] >> RCoord[1] >> RCoord[2] >> nonZero; + RCoordinates.push_back(RCoord); + + std::vector csr_values(nonZero); + std::vector csr_col_ind(nonZero); + std::vector csr_row_ptr(matrixDimension + 1); + + // read CSR values + readLine(); + for (int i = 0; i < nonZero; i++) + { + ss >> csr_values[i]; + } + // read column indices + readLine(); + for (int i = 0; i < nonZero; i++) + { + ss >> csr_col_ind[i]; + } + // read row pointers + readLine(); + for (int i = 0; i < matrixDimension + 1; i++) + { + ss >> csr_row_ptr[i]; + } + + SparseMatrix matrix(matrixDimension, matrixDimension); + matrix.readCSR(csr_values, csr_col_ind, csr_row_ptr); + sparse_matrices.push_back(matrix); + } +} + +// function to get R coordinate +template +std::vector csrFileReader::getRCoordinate(int index) const +{ + if (index < 0 || index >= RCoordinates.size()) + { + ModuleBase::WARNING_QUIT("csrFileReader::getRCoordinate", "Index out of range"); + } + return RCoordinates[index]; +} + +// function to get matrix +template +SparseMatrix csrFileReader::getMatrix(int index) const +{ + if (index < 0 || index >= sparse_matrices.size()) + { + ModuleBase::WARNING_QUIT("csrFileReader::getMatrix", "Index out of range"); + } + return sparse_matrices[index]; +} + +// function to get matrix using R coordinate +template +SparseMatrix csrFileReader::getMatrix(int Rx, int Ry, int Rz) +{ + for (int i = 0; i < RCoordinates.size(); i++) + { + if (RCoordinates[i][0] == Rx && RCoordinates[i][1] == Ry && RCoordinates[i][2] == Rz) + { + return sparse_matrices[i]; + } + } + ModuleBase::WARNING_QUIT("csrFileReader::getMatrix", "R coordinate not found"); +} + +// function to get matrix +template +int csrFileReader::getNumberOfR() const +{ + return numberOfR; +} + +// function to get matrixDimension +template +int csrFileReader::getMatrixDimension() const +{ + return matrixDimension; +} + +// function to get step +template +int csrFileReader::getStep() const +{ + return step; +} + +// T of AtomPair can be double +template class csrFileReader; +// ToDo: T of AtomPair can be std::complex +// template class csrFileReader>; + +} // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/csr_reader.h b/source/module_io/csr_reader.h new file mode 100644 index 00000000000..15adf754279 --- /dev/null +++ b/source/module_io/csr_reader.h @@ -0,0 +1,71 @@ +#ifndef CSR_READER_H +#define CSR_READER_H + +#include + +#include "file_reader.h" +#include "sparse_matrix.h" + +namespace ModuleIO +{ + +/** + * @brief Class to read CSR file + * @details This class is used to read CSR file + * the default format is like: + * ``` + * STEP: 0 + * Matrix Dimension of S(R): 4 + * Matrix number of S(R): 2 + * 0 1 1 2 # (0,1,1) is the R coordinate, 2 is the number of non-zero elements + * 4.00e+00 7.00e+00 # non-zero elements + * 3 2 # column indices + * 0 1 2 2 2 # row pointer + * 0 0 0 3 # the second R coordinate and number of non-zero elements + * 5.00e+00 6.00e+00 1.00e+01 + * 2 3 3 + * 0 0 0 2 3 + * ``` + * It will store the R coordinates and sparse matrices as two vectors. + * One can use getter functions to get the R coordinates and sparse matrices, + * and related info including step, matrix dimension, number of R. + */ +template +class csrFileReader : public FileReader +{ + public: + // Constructor + csrFileReader(const std::string& filename); + + // read all matrices of all R coordinates + void parseFile(); + + // get number of R + int getNumberOfR() const; + + // get sparse matrix of a specific R coordinate + SparseMatrix getMatrix(int Rx, int Ry, int Rz); + + // get matrix by using index + SparseMatrix getMatrix(int index) const; + + // get R coordinate using index + std::vector getRCoordinate(int index) const; + + // get step + int getStep() const; + + // get matrix dimension + int getMatrixDimension() const; + + private: + std::vector> RCoordinates; + std::vector> sparse_matrices; + int step; + int matrixDimension; + int numberOfR; +}; + +} // namespace ModuleIO + +#endif // READ_CSR_H \ No newline at end of file diff --git a/source/module_io/file_reader.cpp b/source/module_io/file_reader.cpp new file mode 100644 index 00000000000..fd7a1e18b5d --- /dev/null +++ b/source/module_io/file_reader.cpp @@ -0,0 +1,50 @@ +#include "file_reader.h" + +#include "module_base/tool_quit.h" + +namespace ModuleIO +{ + +// Constructor +FileReader::FileReader(std::string filename) +{ + ifs.open(filename.c_str()); + if (!ifs.is_open()) + { + ModuleBase::WARNING_QUIT("FileReader::FileReader", "Error opening file"); + } +} + +// Destructor +FileReader::~FileReader() +{ + if (ifs.is_open()) + { + ifs.close(); + } +} + +// Function to check if file is open +bool FileReader::isOpen() const +{ + return ifs.is_open(); +} + +// Function to read a line and return string stream +void FileReader::readLine() +{ + // add warning if file is not open + if (!ifs.eof()) + { + std::string line; + std::getline(ifs, line); + ss.clear(); + ss.str(line); + } + else + { + ModuleBase::WARNING_QUIT("FileReader::readLine", "End of file"); + } +} + +} // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/file_reader.h b/source/module_io/file_reader.h new file mode 100644 index 00000000000..bd2c7481d72 --- /dev/null +++ b/source/module_io/file_reader.h @@ -0,0 +1,39 @@ +#ifndef FILE_READER_H +#define FILE_READER_H + +#include +#include +#include + +namespace ModuleIO +{ + +/** + * @brief A base class of file reader + * @details This class is supposed to be a base class to read a text file. + * it will open a file with a given filename. The function readLine() + * will read a line to a string stream. The function isOpen() check if the file + * is open. The destructor will close the file automatically. + */ +class FileReader +{ + public: + // Default constructor + FileReader(std::string filename); + ~FileReader(); + + // Check if file is open + bool isOpen() const; + + // read a line to string stream + void readLine(); + + std::stringstream ss; + + private: + std::ifstream ifs; +}; + +} // namespace ModuleIO + +#endif \ No newline at end of file diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index 54fd360095b..a43503c196b 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -120,4 +120,16 @@ AddTest( TARGET io_sparse_matrix_test LIBS base ${math_libs} device SOURCES sparse_matrix_test.cpp ../sparse_matrix.cpp +) + +AddTest( + TARGET io_file_reader_test + LIBS base ${math_libs} device + SOURCES file_reader_test.cpp ../file_reader.cpp +) + +AddTest( + TARGET io_csr_reader_test + LIBS base ${math_libs} device + SOURCES csr_reader_test.cpp ../csr_reader.cpp ../file_reader.cpp ../sparse_matrix.cpp ) \ No newline at end of file diff --git a/source/module_io/test/csr_reader_test.cpp b/source/module_io/test/csr_reader_test.cpp new file mode 100644 index 00000000000..ee2dcbd4f82 --- /dev/null +++ b/source/module_io/test/csr_reader_test.cpp @@ -0,0 +1,100 @@ +#include "module_io/csr_reader.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +/************************************************ + * unit test of csr_reader.cpp + ***********************************************/ + +/** + * - Tested Functions: + * - csrFileReader() + * - Constructor + * - parseFile() + * - Read the file and parse it + * - getNumberOfR() + * - Get number of R + * - getMatrix(rx, ry, rz) + * - Get sparse matrix of a specific R coordinate + * - getMatrix(index) + * - Get matrix by using index + * - getRCoordinate() + * - Get R coordinate using index + * - getStep() + * - Get step + * - getMatrixDimension() + * - Get matrix dimension + */ + +class csrFileReaderTest : public testing::Test +{ + protected: + std::string filename = "./support/SR.csr"; +}; + +TEST_F(csrFileReaderTest, CsrReader) +{ + ModuleIO::csrFileReader csr(filename); + // Check if file is open + EXPECT_TRUE(csr.isOpen()); + // get step + EXPECT_EQ(csr.getStep(), 0); + // get matrix dimension + EXPECT_EQ(csr.getMatrixDimension(), 4); + // get number of R + EXPECT_EQ(csr.getNumberOfR(), 2); + // get R coordinate using index + std::vector RCoord; + // get R coordinate using index + RCoord = csr.getRCoordinate(0); + EXPECT_EQ(RCoord[0], 0); + EXPECT_EQ(RCoord[1], 1); + EXPECT_EQ(RCoord[2], 1); + RCoord = csr.getRCoordinate(1); + EXPECT_EQ(RCoord[0], 0); + EXPECT_EQ(RCoord[1], 0); + EXPECT_EQ(RCoord[2], 0); + + // get matrix by using index + ModuleIO::SparseMatrix sparse_matrix; + ModuleIO::SparseMatrix sparse_matrix1; + // the first matrix should be + // 0 0 0 4 + // 0 0 7 0 + // 0 0 0 0 + // 0 0 0 0 + // the second matrix should be + // 0 0 0 0 + // 0 0 0 0 + // 0 0 5 6 + // 0 0 0 10 + sparse_matrix = csr.getMatrix(0); + sparse_matrix1 = csr.getMatrix(0, 1, 1); + for (const auto& element : sparse_matrix.getElements()) + { + auto it = sparse_matrix1.getElements().find(element.first); + EXPECT_EQ(it->first.first, element.first.first); + EXPECT_EQ(it->first.second, element.first.second); + EXPECT_DOUBLE_EQ(it->second, element.second); + //std::cout << "element( " << element.first.first << ", " << element.first.second << " ) = " << element.second << std::endl; + } + EXPECT_DOUBLE_EQ(sparse_matrix(0, 3), 4.0); + EXPECT_DOUBLE_EQ(sparse_matrix(1, 2), 7.0); + EXPECT_DOUBLE_EQ(sparse_matrix(0, 0), 0.0); + // the second R + sparse_matrix = csr.getMatrix(1); + sparse_matrix1 = csr.getMatrix(0, 0, 0); + for (const auto& element : sparse_matrix.getElements()) + { + auto it = sparse_matrix1.getElements().find(element.first); + EXPECT_EQ(it->first.first, element.first.first); + EXPECT_EQ(it->first.second, element.first.second); + EXPECT_DOUBLE_EQ(it->second, element.second); + //std::cout << "element( " << element.first.first << ", " << element.first.second << " ) = " << element.second << std::endl; + } + EXPECT_DOUBLE_EQ(sparse_matrix(2, 2), 5.0); + EXPECT_DOUBLE_EQ(sparse_matrix(2, 3), 6.0); + EXPECT_DOUBLE_EQ(sparse_matrix(3, 3), 10.0); + EXPECT_DOUBLE_EQ(sparse_matrix(0, 0), 0.0); +} \ No newline at end of file diff --git a/source/module_io/test/file_reader_test.cpp b/source/module_io/test/file_reader_test.cpp new file mode 100644 index 00000000000..e1d1e38e00e --- /dev/null +++ b/source/module_io/test/file_reader_test.cpp @@ -0,0 +1,73 @@ +#include "module_io/file_reader.h" + +#include + +#include "gmock/gmock.h" +#include "gtest/gtest.h" + +/************************************************ + * unit test of file_reader.cpp + ***********************************************/ + +/** + * - Tested Functions: + * - FileReader() + * - Constructor + * - isOpen() + * - Check if file is open + * - readLine() + * - Read a line to string stream + */ + +class FileReaderTest : public testing::Test +{ + protected: + std::ofstream ofs; + std::string filename = "test.txt"; + void SetUp() override + { + ofs.open(filename.c_str()); + // write some info to the file + ofs << "step 0" << std::endl; + ofs << "matrix dimension: 4" << std::endl; + ofs << "number of R: 2" << std::endl; + ofs << "R: 0 0 0 2" << std::endl; + ofs << "1 2" << std::endl; + ofs.close(); + } + void TearDown() override + { + remove("test.txt"); + } + std::string output; +}; + +TEST_F(FileReaderTest, Constructor) +{ + ModuleIO::FileReader fr(filename); + // Check if file is open + EXPECT_TRUE(fr.isOpen()); +} + +TEST_F(FileReaderTest, ReadLine) +{ + ModuleIO::FileReader fr(filename); + // Check if file is open + EXPECT_TRUE(fr.isOpen()); + // read a line to string stream + fr.readLine(); + EXPECT_EQ(fr.ss.str(), "step 0"); + fr.readLine(); + EXPECT_EQ(fr.ss.str(), "matrix dimension: 4"); + fr.readLine(); + EXPECT_EQ(fr.ss.str(), "number of R: 2"); + fr.readLine(); + EXPECT_EQ(fr.ss.str(), "R: 0 0 0 2"); + fr.readLine(); + EXPECT_EQ(fr.ss.str(), "1 2"); + fr.readLine(); + testing::internal::CaptureStdout(); + EXPECT_EXIT(fr.readLine(), ::testing::ExitedWithCode(0), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("End of file")); +} \ No newline at end of file diff --git a/source/module_io/test/support/SR.csr b/source/module_io/test/support/SR.csr new file mode 100644 index 00000000000..e86b151277d --- /dev/null +++ b/source/module_io/test/support/SR.csr @@ -0,0 +1,11 @@ +STEP: 0 +Matrix Dimension of S(R): 4 +Matrix number of S(R): 2 +0 1 1 2 + 4.00e+00 7.00e+00 + 3 2 + 0 1 2 2 2 +0 0 0 3 + 5.00e+00 6.00e+00 1.00e+01 + 2 3 3 + 0 0 0 2 3 From 5a8db92ac262909dced2f7cdbf863012236e8b89 Mon Sep 17 00:00:00 2001 From: root Date: Wed, 26 Jul 2023 05:53:26 +0000 Subject: [PATCH 4/5] try to fix error --- .../module_hcontainer/test/output_hcontainer_test.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp b/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp index f4d870c1fd4..c0798b5456c 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp @@ -86,16 +86,16 @@ TEST_F(OutputHContainerTest, Write) remove("output_hcontainer.log"); hamilt::HContainer HR(&ParaV); double correct_array[16] = {1, 2, 0, 4, 5, 0, 7, 0, 3, 0, 5, 6, 7, 8, 0, 10}; + double correct_array1[16] = {1, 2, 0, 4, 5, 0, 7, 0, 3, 0, 5, 6, 7, 8, 0, 10}; // correct_array represent a matrix of // 1 2 0 4 // 5 0 7 0 // 3 0 5 6 // 7 8 0 10 hamilt::AtomPair ap1(0, 1, 0, 1, 1, &ParaV, correct_array); - hamilt::AtomPair ap2(1, 1, 0, 0, 0, &ParaV, correct_array); + hamilt::AtomPair ap2(1, 1, 0, 0, 0, &ParaV, correct_array1); HR.insert_pair(ap1); HR.insert_pair(ap2); - /* for (int ir = 0; ir < HR.size_R_loop(); ++ir) { int rx, ry, rz; @@ -130,7 +130,6 @@ TEST_F(OutputHContainerTest, Write) } HR.unfix_R(); } - */ double sparse_threshold = 0.1; hamilt::Output_HContainer output_HR(&HR, &ParaV, ucell, std::cout, sparse_threshold, 2); // the first R From 97e1f29eab697c92a31a67dd7d7fb2bca4ccdc3c Mon Sep 17 00:00:00 2001 From: root Date: Wed, 26 Jul 2023 07:04:29 +0000 Subject: [PATCH 5/5] fix error --- .../module_hcontainer/test/CMakeLists.txt | 12 ++++++++++-- .../test/output_hcontainer_test.cpp | 13 ++++++++++--- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt b/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt index 3f801ee2973..4b32fca0da2 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt @@ -22,8 +22,16 @@ AddTest( ) AddTest( - TARGET test_output_hcontainer + TARGET test_hcontainer_output LIBS base ${math_libs} device - SOURCES output_hcontainer_test.cpp ../output_hcontainer.cpp ../../../module_basis/module_ao/parallel_2d.cpp ../../../module_basis/module_ao/parallel_orbitals.cpp ../../../module_io/sparse_matrix.cpp ../base_matrix.cpp ../hcontainer.cpp ../atom_pair.cpp tmp_mocks.cpp + SOURCES output_hcontainer_test.cpp + tmp_mocks.cpp + ../output_hcontainer.cpp + ../base_matrix.cpp + ../hcontainer.cpp + ../atom_pair.cpp + ../../../module_basis/module_ao/parallel_2d.cpp + ../../../module_basis/module_ao/parallel_orbitals.cpp + ../../../module_io/sparse_matrix.cpp ) endif() \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp b/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp index c0798b5456c..cda35ddcc73 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/test/output_hcontainer_test.cpp @@ -67,9 +67,16 @@ class OutputHContainerTest : public testing::Test TEST_F(OutputHContainerTest, Write) { Parallel_Orbitals ParaV; - ParaV.atom_begin_row.resize(2); - ParaV.atom_begin_col.resize(2); - for (int i = 0; i < 2; i++) + ParaV.atom_begin_row.resize(3); + ParaV.atom_begin_col.resize(3); + // Here is a bug, "nat" is 2 in this test, but + // the Parallel_Orbitals::get_col_size function + // will run "iat += 1" even if iat == nat - 1 + // resulting in wrong indexing of atom_begin_col + // so it is resized to 3 to avoid this bug temporarily + // @ 2023-07-26 + // this bug is expected to be fixed soon + for (int i = 0; i < 3; i++) { ParaV.atom_begin_row[i] = i * 2; ParaV.atom_begin_col[i] = i * 2;