Skip to content

Commit 586880f

Browse files
loumalouomegaloumalouomega
committed
Implement AuxiliaryMatrixMarket for MatrixMarket I/O and refactor TrilinosSpaceExperimental to use it
1 parent e054f34 commit 586880f

4 files changed

Lines changed: 387 additions & 108 deletions

File tree

‎applications/TrilinosApplication/CMakeLists.txt‎

Lines changed: 7 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -116,18 +116,13 @@ ENDIF(CMAKE_UNITY_BUILD MATCHES ON)
116116

117117
# Exclude from unity build
118118
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/linear_solver_trilinos.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
119-
# The experimental (Tpetra) files include MatrixMarket_Tpetra.hpp → mmio_Tpetra.h, which shares the
120-
# same MM_IO_H include guard as Kratos' mmio.h but declares mm_read_mtx_crd_size with size_t*
121-
# instead of int*. When merged into a unity unit with Epetra-based files that have already pulled
122-
# in Kratos' mmio.h, the guard prevents mmio_Tpetra.h from being re-included, leaving the wrong
123-
# signature in scope. Compiling the experimental files as separate translation units ensures the
124-
# correct include order and avoids free-function name conflicts with the Epetra equivalents.
125-
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/custom_python/add_trilinos_space_experimental_to_python.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
126-
if(${KRATOS_BUILD_TESTING} MATCHES ON)
127-
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/tests/cpp_tests/test_trilinos_space_experimental.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
128-
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/tests/cpp_tests/trilinos_cpp_test_experimental_utilities.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
129-
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/tests/cpp_tests/utilities/test_cpp_test_utilities_experimental.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
130-
endif()
119+
# auxiliary_matrix_market.cpp is the *only* translation unit that includes
120+
# MatrixMarket_Tpetra.hpp (which pulls in mmio_Tpetra.h). mmio_Tpetra.h reuses
121+
# the MM_IO_H include guard that Kratos' own mmio.h already claims, but with a
122+
# different signature for mm_read_mtx_crd_size (size_t* vs int*). Compiling
123+
# this file as its own translation unit (never merged into a unity batch) ensures
124+
# the correct header order with no guard collision.
125+
set_source_files_properties (${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/auxiliary_matrix_market.cpp PROPERTIES SKIP_UNITY_BUILD_INCLUSION TRUE)
131126

132127
# Changing the .dll suffix to .pyd
133128
if(${CMAKE_SYSTEM_NAME} MATCHES "Windows")

‎applications/TrilinosApplication/custom_utilities/auxiliary_matrix_market.cpp‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
// KRATOS _____ _ _ _
2+
// |_ _| __(_) (_)_ __ ___ ___
3+
// | || '__| | | | '_ \ / _ \/ __|
4+
// | || | | | | | | | | (_) \__
5+
// |_||_| |_|_|_|_| |_|\___/|___/ APPLICATION
6+
//
7+
// License: BSD License
8+
// Kratos default license: kratos/license.txt
9+
//
10+
// Main authors: Vicente Mataix Ferrandiz
11+
//
12+
13+
/**
14+
* @file auxiliary_matrix_market.cpp
15+
* @brief Implementation of AuxiliaryMatrixMarket — the *only* file that includes
16+
* MatrixMarket_Tpetra.hpp in the TrilinosApplication.
17+
*
18+
* Keep this file marked SKIP_UNITY_BUILD_INCLUSION in CMakeLists.txt so that
19+
* it is always compiled as its own, isolated translation unit. See
20+
* auxiliary_matrix_market.h for a full explanation of the include-guard conflict
21+
* that makes this isolation necessary.
22+
*/
23+
24+
// Project includes — must come before MatrixMarket_Tpetra.hpp so that the safe
25+
// Tpetra headers (pulled in by auxiliary_matrix_market.h) are processed first.
26+
#include "custom_utilities/auxiliary_matrix_market.h"
27+
#include "trilinos_application.h"
28+
29+
// External includes — the conflicting header, isolated here intentionally.
30+
// MatrixMarket_Tpetra.hpp includes mmio_Tpetra.h, which reuses the MM_IO_H
31+
// include guard. Compiling this file alone guarantees that the guard has not
32+
// already been claimed by Kratos' own mmio.h.
33+
#include <MatrixMarket_Tpetra.hpp>
34+
35+
namespace Kratos
36+
{
37+
38+
///@name AuxiliaryMatrixMarket Implementation
39+
///@{
40+
41+
template<class TMatrixType, class TVectorType>
42+
typename AuxiliaryMatrixMarket<TMatrixType, TVectorType>::MatrixPointerType
43+
AuxiliaryMatrixMarket<TMatrixType, TVectorType>::ReadMatrixMarket(
44+
const std::string& rFileName,
45+
CommunicatorType& rComm
46+
)
47+
{
48+
KRATOS_TRY
49+
50+
// Read into a regular (non-FE) CrsMatrix first — the Tpetra reader requires it.
51+
auto p_comm = Teuchos::rcp(&rComm, false);
52+
auto p_crs = Tpetra::MatrixMarket::Reader<CrsMatrixType>::readSparseFile(rFileName, p_comm);
53+
54+
const auto p_row_map = p_crs->getRowMap();
55+
const auto p_col_map = p_crs->getColMap();
56+
const LO num_local_rows = static_cast<LO>(p_crs->getNodeNumRows());
57+
58+
// Compute the per-row maximum entry count (FECrsGraph takes a scalar bound).
59+
std::size_t max_entries_per_row = 0;
60+
for (LO i = 0; i < num_local_rows; ++i) {
61+
max_entries_per_row = std::max(max_entries_per_row,
62+
static_cast<std::size_t>(p_crs->getNumEntriesInLocalRow(i)));
63+
}
64+
65+
// Build an FECrsGraph mirroring the sparsity pattern of the read matrix.
66+
auto p_graph = Teuchos::rcp(new GraphType(p_row_map, p_col_map, max_entries_per_row));
67+
p_graph->beginAssembly();
68+
for (LO i = 0; i < num_local_rows; ++i) {
69+
const GO global_row = p_row_map->getGlobalElement(i);
70+
typename CrsMatrixType::local_inds_host_view_type local_cols;
71+
typename CrsMatrixType::values_host_view_type vals;
72+
p_crs->getLocalRowView(i, local_cols, vals);
73+
if (local_cols.extent(0) > 0) {
74+
Teuchos::Array<GO> global_cols(local_cols.extent(0));
75+
for (std::size_t j = 0; j < local_cols.extent(0); ++j) {
76+
global_cols[j] = p_col_map->getGlobalElement(local_cols(j));
77+
}
78+
p_graph->insertGlobalIndices(global_row,
79+
Teuchos::ArrayView<const GO>(global_cols.data(), global_cols.size()));
80+
}
81+
}
82+
p_graph->endAssembly();
83+
84+
// Construct FECrsMatrix from the closed graph.
85+
auto p_matrix = Teuchos::rcp(new MatrixType(
86+
Teuchos::rcp_const_cast<const GraphType>(p_graph)));
87+
88+
// Normalise fill state, then populate values from the read CrsMatrix.
89+
if (p_matrix->isFillActive()) {
90+
p_matrix->fillComplete();
91+
}
92+
p_matrix->beginAssembly();
93+
p_matrix->setAllToScalar(static_cast<ST>(0));
94+
for (LO i = 0; i < num_local_rows; ++i) {
95+
const GO global_row = p_row_map->getGlobalElement(i);
96+
typename CrsMatrixType::local_inds_host_view_type local_cols;
97+
typename CrsMatrixType::values_host_view_type vals;
98+
p_crs->getLocalRowView(i, local_cols, vals);
99+
if (vals.extent(0) > 0) {
100+
Teuchos::Array<GO> global_cols(local_cols.extent(0));
101+
for (std::size_t j = 0; j < local_cols.extent(0); ++j) {
102+
global_cols[j] = p_col_map->getGlobalElement(local_cols(j));
103+
}
104+
p_matrix->sumIntoGlobalValues(global_row, static_cast<LO>(global_cols.size()),
105+
vals.data(), global_cols.data());
106+
}
107+
}
108+
p_matrix->endAssembly();
109+
110+
return p_matrix;
111+
112+
KRATOS_CATCH("")
113+
}
114+
115+
template<class TMatrixType, class TVectorType>
116+
typename AuxiliaryMatrixMarket<TMatrixType, TVectorType>::VectorPointerType
117+
AuxiliaryMatrixMarket<TMatrixType, TVectorType>::ReadMatrixMarketVector(
118+
const std::string& rFileName,
119+
CommunicatorPointerType pComm,
120+
const int N
121+
)
122+
{
123+
KRATOS_TRY
124+
125+
// Create a uniform contiguous map over all N global DOFs.
126+
MapPointerType p_map = Teuchos::rcp(new MapType(static_cast<GO>(N), 0, pComm));
127+
128+
// Read dense file into a plain (non-FE) MultiVector.
129+
auto p_mv = Tpetra::MatrixMarket::Reader<CrsMatrixType>::readDenseFile(rFileName, pComm, p_map);
130+
131+
// Create the FEMultiVector (whose copy constructor is deleted) and copy via update().
132+
VectorPointerType p_vector;
133+
if constexpr (std::is_same_v<VectorType, Tpetra::FEMultiVector<ST, LO, GO, NT>>) {
134+
p_vector = Teuchos::rcp(new VectorType(p_map, Teuchos::null, 1));
135+
} else {
136+
p_vector = Teuchos::rcp(new VectorType(p_map));
137+
}
138+
p_vector->update(static_cast<ST>(1), *p_mv, static_cast<ST>(0));
139+
140+
return p_vector;
141+
142+
KRATOS_CATCH("")
143+
}
144+
145+
template<class TMatrixType, class TVectorType>
146+
void AuxiliaryMatrixMarket<TMatrixType, TVectorType>::WriteMatrixMarketMatrix(
147+
const char* pFileName,
148+
const MatrixType& rA,
149+
const bool Symmetric
150+
)
151+
{
152+
KRATOS_TRY
153+
154+
// FECrsMatrix inherits from CrsMatrix; cast so Writer can accept it.
155+
auto p_crs = Teuchos::rcp_static_cast<const CrsMatrixType>(Teuchos::rcpFromRef(rA));
156+
Tpetra::MatrixMarket::Writer<CrsMatrixType>::writeSparseFile(std::string(pFileName), p_crs);
157+
158+
KRATOS_CATCH("")
159+
}
160+
161+
template<class TMatrixType, class TVectorType>
162+
void AuxiliaryMatrixMarket<TMatrixType, TVectorType>::WriteMatrixMarketVector(
163+
const char* pFileName,
164+
const VectorType& rV
165+
)
166+
{
167+
KRATOS_TRY
168+
169+
// FEMultiVector inherits from MultiVector; cast so Writer can accept it.
170+
using MultiVectorType = Tpetra::MultiVector<ST, LO, GO, NT>;
171+
auto p_mv = Teuchos::rcp_static_cast<const MultiVectorType>(Teuchos::rcpFromRef(rV));
172+
Tpetra::MatrixMarket::Writer<CrsMatrixType>::writeDenseFile(std::string(pFileName), p_mv);
173+
174+
KRATOS_CATCH("")
175+
}
176+
177+
///@}
178+
179+
// ---------------------------------------------------------------------------
180+
// Explicit instantiation — provides the compiled symbols for the only
181+
// concrete type pair used in the TrilinosApplication. The declaration-side
182+
// "extern template" in the header prevents the compiler from re-generating
183+
// these symbols in every translation unit that includes the header.
184+
// ---------------------------------------------------------------------------
185+
template class AuxiliaryMatrixMarket<Tpetra::FECrsMatrix<>, Tpetra::FEMultiVector<>>;
186+
187+
} // namespace Kratos

0 commit comments

Comments
 (0)