Dear all,
In the attached example, I deliberately assign NaN values to a
TrilinosWrappers::SparseMatrix to test how the Amesos_SuperluDist direct
solver handles such inputs. The code is run using deal.II version 9.4.0.
Unfortunately, the program results in a segmentation fault. Specifically,
no exception (TrilinosWrappers::SolverDirect::ExcTrilinosError) is thrown
when the matrix containing NaN values is passed to the solver's solve
function.
In contrast, when running a sequential program with SparseDirectUMFPACK,
the solver correctly handles the issue, and I am able to catch the
exception SparseDirectUMFPACK::ExcUMFPACKError.
Can this be considered a Trilinos bug? Has this issue been addressed in
newer versions of deal.II or Trilinos?
Clearly, assigning NaN values to a matrix does not make much sense.
However, in the optimized release version of my program, bad inputs can
occasionally lead to cases where the cell matrix contains NaN values. In
such situations, I rely on catching exceptions (e.g., from
SparseDirectUMFPACK) to handle the error gracefully.
Best,
Simon
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/8a68096d-c0ba-4edc-976c-17788993c18an%40googlegroups.com.
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/trilinos_solver.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/mpi.h>
#include <iostream>
#include <limits>
int main(int argc, char *argv[])
{
// Initialize MPI
dealii::Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
const MPI_Comm mpi_communicator = MPI_COMM_WORLD;
// Create a 1D triangulation with 2 cells
dealii::Triangulation<1> triangulation;
dealii::GridGenerator::hyper_cube(triangulation, 0, 1);
triangulation.refine_global(1); // Create 2 cells
// Finite element and DoF handler
const dealii::FE_Q<1> fe(1); // Linear elements
dealii::DoFHandler<1> dof_handler(triangulation);
dof_handler.distribute_dofs(fe);
// IndexSet for locally owned DoFs
dealii::IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
dealii::IndexSet locally_relevant_dofs;
dealii::DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
// Trilinos sparse matrix
dealii::TrilinosWrappers::SparseMatrix system_matrix;
dealii::DynamicSparsityPattern dsp(locally_owned_dofs);
dealii::DoFTools::make_sparsity_pattern(dof_handler, dsp, dealii::AffineConstraints{}, true);
system_matrix.reinit(locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator);
dealii::SparsityTools::distribute_sparsity_pattern(dsp,
locally_owned_dofs,
mpi_communicator,
locally_relevant_dofs);
// Local matrix
const unsigned int dofs_per_cell = fe.dofs_per_cell;
dealii::FullMatrix<double> local_matrix(dofs_per_cell, dofs_per_cell);
std::vector<dealii::types::global_dof_index> local_dof_indices(dofs_per_cell);
// Assign NaN to local matrix and distribute it to global system matrix
for (const auto &cell : dof_handler.active_cell_iterators())
{
cell->get_dof_indices(local_dof_indices);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
local_matrix(i, j) = std::numeric_limits<double>::quiet_NaN(); // Fill local matrix with NaN values
cell->distribute_local_to_global(local_matrix, system_matrix);
}
// Finalize the global matrix assembly
system_matrix.compress(dealii::VectorOperation::add);
// Right-hand-side vector and solution vector
dealii::TrilinosWrappers::MPI::Vector rhs(locally_owned_dofs, mpi_communicator);
dealii::TrilinosWrappers::MPI::Vector solution(locally_owned_dofs, mpi_communicator);
// Solver
dealii::TrilinosWrappers::SolverDirect::AdditionalData additionalData;
additionalData.solver_type = "Amesos_Superludist";
dealii::SolverControl solverControl;
dealii::TrilinosWrappers::SolverDirect directSolver(solverControl, additionalData);
// Attempt to solve the system
directSolver.solve(system_matrix, solution, rhs);
// Output the solution
solution.print(std::cout);
return 0;
}
