8000 Power iteration: check norm of initial vector by peterrum · Pull Request #18178 · dealii/dealii · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Power iteration: check norm of initial vector #18178

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Power iteration: check norm of initial vector #18178

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 46 additions & 36 deletions include/deal.II/lac/precondition.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2454,44 +2454,54 @@ namespace internal
internal::set_initial_guess(temp_vector1);
data.constraints.set_zero(temp_vector1);

if (data.eigenvalue_algorithm == internal::EigenvalueAlgorithm::lanczos)
/// We want to run the power iteration starting with temp_vector1, but
/// for that it needs to be a nonzero vector. However, there are
// situations where it only contains zeroes (e.g., if all dofs are
// constrained). If that is the case, there is nothing to do.
if (temp_vector1.l2_norm() != 0.0)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it possible to use all_zero() here? We do require that for the vector space vector concept (which this template doesn't check for).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed, that should be cheaper.

{
// set a very strict tolerance to force at least two iterations
IterationNumberControl control(data.eig_cg_n_iterations,
1e-10,
false,
false);

dealii::SolverCG<VectorType> solver(control);
solver.connect_eigenvalues_slot(
[&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
eigenvalue_tracker.slot(eigenvalues);
});

solver.solve(*matrix_ptr,
solution_old,
temp_vector1,
*data.preconditioner);

info.cg_iterations = control.last_step();
}
else if (data.eigenvalue_algorithm ==
internal::EigenvalueAlgorithm::power_iteration)
{
(void)degree;

Assert(degree != numbers::invalid_unsigned_int,
ExcMessage("Cannot estimate the minimal eigenvalue with the "
"power iteration"));

eigenvalue_tracker.values.push_back(
internal::power_iteration(*matrix_ptr,
temp_vector1,
*data.preconditioner,
data.eig_cg_n_iterations));
if (data.eigenvalue_algorithm ==
internal::EigenvalueAlgorithm::lanczos)
{
// set a very strict tolerance to force at least two iterations
IterationNumberControl control(data.eig_cg_n_iterations,
1e-10,
false,
false);

dealii::SolverCG<VectorType> solver(control);
solver.connect_eigenvalues_slot(
[&eigenvalue_tracker](
const std::vector<double> &eigenvalues) {
eigenvalue_tracker.slot(eigenvalues);
});

solver.solve(*matrix_ptr,
solution_old,
temp_vector1,
*data.preconditioner);

info.cg_iterations = control.last_step();
}
else if (data.eigenvalue_algorithm ==
internal::EigenvalueAlgorithm::power_iteration)
{
(void)degree;

Assert(degree != numbers::invalid_unsigned_int,
ExcMessage(
"Cannot estimate the minimal eigenvalue with the "
"power iteration"));

eigenvalue_tracker.values.push_back(
internal::power_iteration(*matrix_ptr,
temp_vector1,
*data.preconditioner,
data.eig_cg_n_iterations));
}
else
DEAL_II_NOT_IMPLEMENTED();
}
else
DEAL_II_NOT_IMPLEMENTED();

// read the eigenvalues from the attached eigenvalue tracker
if (eigenvalue_tracker.values.empty())
Expand Down
64 changes: 64 additions & 0 deletions tests/lac/precondition_chebyshev_08.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------


// Test PreconditionChebyshev::estimate_eigenvalues() returns 1 as
// as minimum and maximum eigenvalues for fully constrained systems.

#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>

#include "../tests.h"

void
test(const internal::EigenvalueAlgorithm &eigenvalue_algorithm)
{
FullMatrix<double> matrix(2, 2);
m C09C atrix(0, 0) = 1.0;
matrix(1, 1) = 1.0;

typename PreconditionChebyshev<FullMatrix<double>,
Vector<double>,
PreconditionIdentity>::AdditionalData ad;

ad.preconditioner = std::make_shared<PreconditionIdentity>();
ad.constraints.constrain_dof_to_zero(0);
ad.constraints.constrain_dof_to_zero(1);

ad.eigenvalue_algorithm = eigenvalue_algorithm;

PreconditionChebyshev<FullMatrix<double>,
Vector<double>,
PreconditionIdentity>
precon;

precon.initialize(matrix, ad);

Vector<double> vec(2);
const auto ev = precon.estimate_eigenvalues(vec);

AssertDimension(ev.min_eigenvalue_estimate, 1.0);
AssertDimension(ev.max_eigenvalue_estimate, 1.0);
}

int
main()
{
initlog();

test(internal::EigenvalueAlgorithm::lanczos);
test(internal::EigenvalueAlgorithm::power_iteration);

deallog << "OK!" << std::endl;
}
2 changes: 2 additions & 0 deletions tests/lac/precondition_chebyshev_08.output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@

DEAL::OK!
Loading
0