deal.II version 9.8.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
matrix_scaling.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2001 - 2026 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/memory_consumption.h>
14#include <deal.II/base/numbers.h>
15
16#include <deal.II/lac/block_sparse_matrix.h>
17#include <deal.II/lac/block_sparse_matrix_ez.h>
18#include <deal.II/lac/block_vector.h>
19#include <deal.II/lac/full_matrix.h>
20#include <deal.II/lac/matrix_scaling.h>
21#include <deal.II/lac/sparse_matrix.h>
22#include <deal.II/lac/sparse_matrix_ez.h>
23#include <deal.II/lac/trilinos_tpetra_to_trilinos_wrappers.h>
24#include <deal.II/lac/vector.h>
25
26#include <boost/serialization/utility.hpp>
27
28DEAL_II_NAMESPACE_OPEN
29
30namespace PETScWrappers
31{
32 namespace MPI
33 {
34 class SparseMatrix;
35 class Vector;
36 } // namespace MPI
37} // namespace PETScWrappers
38
39#include "lac/matrix_scaling.inst"
40
41// Check if the matrix type is sequential
42template <typename Matrix>
43constexpr bool
44is_sequential_matrix()
45{
46 return std::is_same_v<Matrix, ::SparseMatrix<double>> ||
47 std::is_same_v<Matrix, ::SparseMatrix<float>> ||
48 std::is_same_v<Matrix, ::FullMatrix<double>> ||
49 std::is_same_v<Matrix, ::FullMatrix<float>> ||
50 std::is_same_v<Matrix, ::BlockSparseMatrix<double>> ||
51 std::is_same_v<Matrix, ::BlockSparseMatrix<float>> ||
52 std::is_same_v<Matrix, ::SparseMatrixEZ<double>> ||
53 std::is_same_v<Matrix, ::SparseMatrixEZ<float>>;
54}
55
56
57
64
65
66
75
76
77
78MatrixScaling::AdditionalData::AdditionalData(
79 const double scaling_tolerance,
80 const ScalingAlgorithm alg,
81 const SKParameters sk_params,
82 const l1linfParameters l1linf_params)
83 : scaling_tolerance(scaling_tolerance)
84 , algorithm(alg)
85 , sinkhorn_knopp_parameters(sk_params)
86 , l1linf_parameters(l1linf_params)
87{}
88
89
90
96
97
98
99template <class Matrix>
100bool
101MatrixScaling::find_scaling_and_scale_matrix(Matrix &matrix)
102{
103 bool converged = false;
104
105 if constexpr (is_sequential_matrix<Matrix>())
106 {
107 const auto n_rows = matrix.m();
108 const auto n_cols = matrix.n();
109
110 row_scaling.reinit(n_rows);
111 column_scaling.reinit(n_cols);
112
113 row_scaling = 1.0;
114 column_scaling = 1.0;
115
116 switch (control.algorithm)
117 {
118 case MatrixScaling::AdditionalData::ScalingAlgorithm::sinkhorn_knopp:
119 {
120 converged =
121 do_sk_scaling(matrix,
122 control.sinkhorn_knopp_parameters.max_iterations);
123
124 break;
125 }
126
127 case MatrixScaling::AdditionalData::ScalingAlgorithm::
128 l1_linf_symmetry_preserving:
129 {
130 if (control.l1linf_parameters.start_inf_norm_steps > 0 &&
131 !converged)
132 converged = do_linfty_scaling(
133 matrix, control.l1linf_parameters.start_inf_norm_steps);
134 if (control.l1linf_parameters.l1_norm_steps > 0 && !converged)
135 converged =
136 do_l1_scaling(matrix,
137 control.l1linf_parameters.l1_norm_steps);
138 if (control.l1linf_parameters.end_inf_norm_steps > 0 &&
139 !converged)
140 converged = do_linfty_scaling(
141 matrix, control.l1linf_parameters.end_inf_norm_steps);
142
143 break;
144 }
145
146 default:
147 DEAL_II_ASSERT_UNREACHABLE();
148 }
149 }
150 else if constexpr (
151#ifdef DEAL_II_WITH_TRILINOS
152 std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix> ||
153#endif
154#ifdef DEAL_II_WITH_PETSC
155 std::is_same_v<Matrix, PETScWrappers::MPI::SparseMatrix> ||
156#endif
157 false)
158 {
159 locally_owned_rows = matrix.locally_owned_range_indices();
160
161 // If the matrix is square then use the same partitioning for
162 // columns and rows to easily scale a linear system, otherwise create a
163 // new balanced partitioning for the columns. We do this because
164 // distributed matrices do not really distribute columns, but only rows:
165 // we want to avoid saving the whole column scaling on each MPI rank
166 if (matrix.n() == matrix.m())
167 locally_owned_cols = locally_owned_rows;
168 else
169 locally_owned_cols =
170 Utilities::MPI::create_evenly_distributed_partitioning(
171 matrix.get_mpi_communicator(), matrix.n());
172
173 ghost_columns.clear();
174 ghost_columns.set_size(matrix.n());
175
176 row_scaling.reinit(locally_owned_rows.n_elements());
177 column_scaling.reinit(locally_owned_cols.n_elements());
178
179 row_scaling = 1.0;
180 column_scaling = 1.0;
181
182 // Identify ghost columns
183 for (const auto row : locally_owned_rows)
184 {
185 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
186 {
187 if (!locally_owned_cols.is_element(it->column()))
188 ghost_columns.add_index(it->column());
189 }
190 }
191 partitioner.reinit(locally_owned_cols,
192 ghost_columns,
193 matrix.get_mpi_communicator());
194
195 ghost_column_owners =
196 Utilities::MPI::compute_index_owner(locally_owned_cols,
197 ghost_columns,
198 matrix.get_mpi_communicator());
199
200 switch (control.algorithm)
201 {
202 case MatrixScaling::AdditionalData::ScalingAlgorithm::sinkhorn_knopp:
203 {
204 converged =
205 do_sk_scaling(matrix,
206 control.sinkhorn_knopp_parameters.max_iterations);
207
208 break;
209 }
210
211 case MatrixScaling::AdditionalData::ScalingAlgorithm::
212 l1_linf_symmetry_preserving:
213 {
214 if (control.l1linf_parameters.start_inf_norm_steps > 0 &&
215 !converged)
216 converged = do_linfty_scaling(
217 matrix, control.l1linf_parameters.start_inf_norm_steps);
218 if (control.l1linf_parameters.l1_norm_steps > 0 && !converged)
219 converged =
220 do_l1_scaling(matrix,
221 control.l1linf_parameters.l1_norm_steps);
222 if (control.l1linf_parameters.end_inf_norm_steps > 0 &&
223 !converged)
224 converged = do_linfty_scaling(
225 matrix, control.l1linf_parameters.end_inf_norm_steps);
226
227 break;
228 }
229
230 default:
231 DEAL_II_ASSERT_UNREACHABLE();
232 }
233 }
234 else
235 {
236 (void)matrix;
237 Assert(false, ExcNotImplemented());
238 }
239
240 return converged;
241}
242
243
244
245template <class Matrix, class VectorType>
246bool
247MatrixScaling::find_scaling_and_scale_linear_system(Matrix &matrix,
248 VectorType &rhs)
249{
250 AssertDimension(matrix.m(), rhs.size());
251 AssertDimension(matrix.m(), matrix.n());
252
253 bool converged = false;
254
255 converged = find_scaling_and_scale_matrix(matrix);
256
257 if constexpr (std::is_same_v<VectorType, ::Vector<double>> ||
258 std::is_same_v<VectorType, ::Vector<float>>)
259 rhs.scale(row_scaling);
260 else if constexpr (std::is_same_v<VectorType, ::BlockVector<double>> ||
261 std::is_same_v<VectorType, ::BlockVector<float>>)
262 {
263 for (unsigned int i = 0; i < rhs.size(); ++i)
264 rhs[i] *= row_scaling[i];
265 }
266 else if constexpr (
267#ifdef DEAL_II_WITH_TRILINOS
268 std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector> ||
269#endif
270#ifdef DEAL_II_WITH_PETSC
271 std::is_same_v<VectorType, PETScWrappers::MPI::Vector> ||
272#endif
273 false)
274 {
275 Assert(matrix.locally_owned_range_indices() ==
276 rhs.locally_owned_elements(),
277 ExcMessage("Matrix and vector must have the same partitioning"));
278 std::vector<double> local_updates(row_scaling.size());
279 for (auto i : locally_owned_rows)
280 {
281 auto local_idx = partitioner.global_to_local(i);
282 local_updates[local_idx] = rhs[i] * row_scaling[local_idx];
283 }
284 rhs.set(locally_owned_rows.get_index_vector(), local_updates);
285 rhs.compress(VectorOperation::insert);
286 }
287 else
288 Assert(false, ExcNotImplemented());
289
290 return converged;
291}
292
293
294
295template <class VectorType>
296void
297MatrixScaling::scale_system_solution(VectorType &sol) const
298{
299 if constexpr (std::is_same_v<VectorType, ::Vector<double>> ||
300 std::is_same_v<VectorType, ::Vector<float>>)
301 {
302 AssertDimension(sol.size(), column_scaling.size());
303 sol.scale(column_scaling);
304 }
305 else if constexpr (std::is_same_v<VectorType, ::BlockVector<double>> ||
306 std::is_same_v<VectorType, ::BlockVector<float>>)
307 {
308 AssertDimension(sol.size(), column_scaling.size());
309 for (unsigned int i = 0; i < sol.size(); ++i)
310 sol[i] *= column_scaling[i];
311 }
312 else if constexpr (
313#ifdef DEAL_II_WITH_TRILINOS
314 std::is_same_v<VectorType, TrilinosWrappers::MPI::Vector> ||
315#endif
316#ifdef DEAL_II_WITH_PETSC
317 std::is_same_v<VectorType, PETScWrappers::MPI::Vector> ||
318#endif
319 false)
320 {
321 Assert(locally_owned_cols == sol.locally_owned_elements(),
322 ExcMessage("Matrix and vector must have the same partitioning"));
323 std::vector<double> local_updates(column_scaling.size());
324 for (auto i : locally_owned_cols)
325 {
326 auto local_idx = partitioner.global_to_local(i);
327 local_updates[local_idx] = sol[i] * column_scaling[local_idx];
328 }
329 sol.set(locally_owned_cols.get_index_vector(), local_updates);
330 sol.compress(VectorOperation::insert);
331 }
332 else
333 Assert(false, ExcNotImplemented());
334}
335
336
337
338const Vector<double> &
339MatrixScaling::get_row_scaling() const
340{
341 return row_scaling;
342}
343
344
345
346const Vector<double> &
351
352
353
354template <typename Number>
355bool
356MatrixScaling::check_convergence(
357 const Vector<Number> &row_col_norm,
358 const MatrixScaling::ConvergenceNormType &norm_type) const
359{
360 Number convergence_norm = 0;
361 switch (norm_type)
362 {
363 case ConvergenceNormType::l1:
364 {
365 for (const auto &val : row_col_norm)
366 convergence_norm += std::abs(val - Number(1.0));
367
368 return convergence_norm < control.scaling_tolerance;
369 }
370
371 case ConvergenceNormType::l_infty:
372 {
373 for (const auto &val : row_col_norm)
374 convergence_norm =
375 std::max(convergence_norm, std::abs(val - Number(1.0)));
376
377 return convergence_norm < control.scaling_tolerance;
378 }
379
380 default:
381 {
382 DEAL_II_ASSERT_UNREACHABLE();
383 return false;
384 }
385 }
386}
387
388
389
390template <typename Number>
391bool
392MatrixScaling::check_convergence(
393 const Vector<Number> &row_norm,
394 const Vector<Number> &col_norm,
395 const MatrixScaling::ConvergenceNormType &norm_type) const
396{
397 Number convergence_row_norm = 0;
398 Number convergence_col_norm = 0;
399 switch (norm_type)
400 {
401 case ConvergenceNormType::l1:
402 {
403 for (const auto &val : row_norm)
404 convergence_row_norm += std::abs(val - Number(1.0));
405 for (const auto &val : col_norm)
406 convergence_col_norm += std::abs(val - Number(1.0));
407
408 return (convergence_row_norm < control.scaling_tolerance &&
409 convergence_col_norm < control.scaling_tolerance);
410 }
411
412 case ConvergenceNormType::l_infty:
413 {
414 for (const auto &val : row_norm)
415 convergence_row_norm =
416 std::max(convergence_row_norm, std::abs(val - Number(1.0)));
417 for (const auto &val : col_norm)
418 convergence_col_norm =
419 std::max(convergence_col_norm, std::abs(val - Number(1.0)));
420
421 return (convergence_row_norm < control.scaling_tolerance &&
422 convergence_col_norm < control.scaling_tolerance);
423 }
424
425 default:
426 {
427 DEAL_II_ASSERT_UNREACHABLE();
428 return false;
429 }
430 }
431}
432
433
434
435bool
436MatrixScaling::check_convergence(
437 const Vector<double> &local_row_col_norm,
438 const MatrixScaling::ConvergenceNormType &norm_type,
439 const MPI_Comm mpi_communicator) const
440{
441 double convergence_norm = 0;
442 bool local_not_converged;
443
444 switch (norm_type)
445 {
446 case ConvergenceNormType::l1:
447 {
448 for (const auto &val : local_row_col_norm)
449 convergence_norm += std::abs(val - 1.0);
450
451 local_not_converged = !(convergence_norm < control.scaling_tolerance);
452 break;
453 }
454
455 case ConvergenceNormType::l_infty:
456 {
457 for (const auto &val : local_row_col_norm)
458 convergence_norm = std::max(convergence_norm, std::abs(val - 1.0));
459
460 local_not_converged = !(convergence_norm < control.scaling_tolerance);
461 break;
462 }
463
464 default:
465 {
466 DEAL_II_ASSERT_UNREACHABLE();
467 return false;
468 }
469 }
470
471 const bool any_not_converged =
472 Utilities::MPI::logical_or(local_not_converged, mpi_communicator);
473
474 return !any_not_converged;
475}
476
477
478
479bool
480MatrixScaling::check_convergence(
481 const Vector<double> &local_row_norm,
482 const Vector<double> &local_col_norm,
483 const MatrixScaling::ConvergenceNormType &norm_type,
484 const MPI_Comm mpi_communicator) const
485{
486 double convergence_row_norm = 0;
487 double convergence_col_norm = 0;
488 bool local_not_converged;
489
490 switch (norm_type)
491 {
492 case ConvergenceNormType::l1:
493 {
494 for (const auto &val : local_row_norm)
495 convergence_row_norm += std::abs(val - 1.0);
496 for (const auto &val : local_col_norm)
497 convergence_col_norm += std::abs(val - 1.0);
498
499 local_not_converged =
500 !(convergence_row_norm < control.scaling_tolerance &&
501 convergence_col_norm < control.scaling_tolerance);
502 break;
503 }
504
505 case ConvergenceNormType::l_infty:
506 {
507 for (const auto &val : local_row_norm)
508 convergence_row_norm =
509 std::max(convergence_row_norm, std::abs(val - 1.0));
510 for (const auto &val : local_col_norm)
511 convergence_col_norm =
512 std::max(convergence_col_norm, std::abs(val - 1.0));
513
514 local_not_converged =
515 !(convergence_row_norm < control.scaling_tolerance &&
516 convergence_col_norm < control.scaling_tolerance);
517 break;
518 }
519
520 default:
521 {
522 DEAL_II_ASSERT_UNREACHABLE();
523 return false;
524 }
525 }
526
527 const bool any_not_converged =
528 Utilities::MPI::logical_or(local_not_converged, mpi_communicator);
529
530 return !any_not_converged;
531}
532
533
534
535void
536MatrixScaling::send_prepare_col_norms(
537 const std::map<types::global_dof_index, double> &partial_column_norms,
538 std::map<unsigned int,
539 std::vector<std::pair<types::global_dof_index, double>>> &send_data,
540 Vector<double> &local_col_norms)
541{
542 for (const auto &[col_idx, norm_value] : partial_column_norms)
543 {
544 if (locally_owned_cols.is_element(col_idx))
545 {
546 unsigned int local_idx = partitioner.global_to_local(col_idx);
547 local_col_norms[local_idx] = norm_value;
548 }
549 else
550 {
551 auto ghost_pos = ghost_columns.index_within_set(col_idx);
552 unsigned int target_rank = ghost_column_owners[ghost_pos];
553
554 send_data[target_rank].emplace_back(col_idx, norm_value);
555 }
556 }
557}
558
559
560
561void
562MatrixScaling::send_prepare_updated_col_norms(
563 const std::map<unsigned int,
564 std::vector<std::pair<types::global_dof_index, double>>>
565 &received_data,
566 const Vector<double> &local_col_norms,
567 std::map<unsigned int,
568 std::vector<std::pair<types::global_dof_index, double>>>
569 &send_column_norms)
570{
571 // For each rank that requested column norm data from me,
572 // send them the corresponding column norms updated
573 for (const auto &[sender_rank, pairs] : received_data)
574 {
575 for (const auto &[global_col, norm_contribution] : pairs)
576 {
577 // I received norm data for global_col, so sender_rank needs
578 // my scaling value for global_col
579 if (locally_owned_cols.is_element(global_col))
580 {
581 unsigned int local_idx = partitioner.global_to_local(global_col);
582 send_column_norms[sender_rank].emplace_back(
583 global_col, local_col_norms[local_idx]);
584 }
585 }
586 }
587}
588
589
590
591template <class Matrix>
592bool
593MatrixScaling::do_l1_scaling(Matrix &matrix, const unsigned int nsteps)
594{
595 if constexpr (is_sequential_matrix<Matrix>())
596 {
597 using Number = typename Matrix::value_type;
598
599 Vector<Number> row_norms(matrix.m());
600 Vector<Number> col_norms(matrix.n());
601
602 for (unsigned int i = 0; i < nsteps; i++)
603 {
604 row_norms = 0;
605 col_norms = 0;
606
607 for (unsigned int row = 0; row < matrix.m(); ++row)
608 {
609 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
610 {
611 row_norms[row] += std::abs(it->value());
612 col_norms[it->column()] += std::abs(it->value());
613 }
614 }
615
616 if (check_convergence(row_norms,
617 col_norms,
618 MatrixScaling::ConvergenceNormType::l1))
619 return true;
620
621 for (unsigned int row = 0; row < matrix.m(); ++row)
622 {
623 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
624 {
625 matrix.set(row,
626 it->column(),
627 it->value() / std::sqrt(row_norms[row] *
628 col_norms[it->column()]));
629 }
630 row_scaling[row] /= std::sqrt(row_norms[row]);
631 }
632 for (unsigned int col = 0; col < matrix.n(); ++col)
633 {
634 column_scaling[col] /= std::sqrt(col_norms[col]);
635 }
636 }
637 }
638 else if constexpr (
639#ifdef DEAL_II_WITH_TRILINOS
640 std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix> ||
641#endif
642#ifdef DEAL_II_WITH_PETSC
643 std::is_same_v<Matrix, PETScWrappers::MPI::SparseMatrix> ||
644#endif
645 false)
646 {
647 Vector<double> local_row_norms(locally_owned_rows.n_elements());
648 Vector<double> local_col_norms(locally_owned_cols.n_elements());
649 std::map<types::global_dof_index, double>
650 partial_column_norms; // column -> local norm
651
652 for (unsigned int i = 0; i < nsteps; i++)
653 {
654 local_row_norms = 0;
655 local_col_norms = 0;
656 partial_column_norms.clear();
657
658 for (const auto row : locally_owned_rows)
659 {
660 auto local_row_idx = partitioner.global_to_local(row);
661 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
662 {
663 local_row_norms[local_row_idx] += std::abs(it->value());
664 partial_column_norms[it->column()] += std::abs(it->value());
665 }
666 }
667 // Communicate partial column norms
668 std::map<unsigned int,
669 std::vector<std::pair<types::global_dof_index, double>>>
670 send_data;
671
672 send_prepare_col_norms(partial_column_norms,
673 send_data,
674 local_col_norms);
675
676 auto received_data =
677 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
678 send_data);
679
680 // Process received data and fill local_col_norms
681 for (const auto &[sender_rank, pairs] : received_data)
682 {
683 for (const auto &[global_col, contribution] : pairs)
684 {
685 unsigned int local_idx =
686 partitioner.global_to_local(global_col);
687 local_col_norms[local_idx] += contribution;
688 }
689 }
690
691 if (check_convergence(local_row_norms,
692 local_col_norms,
693 MatrixScaling::ConvergenceNormType::l1,
694 matrix.get_mpi_communicator()))
695 return true;
696
697
698 for (unsigned int i = 0; i < local_row_norms.size(); ++i)
699 row_scaling[i] /= std::sqrt(local_row_norms[i]);
700 for (unsigned int i = 0; i < local_col_norms.size(); ++i)
701 column_scaling[i] /= std::sqrt(local_col_norms[i]);
702
703 // Communicate column norm values to all ranks that need them
704 std::map<unsigned int,
705 std::vector<std::pair<types::global_dof_index, double>>>
706 send_column_norms;
707
708 send_prepare_updated_col_norms(received_data,
709 local_col_norms,
710 send_column_norms);
711
712 auto received_column_norms =
713 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
714 send_column_norms);
715
716 std::map<types::global_dof_index, double> ghost_column_norms_lookup;
717 for (const auto &[sender_rank, pairs] : received_column_norms)
718 {
719 for (const auto &[col_id, norm_val] : pairs)
720 {
721 ghost_column_norms_lookup[col_id] = norm_val;
722 }
723 }
724
725 // scale the matrix
726 for (const auto row : locally_owned_rows)
727 {
728 auto local_row_idx = partitioner.global_to_local(row);
729 auto row_size = matrix.row_length(row);
730 std::vector<double> values(row_size);
731 std::vector<types::global_dof_index> columns(row_size);
732 unsigned int idx = 0;
733 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
734 {
735 if (locally_owned_cols.is_element(it->column()))
736 {
737 unsigned int local_col_idx =
738 partitioner.global_to_local(it->column());
739 columns[idx] = it->column();
740 values[idx] =
741 it->value() / std::sqrt(local_col_norms[local_col_idx] *
742 local_row_norms[local_row_idx]);
743 ++idx;
744 }
745 else
746 {
747 double col_norms =
748 ghost_column_norms_lookup[it->column()];
749 columns[idx] = it->column();
750 values[idx] =
751 it->value() /
752 std::sqrt(col_norms * local_row_norms[local_row_idx]);
753 ++idx;
754 }
755 }
756 matrix.set(row, columns, values);
757 if constexpr (std::is_same_v<Matrix,
758 PETScWrappers::MPI::SparseMatrix>)
759 matrix.compress(VectorOperation::insert);
760 }
761#ifdef DEAL_II_WITH_TRILINOS
762 if constexpr (std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix>)
763 matrix.compress(VectorOperation::insert);
764#endif
765 }
766 }
767 return false;
768}
769
770
771
772template <class Matrix>
773bool
774MatrixScaling::do_linfty_scaling(Matrix &matrix, const unsigned int nsteps)
775{
776 if constexpr (is_sequential_matrix<Matrix>())
777 {
778 using Number = typename Matrix::value_type;
779
780 Vector<Number> row_norms(matrix.m());
781 Vector<Number> col_norms(matrix.n());
782
783 for (unsigned int i = 0; i < nsteps; i++)
784 {
785 row_norms = 0;
786 col_norms = 0;
787
788 for (unsigned int row = 0; row < matrix.m(); ++row)
789 {
790 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
791 {
792 row_norms[row] =
793 std::max(row_norms[row], std::abs(it->value()));
794 col_norms[it->column()] =
795 std::max(col_norms[it->column()], std::abs(it->value()));
796 }
797 }
798
799 if (check_convergence(row_norms,
800 col_norms,
801 MatrixScaling::ConvergenceNormType::l_infty))
802 return true;
803
804
805 for (unsigned int row = 0; row < matrix.m(); ++row)
806 {
807 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
808 {
809 matrix.set(row,
810 it->column(),
811 it->value() / std::sqrt(row_norms[row] *
812 col_norms[it->column()]));
813 }
814 row_scaling[row] /= std::sqrt(row_norms[row]);
815 }
816 for (unsigned int col = 0; col < matrix.n(); ++col)
817 {
818 column_scaling[col] /= std::sqrt(col_norms[col]);
819 }
820 }
821 }
822 else if constexpr (
823#ifdef DEAL_II_WITH_TRILINOS
824 std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix> ||
825#endif
826#ifdef DEAL_II_WITH_PETSC
827 std::is_same_v<Matrix, PETScWrappers::MPI::SparseMatrix> ||
828#endif
829 false)
830 {
831 Vector<double> local_row_norms(locally_owned_rows.n_elements());
832 Vector<double> local_col_norms(locally_owned_cols.n_elements());
833 std::map<types::global_dof_index, double>
834 partial_column_norms; // column -> local norm
835
836 for (unsigned int i = 0; i < nsteps; i++)
837 {
838 local_row_norms = 0;
839 local_col_norms = 0;
840 partial_column_norms.clear();
841
842 for (const auto row : locally_owned_rows)
843 {
844 auto local_row_idx = partitioner.global_to_local(row);
845 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
846 {
847 local_row_norms[local_row_idx] =
848 std::max(local_row_norms[local_row_idx],
849 std::abs(it->value()));
850 partial_column_norms[it->column()] =
851 std::max(partial_column_norms[it->column()],
852 std::abs(it->value()));
853 }
854 }
855 // Communicate partial column norms
856 std::map<unsigned int,
857 std::vector<std::pair<types::global_dof_index, double>>>
858 send_data;
859
860 send_prepare_col_norms(partial_column_norms,
861 send_data,
862 local_col_norms);
863
864 auto received_data =
865 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
866 send_data);
867
868 // Process received data and fill local_col_norms
869 for (const auto &[sender_rank, pairs] : received_data)
870 {
871 for (const auto &[global_col, contribution] : pairs)
872 {
873 unsigned int local_idx =
874 partitioner.global_to_local(global_col);
875 local_col_norms[local_idx] =
876 std::max(local_col_norms[local_idx], contribution);
877 }
878 }
879
880 if (check_convergence(local_row_norms,
881 local_col_norms,
882 MatrixScaling::ConvergenceNormType::l_infty,
883 matrix.get_mpi_communicator()))
884 return true;
885
886 for (unsigned int i = 0; i < local_row_norms.size(); ++i)
887 row_scaling[i] /= std::sqrt(local_row_norms[i]);
888 for (unsigned int i = 0; i < local_col_norms.size(); ++i)
889 column_scaling[i] /= std::sqrt(local_col_norms[i]);
890
891 // Communicate column norm values to all ranks that need them
892 std::map<unsigned int,
893 std::vector<std::pair<types::global_dof_index, double>>>
894 send_column_norms;
895
896 send_prepare_updated_col_norms(received_data,
897 local_col_norms,
898 send_column_norms);
899
900 auto received_column_norms =
901 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
902 send_column_norms);
903
904 std::map<types::global_dof_index, double> ghost_column_norms_lookup;
905 for (const auto &[sender_rank, pairs] : received_column_norms)
906 {
907 for (const auto &[col_id, norm_val] : pairs)
908 {
909 ghost_column_norms_lookup[col_id] = norm_val;
910 }
911 }
912
913 // scale the matrix
914 for (const auto row : locally_owned_rows)
915 {
916 auto local_row_idx = partitioner.global_to_local(row);
917 auto row_size = matrix.row_length(row);
918 std::vector<double> values(row_size);
919 std::vector<types::global_dof_index> columns(row_size);
920 unsigned int idx = 0;
921 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
922 {
923 if (locally_owned_cols.is_element(it->column()))
924 {
925 unsigned int local_col_idx =
926 partitioner.global_to_local(it->column());
927 columns[idx] = it->column();
928 values[idx] =
929 it->value() / std::sqrt(local_col_norms[local_col_idx] *
930 local_row_norms[local_row_idx]);
931 ++idx;
932 }
933 else
934 {
935 double col_norms =
936 ghost_column_norms_lookup[it->column()];
937 columns[idx] = it->column();
938 values[idx] =
939 it->value() /
940 std::sqrt(col_norms * local_row_norms[local_row_idx]);
941 ++idx;
942 }
943 }
944 matrix.set(row, columns, values);
945 if constexpr (std::is_same_v<Matrix,
946 PETScWrappers::MPI::SparseMatrix>)
947 matrix.compress(VectorOperation::insert);
948 }
949#ifdef DEAL_II_WITH_TRILINOS
950 if constexpr (std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix>)
951 matrix.compress(VectorOperation::insert);
952#endif
953 }
954 }
955 return false;
956}
957
958
959
960template <class Matrix>
961bool
962MatrixScaling::do_sk_scaling(Matrix &matrix, const unsigned int nsteps)
963{
964 if constexpr (is_sequential_matrix<Matrix>())
965 {
966 using Number = typename Matrix::value_type;
967
968 Vector<Number> row_norms(matrix.m());
969 Vector<Number> col_norms(matrix.n());
970
971 switch (control.sinkhorn_knopp_parameters.norm_type)
972 {
973 case MatrixScaling::AdditionalData::SKParameters::NormType::l1:
974 {
975 // Row_norms_0 to start the procedure
976 row_norms = 0;
977 for (unsigned int row = 0; row < matrix.m(); ++row)
978 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
979 row_norms[row] += std::abs(it->value());
980
981 for (unsigned int i = 0; i < nsteps; i++)
982 {
983 // Row step
984 col_norms = 0;
985 for (unsigned int row = 0; row < matrix.m(); ++row)
986 {
987 for (auto it = matrix.begin(row); it != matrix.end(row);
988 ++it)
989 {
990 matrix.set(row,
991 it->column(),
992 it->value() / row_norms[row]);
993 col_norms[it->column()] += std::abs(it->value());
994 }
995 row_scaling[row] /= row_norms[row];
996 }
997
998 if (check_convergence(col_norms,
999 MatrixScaling::ConvergenceNormType::l1))
1000 return true;
1001
1002 // Column step
1003 row_norms = 0;
1004 for (unsigned int row = 0; row < matrix.m(); ++row)
1005 {
1006 for (auto it = matrix.begin(row); it != matrix.end(row);
1007 ++it)
1008 {
1009 matrix.set(row,
1010 it->column(),
1011 it->value() / col_norms[it->column()]);
1012 row_norms[row] += std::abs(it->value());
1013 }
1014 }
1015 for (unsigned int col = 0; col < matrix.n(); ++col)
1016 {
1017 column_scaling[col] /= col_norms[col];
1018 }
1019
1020 if (check_convergence(row_norms,
1021 MatrixScaling::ConvergenceNormType::l1))
1022 return true;
1023 }
1024 break;
1025 }
1026
1027 case MatrixScaling::AdditionalData::SKParameters::NormType::l_infty:
1028 {
1029 // Row_norms_0 to start the procedure
1030 row_norms = 0;
1031 for (unsigned int row = 0; row < matrix.m(); ++row)
1032 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
1033 row_norms[row] +=
1034 std::max(row_norms[row], std::abs(it->value()));
1035
1036 for (unsigned int i = 0; i < nsteps; i++)
1037 {
1038 // Row step
1039 col_norms = 0;
1040 for (unsigned int row = 0; row < matrix.m(); ++row)
1041 {
1042 for (auto it = matrix.begin(row); it != matrix.end(row);
1043 ++it)
1044 {
1045 matrix.set(row,
1046 it->column(),
1047 it->value() / row_norms[row]);
1048 col_norms[it->column()] =
1049 std::max(col_norms[it->column()],
1050 std::abs(it->value()));
1051 }
1052 row_scaling[row] /= row_norms[row];
1053 }
1054
1055 if (check_convergence(
1056 col_norms, MatrixScaling::ConvergenceNormType::l_infty))
1057 return true;
1058
1059 // Column step
1060 row_norms = 0;
1061 for (unsigned int row = 0; row < matrix.m(); ++row)
1062 {
1063 for (auto it = matrix.begin(row); it != matrix.end(row);
1064 ++it)
1065 {
1066 matrix.set(row,
1067 it->column(),
1068 it->value() / col_norms[it->column()]);
1069 row_norms[row] +=
1070 std::max(row_norms[row], std::abs(it->value()));
1071 }
1072 }
1073 for (unsigned int col = 0; col < matrix.n(); ++col)
1074 {
1075 column_scaling[col] /= col_norms[col];
1076 }
1077
1078 if (check_convergence(
1079 row_norms, MatrixScaling::ConvergenceNormType::l_infty))
1080 return true;
1081 }
1082 break;
1083 }
1084
1085 default:
1086 DEAL_II_ASSERT_UNREACHABLE();
1087 }
1088 }
1089 else if constexpr (
1090#ifdef DEAL_II_WITH_TRILINOS
1091 std::is_same_v<Matrix, TrilinosWrappers::SparseMatrix> ||
1092#endif
1093#ifdef DEAL_II_WITH_PETSC
1094 std::is_same_v<Matrix, PETScWrappers::MPI::SparseMatrix> ||
1095#endif
1096 false)
1097 {
1098 Vector<double> local_row_norms(locally_owned_rows.n_elements());
1099 Vector<double> local_col_norms(locally_owned_cols.n_elements());
1100 std::map<types::global_dof_index, double>
1101 partial_column_norms; // column -> local norm
1102
1103 switch (control.sinkhorn_knopp_parameters.norm_type)
1104 {
1105 case MatrixScaling::AdditionalData::SKParameters::NormType::l1:
1106 {
1107 // Row_norms_0 to start the procedure
1108 local_row_norms = 0;
1109 for (const auto row : locally_owned_rows)
1110 {
1111 auto local_row_idx = partitioner.global_to_local(row);
1112 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
1113 local_row_norms[local_row_idx] += std::abs(it->value());
1114 }
1115
1116 for (unsigned int i = 0; i < nsteps; i++)
1117 {
1118 // Row step
1119 local_col_norms = 0;
1120 partial_column_norms.clear();
1121
1122 for (const auto row : locally_owned_rows)
1123 {
1124 auto local_row_idx = partitioner.global_to_local(row);
1125 auto row_size = matrix.row_length(row);
1126 std::vector<double> values(row_size);
1127 std::vector<types::global_dof_index> columns(row_size);
1128 unsigned int idx = 0;
1129 for (auto it = matrix.begin(row); it != matrix.end(row);
1130 ++it)
1131 {
1132 partial_column_norms[it->column()] += std::abs(
1133 it->value() / local_row_norms[local_row_idx]);
1134 columns[idx] = it->column();
1135 values[idx] =
1136 it->value() / local_row_norms[local_row_idx];
1137 ++idx;
1138 }
1139 matrix.set(row, columns, values);
1140 if constexpr (std::is_same_v<
1141 Matrix,
1142 PETScWrappers::MPI::SparseMatrix>)
1143 matrix.compress(VectorOperation::insert);
1144
1145 row_scaling[local_row_idx] /=
1146 local_row_norms[local_row_idx];
1147 }
1148#ifdef DEAL_II_WITH_TRILINOS
1149 if constexpr (std::is_same_v<Matrix,
1150 TrilinosWrappers::SparseMatrix>)
1151 matrix.compress(VectorOperation::insert);
1152#endif
1153
1154 // Communicate partial column norms
1155 std::map<
1156 unsigned int,
1157 std::vector<std::pair<types::global_dof_index, double>>>
1158 send_data;
1159
1160 send_prepare_col_norms(partial_column_norms,
1161 send_data,
1162 local_col_norms);
1163
1164 auto received_data =
1165 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
1166 send_data);
1167
1168 // Process received data and fill local_col_norms
1169 for (const auto &[sender_rank, pairs] : received_data)
1170 {
1171 for (const auto &[global_col, contribution] : pairs)
1172 {
1173 unsigned int local_idx =
1174 partitioner.global_to_local(global_col);
1175 local_col_norms[local_idx] += contribution;
1176 }
1177 }
1178
1179 // Convergence check only on columns
1180 if (check_convergence(local_col_norms,
1181 MatrixScaling::ConvergenceNormType::l1,
1182 matrix.get_mpi_communicator()))
1183 return true;
1184
1185 // Column step
1186 local_row_norms = 0;
1187
1188 // Communicate column norm values to all ranks that need them
1189 std::map<
1190 unsigned int,
1191 std::vector<std::pair<types::global_dof_index, double>>>
1192 send_column_norms;
1193
1194 send_prepare_updated_col_norms(received_data,
1195 local_col_norms,
1196 send_column_norms);
1197
1198 auto received_column_norms =
1199 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
1200 send_column_norms);
1201
1202 std::map<types::global_dof_index, double>
1203 ghost_column_norms_lookup;
1204 for (const auto &[sender_rank, pairs] : received_column_norms)
1205 {
1206 for (const auto &[col_id, norm_val] : pairs)
1207 {
1208 ghost_column_norms_lookup[col_id] = norm_val;
1209 }
1210 }
1211
1212 for (const auto row : locally_owned_rows)
1213 {
1214 auto local_row_idx = partitioner.global_to_local(row);
1215 auto row_size = matrix.row_length(row);
1216 std::vector<double> values(row_size);
1217 std::vector<types::global_dof_index> columns(row_size);
1218 unsigned int idx = 0;
1219 for (auto it = matrix.begin(row); it != matrix.end(row);
1220 ++it)
1221 {
1222 if (locally_owned_cols.is_element(it->column()))
1223 {
1224 unsigned int local_col_idx =
1225 partitioner.global_to_local(it->column());
1226
1227 local_row_norms[local_row_idx] += std::abs(
1228 it->value() / local_col_norms[local_col_idx]);
1229
1230 columns[idx] = it->column();
1231 values[idx] =
1232 it->value() / local_col_norms[local_col_idx];
1233 ++idx;
1234 }
1235 else
1236 {
1237 double col_norms =
1238 ghost_column_norms_lookup[it->column()];
1239
1240 local_row_norms[local_row_idx] +=
1241 std::abs(it->value() / col_norms);
1242
1243 columns[idx] = it->column();
1244 values[idx] = it->value() / col_norms;
1245 ++idx;
1246 }
1247 }
1248 matrix.set(row, columns, values);
1249 if constexpr (std::is_same_v<
1250 Matrix,
1251 PETScWrappers::MPI::SparseMatrix>)
1252 matrix.compress(VectorOperation::insert);
1253 }
1254#ifdef DEAL_II_WITH_TRILINOS
1255 if constexpr (std::is_same_v<Matrix,
1256 TrilinosWrappers::SparseMatrix>)
1257 matrix.compress(VectorOperation::insert);
1258#endif
1259 for (unsigned int i = 0; i < local_col_norms.size(); ++i)
1260 column_scaling[i] /= local_col_norms[i];
1261
1262 if (check_convergence(local_row_norms,
1263 MatrixScaling::ConvergenceNormType::l1,
1264 matrix.get_mpi_communicator()))
1265 return true;
1266 }
1267 break;
1268 }
1269
1270 case MatrixScaling::AdditionalData::SKParameters::NormType::l_infty:
1271 {
1272 // Row_norms_0 to start the procedure
1273 local_row_norms = 0;
1274 for (const auto row : locally_owned_rows)
1275 {
1276 auto local_row_idx = partitioner.global_to_local(row);
1277 for (auto it = matrix.begin(row); it != matrix.end(row); ++it)
1278 local_row_norms[local_row_idx] =
1279 std::max(local_row_norms[local_row_idx],
1280 std::abs(it->value()));
1281 }
1282
1283 for (unsigned int i = 0; i < nsteps; i++)
1284 {
1285 // Row step
1286 local_col_norms = 0;
1287 partial_column_norms.clear();
1288
1289 for (const auto row : locally_owned_rows)
1290 {
1291 auto local_row_idx = partitioner.global_to_local(row);
1292 auto row_size = matrix.row_length(row);
1293 std::vector<double> values(row_size);
1294 std::vector<types::global_dof_index> columns(row_size);
1295 unsigned int idx = 0;
1296 for (auto it = matrix.begin(row); it != matrix.end(row);
1297 ++it)
1298 {
1299 partial_column_norms[it->column()] =
1300 std::max(partial_column_norms[it->column()],
1301 std::abs(it->value() /
1302 local_row_norms[local_row_idx]));
1303 columns[idx] = it->column();
1304 values[idx] =
1305 it->value() / local_row_norms[local_row_idx];
1306 ++idx;
1307 }
1308 matrix.set(row, columns, values);
1309 if constexpr (std::is_same_v<
1310 Matrix,
1311 PETScWrappers::MPI::SparseMatrix>)
1312 matrix.compress(VectorOperation::insert);
1313
1314 row_scaling[local_row_idx] /=
1315 local_row_norms[local_row_idx];
1316 }
1317#ifdef DEAL_II_WITH_TRILINOS
1318 if constexpr (std::is_same_v<Matrix,
1319 TrilinosWrappers::SparseMatrix>)
1320 matrix.compress(VectorOperation::insert);
1321#endif
1322
1323 // Communicate partial column norms
1324 std::map<
1325 unsigned int,
1326 std::vector<std::pair<types::global_dof_index, double>>>
1327 send_data;
1328
1329 send_prepare_col_norms(partial_column_norms,
1330 send_data,
1331 local_col_norms);
1332
1333 auto received_data =
1334 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
1335 send_data);
1336
1337 // Process received data and fill local_col_norms
1338 for (const auto &[sender_rank, pairs] : received_data)
1339 {
1340 for (const auto &[global_col, contribution] : pairs)
1341 {
1342 unsigned int local_idx =
1343 partitioner.global_to_local(global_col);
1344 local_col_norms[local_idx] =
1345 std::max(local_col_norms[local_idx], contribution);
1346 }
1347 }
1348
1349 if (check_convergence(
1350 local_col_norms,
1351 MatrixScaling::ConvergenceNormType::l_infty,
1352 matrix.get_mpi_communicator()))
1353 return true;
1354
1355 // Column step
1356 local_row_norms = 0;
1357
1358 // Communicate column norm values to all ranks that need them
1359 std::map<
1360 unsigned int,
1361 std::vector<std::pair<types::global_dof_index, double>>>
1362 send_column_norms;
1363
1364 send_prepare_updated_col_norms(received_data,
1365 local_col_norms,
1366 send_column_norms);
1367
1368 auto received_column_norms =
1369 Utilities::MPI::some_to_some(matrix.get_mpi_communicator(),
1370 send_column_norms);
1371
1372 std::map<types::global_dof_index, double>
1373 ghost_column_norms_lookup;
1374 for (const auto &[sender_rank, pairs] : received_column_norms)
1375 {
1376 for (const auto &[col_id, norm_val] : pairs)
1377 {
1378 ghost_column_norms_lookup[col_id] = norm_val;
1379 }
1380 }
1381
1382 for (const auto row : locally_owned_rows)
1383 {
1384 auto local_row_idx = partitioner.global_to_local(row);
1385 auto row_size = matrix.row_length(row);
1386 std::vector<double> values(row_size);
1387 std::vector<types::global_dof_index> columns(row_size);
1388 unsigned int idx = 0;
1389 for (auto it = matrix.begin(row); it != matrix.end(row);
1390 ++it)
1391 {
1392 if (locally_owned_cols.is_element(it->column()))
1393 {
1394 unsigned int local_col_idx =
1395 partitioner.global_to_local(it->column());
1396
1397 local_row_norms[local_row_idx] =
1398 std::max(local_row_norms[local_row_idx],
1399 std::abs(
1400 it->value() /
1401 local_col_norms[local_col_idx]));
1402
1403 columns[idx] = it->column();
1404 values[idx] =
1405 it->value() / local_col_norms[local_col_idx];
1406 ++idx;
1407 }
1408 else
1409 {
1410 double col_norms =
1411 ghost_column_norms_lookup[it->column()];
1412
1413 local_row_norms[local_row_idx] =
1414 std::max(local_row_norms[local_row_idx],
1415 std::abs(it->value() / col_norms));
1416
1417 columns[idx] = it->column();
1418 values[idx] = it->value() / col_norms;
1419 ++idx;
1420 }
1421 }
1422 matrix.set(row, columns, values);
1423 if constexpr (std::is_same_v<
1424 Matrix,
1425 PETScWrappers::MPI::SparseMatrix>)
1426 matrix.compress(VectorOperation::insert);
1427 }
1428#ifdef DEAL_II_WITH_TRILINOS
1429 if constexpr (std::is_same_v<Matrix,
1430 TrilinosWrappers::SparseMatrix>)
1431 matrix.compress(VectorOperation::insert);
1432#endif
1433 for (unsigned int i = 0; i < local_col_norms.size(); ++i)
1434 column_scaling[i] /= local_col_norms[i];
1435
1436 if (check_convergence(
1437 local_row_norms,
1438 MatrixScaling::ConvergenceNormType::l_infty,
1439 matrix.get_mpi_communicator()))
1440 return true;
1441 }
1442 break;
1443 }
1444
1445 default:
1446 DEAL_II_ASSERT_UNREACHABLE();
1447 }
1448 }
1449 return false;
1450}
1451
1452
1453
1454#define InstantiateMatrixScaling(MATRIX, VECTOR) \
1455 template bool MatrixScaling::find_scaling_and_scale_matrix(MATRIX &); \
1456 template bool MatrixScaling::find_scaling_and_scale_linear_system(MATRIX &, \
1457 VECTOR &); \
1458 template void MatrixScaling::scale_system_solution(VECTOR &) const;
1459
1460#ifdef DEAL_II_TRILINOS_WITH_EPETRA
1461InstantiateMatrixScaling(TrilinosWrappers::SparseMatrix,
1462 TrilinosWrappers::MPI::Vector)
1463#endif
1464
1465#ifdef DEAL_II_WITH_PETSC
1466 InstantiateMatrixScaling(PETScWrappers::MPI::SparseMatrix,
1467 PETScWrappers::MPI::Vector)
1468#endif
1469
1470
1471 DEAL_II_NAMESPACE_CLOSE
block_sparse_matrix.h
block_sparse_matrix_ez.h
block_vector.h
MPI_Comm
MatrixScaling::partitioner
Utilities::MPI::Partitioner partitioner
Definition matrix_scaling.h:319
MatrixScaling::MatrixScaling
MatrixScaling(const AdditionalData &control=AdditionalData())
Definition matrix_scaling.cc:91
MatrixScaling::get_row_scaling
const Vector< double > & get_row_scaling() const
Definition matrix_scaling.cc:339
MatrixScaling::column_scaling
Vector< double > column_scaling
Definition matrix_scaling.h:292
MatrixScaling::do_linfty_scaling
bool do_linfty_scaling(Matrix &matrix, const unsigned int nsteps)
Definition matrix_scaling.cc:774
MatrixScaling::do_sk_scaling
bool do_sk_scaling(Matrix &matrix, const unsigned int nsteps)
Definition matrix_scaling.cc:962
MatrixScaling::locally_owned_cols
IndexSet locally_owned_cols
Definition matrix_scaling.h:302
MatrixScaling::control
AdditionalData control
Definition matrix_scaling.h:280
MatrixScaling::row_scaling
Vector< double > row_scaling
Definition matrix_scaling.h:286
MatrixScaling::find_scaling_and_scale_linear_system
bool find_scaling_and_scale_linear_system(Matrix &matrix, VectorType &rhs)
Definition matrix_scaling.cc:247
MatrixScaling::check_convergence
bool check_convergence(const Vector< Number > &row_col_norm, const ConvergenceNormType &norm_type) const
Definition matrix_scaling.cc:356
MatrixScaling::get_column_scaling
const Vector< double > & get_column_scaling() const
Definition matrix_scaling.cc:347
MatrixScaling::send_prepare_updated_col_norms
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)
Definition matrix_scaling.cc:562
MatrixScaling::ConvergenceNormType::l_infty
@ l_infty
l_infinity vector norm
Definition matrix_scaling.h:75
MatrixScaling::send_prepare_col_norms
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)
Definition matrix_scaling.cc:536
MatrixScaling::ghost_columns
IndexSet ghost_columns
Definition matrix_scaling.h:309
MatrixScaling::locally_owned_rows
IndexSet locally_owned_rows
Definition matrix_scaling.h:297
MatrixScaling::ghost_column_owners
std::vector< unsigned int > ghost_column_owners
Definition matrix_scaling.h:314
MatrixScaling::do_l1_scaling
bool do_l1_scaling(Matrix &matrix, const unsigned int nsteps)
Definition matrix_scaling.cc:593
MatrixScaling::scale_system_solution
void scale_system_solution(VectorType &sol) const
Definition matrix_scaling.cc:297
MatrixScaling::find_scaling_and_scale_matrix
bool find_scaling_and_scale_matrix(Matrix &matrix)
Definition matrix_scaling.cc:101
Vector::size
virtual size_type size() const override
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
DEAL_II_ASSERT_UNREACHABLE
#define DEAL_II_ASSERT_UNREACHABLE()
Definition exception_macros.h:781
full_matrix.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
is_sequential_matrix
constexpr bool is_sequential_matrix()
Definition matrix_scaling.cc:44
InstantiateMatrixScaling
#define InstantiateMatrixScaling(MATRIX, VECTOR)
Definition matrix_scaling.cc:1454
matrix_scaling.h
memory_consumption.h
Utilities::MPI::logical_or
T logical_or(const T &t, const MPI_Comm mpi_communicator)
Utilities::MPI::some_to_some
std::map< unsigned int, T > some_to_some(const MPI_Comm comm, const std::map< unsigned int, T > &objects_to_send)
Utilities::MPI::compute_index_owner
std::vector< unsigned int > compute_index_owner(const IndexSet &owned_indices, const IndexSet &indices_to_look_up, const MPI_Comm comm)
Definition mpi.cc:1820
Utilities::MPI::create_evenly_distributed_partitioning
IndexSet create_evenly_distributed_partitioning(const MPI_Comm comm, const types::global_dof_index total_size)
Definition mpi.cc:204
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6741
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6667
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6725
sparse_matrix.h
sparse_matrix_ez.h
MatrixScaling::AdditionalData::SKParameters::SKParameters
SKParameters(const NormType norm_type=NormType::l1, const unsigned int max_iterations=20)
Definition matrix_scaling.cc:58
MatrixScaling::AdditionalData::l1linfParameters::l1linfParameters
l1linfParameters(const unsigned int start_inf_norm_steps=1, const unsigned int l1_norm_steps=3, const unsigned int end_inf_norm_steps=1)
Definition matrix_scaling.cc:67
MatrixScaling::AdditionalData::ScalingAlgorithm::sinkhorn_knopp
@ sinkhorn_knopp
Sinkhorn-Knopp scaling algorithm.
Definition matrix_scaling.h:96
MatrixScaling::AdditionalData::AdditionalData
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())
Definition matrix_scaling.cc:78
trilinos_tpetra_to_trilinos_wrappers.h
Generated by doxygen 1.16.1