matrix_scaling.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2001 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_matrix_scaling_h
#define dealii_matrix_scaling_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/partitioner.h>
 
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
 
#ifdef DEAL_II_WITH_TRILINOS
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_vector.h>
#endif
 
#ifdef DEAL_II_WITH_PETSC
#  include <deal.II/lac/petsc_sparse_matrix.h>
#  include <deal.II/lac/petsc_vector.h>
#endif
 
DEAL_II_NAMESPACE_OPEN
 
class MatrixScaling : public EnableObserverPointer
{
public:
  enum class ConvergenceNormType
  {
    l1,
 
    l_infty
  };
 
  struct AdditionalData
  {
    enum class ScalingAlgorithm
    {
 
      sinkhorn_knopp,
 
 
      l1_linf_symmetry_preserving
    };
 
    struct SKParameters
    {
      enum class NormType
      {
        l1,
 
        l_infty
      };
 
      explicit SKParameters(const NormType     norm_type      = NormType::l1,
                            const unsigned int max_iterations = 20);
 
      unsigned int max_iterations;
 
      NormType norm_type;
    };
 
    struct l1linfParameters
    {
      explicit l1linfParameters(const unsigned int start_inf_norm_steps = 1,
                                const unsigned int l1_norm_steps        = 3,
                                const unsigned int end_inf_norm_steps   = 1);
 
      unsigned int start_inf_norm_steps;
 
      unsigned int l1_norm_steps;
 
      unsigned int end_inf_norm_steps;
    };
 
    explicit AdditionalData(
      const double           scaling_tolerance = 1e-5,
      const ScalingAlgorithm alg =
        ScalingAlgorithm::l1_linf_symmetry_preserving,
      const SKParameters     sk_params     = SKParameters(),
      const l1linfParameters l1linf_params = l1linfParameters());
 
 
    double scaling_tolerance;
 
    ScalingAlgorithm algorithm;
 
    SKParameters sinkhorn_knopp_parameters;
 
    l1linfParameters l1linf_parameters;
  };
 
  explicit MatrixScaling(const AdditionalData &control = AdditionalData());
 
  ~MatrixScaling() = default;
 
  template <class Matrix>
  bool
  find_scaling_and_scale_matrix(Matrix &matrix);
 
  template <class Matrix, class VectorType>
  bool
  find_scaling_and_scale_linear_system(Matrix &matrix, VectorType &rhs);
 
  template <class VectorType>
  void
  scale_system_solution(VectorType &sol) const;
 
  const Vector<double> &
  get_row_scaling() const;
 
  const Vector<double> &
  get_column_scaling() const;
 
private:
  AdditionalData control;
 
  Vector<double> row_scaling;
 
  Vector<double> column_scaling;
 
  IndexSet locally_owned_rows;
 
  IndexSet locally_owned_cols;
 
  IndexSet ghost_columns;
 
  std::vector<unsigned int> ghost_column_owners;
 
  Utilities::MPI::Partitioner partitioner;
 
  template <class Matrix>
  bool
  do_l1_scaling(Matrix &matrix, const unsigned int nsteps);
 
  template <class Matrix>
  bool
  do_linfty_scaling(Matrix &matrix, const unsigned int nsteps);
 
  template <class Matrix>
  bool
  do_sk_scaling(Matrix &matrix, const unsigned int nsteps);
 
  /*
   * Check convergence of the sequential scaling algorithms. Only on either row
   * or columns for SK scaling.
   */
  template <typename Number>
  bool
  check_convergence(const Vector<Number>      &row_col_norm,
                    const ConvergenceNormType &norm_type) const;
 
  /*
   * Check convergence of the sequential scaling algorithms. On both row
   * and columns for symmetry preserving scaling.
   */
  template <typename Number>
  bool
  check_convergence(const Vector<Number>      &row_norm,
                    const Vector<Number>      &col_norm,
                    const ConvergenceNormType &norm_type) const;
 
  /*
   * Check convergence of the parallel scaling algorithm. Only on either row
   * or columns for SK scaling.
   */
  bool
  check_convergence(const Vector<double>      &local_row_col_norm,
                    const ConvergenceNormType &norm_type,
                    const MPI_Comm             mpi_communicator) const;
 
  /*
   * Check convergence of the parallel scaling algorithm. Only on both row
   * and columns for symmetry preserving scaling.
   */
  bool
  check_convergence(const Vector<double>      &local_row_norm,
                    const Vector<double>      &local_col_norm,
                    const ConvergenceNormType &norm_type,
                    const MPI_Comm             mpi_communicator) const;
 
  /*
   * Fill the send_data map with the partial column norms that need to be sent
   * to other MPI processes. Prepare the local_col_norms vector. This is a
   * helper function for the distributed implementation of the scaling
   * algorithms.
   */
  void
  send_prepare_col_norms(
    const std::map<types::global_dof_index, double> &partial_column_norms,
    std::map<unsigned int,
             std::vector<std::pair<types::global_dof_index, double>>>
                   &send_data,
    Vector<double> &local_col_norms);
 
  /*
   * Fill the send_column_norms map with the updated column norms that need to
   * be sent to other MPI processes. This is a helper function for the
   * distributed implementation of the scaling algorithms.
   */
  void
  send_prepare_updated_col_norms(
    const std::map<unsigned int,
                   std::vector<std::pair<types::global_dof_index, double>>>
                         &received_data,
    const Vector<double> &local_col_norms,
    std::map<unsigned int,
             std::vector<std::pair<types::global_dof_index, double>>>
      &send_column_norms);
};
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // dealii_matrix_scaling_h