Commit 07018cfd authored by webmanue's avatar webmanue
Browse files

Merge branch '89-add-eigenpair-extraction' into 'master'

Resolve "Add eigenpair extraction"

Closes #89

See merge request mechanics-and-materials/ae108!78
parents a7eda3c9 542c8fc3
Pipeline #139034 passed with stages
in 6 minutes and 50 seconds
......@@ -64,6 +64,7 @@ add_library(${PROJECT_NAME}-cpppetsc
./nestMatrices.cc
./nestVectors.cc
./readFromViewer.cc
./scalarProduct.cc
./setName.cc
./vertexDataOffsets.cc
./writeToViewer.cc
......
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#pragma once
#include "ae108/cpppetsc/ParallelComputePolicy_fwd.h"
#include "ae108/cpppetsc/SequentialComputePolicy_fwd.h"
#include "ae108/cpppetsc/TaggedVector.h"
#include "ae108/cpppetsc/Vector.h"
namespace ae108 {
namespace cpppetsc {
/**
* @brief Computes the scalar product of the two vectors `x` and `y`.
*/
template <class Policy, class Tag>
typename Vector<Policy>::value_type
scalarProduct(const TaggedEntity<Vector<Policy>, Tag> &x,
const TaggedEntity<Vector<Policy>, Tag> &y) {
auto result = typename Vector<Policy>::value_type();
Policy::handleError(VecDot(x.unwrap().data(), y.unwrap().data(), &result));
return result;
}
} // namespace cpppetsc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cpppetsc/scalarProduct.h"
#include "ae108/cpppetsc/ParallelComputePolicy.h"
#include "ae108/cpppetsc/SequentialComputePolicy.h"
namespace ae108 {
namespace cpppetsc {
template typename Vector<SequentialComputePolicy>::value_type
scalarProduct(const TaggedEntity<Vector<SequentialComputePolicy>, LocalTag> &,
const TaggedEntity<Vector<SequentialComputePolicy>, LocalTag> &);
template typename Vector<SequentialComputePolicy>::value_type scalarProduct(
const TaggedEntity<Vector<SequentialComputePolicy>, DistributedTag> &,
const TaggedEntity<Vector<SequentialComputePolicy>, DistributedTag> &);
template typename Vector<SequentialComputePolicy>::value_type
scalarProduct(const TaggedEntity<Vector<SequentialComputePolicy>, GlobalTag> &,
const TaggedEntity<Vector<SequentialComputePolicy>, GlobalTag> &);
template typename Vector<ParallelComputePolicy>::value_type
scalarProduct(const TaggedEntity<Vector<ParallelComputePolicy>, LocalTag> &,
const TaggedEntity<Vector<ParallelComputePolicy>, LocalTag> &);
template typename Vector<ParallelComputePolicy>::value_type scalarProduct(
const TaggedEntity<Vector<ParallelComputePolicy>, DistributedTag> &,
const TaggedEntity<Vector<ParallelComputePolicy>, DistributedTag> &);
template typename Vector<ParallelComputePolicy>::value_type
scalarProduct(const TaggedEntity<Vector<ParallelComputePolicy>, GlobalTag> &,
const TaggedEntity<Vector<ParallelComputePolicy>, GlobalTag> &);
} // namespace cpppetsc
} // namespace ae108
\ No newline at end of file
......@@ -46,6 +46,7 @@ add_executable(${PROJECT_NAME}-CppPetscTests
./multiply_Test.cc
./nestMatrices_Test.cc
./nestVectors_Test.cc
./scalarProduct_Test.cc
./setName_Test.cc
./vertexDataOffsets_Test.cc
)
......
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cpppetsc/ParallelComputePolicy.h"
#include "ae108/cpppetsc/SequentialComputePolicy.h"
#include "ae108/cpppetsc/Vector.h"
#include "ae108/cpppetsc/scalarProduct.h"
#include "ae108/cppptest/Matchers.h"
#include <gmock/gmock.h>
using ae108::cppptest::ScalarEq;
using testing::Test;
using testing::Types;
namespace ae108 {
namespace cpppetsc {
template <class Policy> struct scalarProduct_Test : Test {
using vector_type = Vector<Policy>;
};
using Policies = Types<SequentialComputePolicy, ParallelComputePolicy>;
TYPED_TEST_CASE(scalarProduct_Test, Policies);
TYPED_TEST(scalarProduct_Test, computes_scalar_product_of_vectors) {
using vector_type = typename TestFixture::vector_type;
const auto x = tag<DistributedTag>(vector_type::fromList({3., 7.}));
const auto y = tag<DistributedTag>(vector_type::fromList({2., -4.}));
EXPECT_THAT(scalarProduct(x, y), ScalarEq(6. - 28.));
}
} // namespace cpppetsc
} // namespace ae108
\ No newline at end of file
......@@ -16,8 +16,12 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
add_library(${PROJECT_NAME}-cppslepc
Context.cc
EigenPair.cc
EigenvalueSolverDivergedException.cc
InvalidEigenvalueIndexException.cc
InvalidProblemTypeException.cc
LinearEigenvalueProblemSolver.cc
NoOperatorsSetException.cc
computeEigenvalues.cc
)
......
// © 2021 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cppslepc/EigenPair.h"
#include "ae108/cpppetsc/ParallelComputePolicy.h"
#include "ae108/cpppetsc/SequentialComputePolicy.h"
namespace ae108 {
namespace cppslepc {
template class EigenPair<cpppetsc::SequentialComputePolicy>;
template class EigenPair<cpppetsc::ParallelComputePolicy>;
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cppslepc/InvalidEigenvalueIndexException.h"
namespace ae108 {
namespace cppslepc {
const char *InvalidEigenvalueIndexException::what() const noexcept {
return "The eigenvalue with this index does not exist.";
}
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cppslepc/InvalidProblemTypeException.h"
namespace ae108 {
namespace cppslepc {
const char *InvalidProblemTypeException::what() const noexcept {
return "The problem type does not match the number of operators provided.";
}
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ae108/cppslepc/NoOperatorsSetException.h"
namespace ae108 {
namespace cppslepc {
const char *NoOperatorsSetException::what() const noexcept {
return "The matrices defining the eigenvalue problem have not been specified "
"yet.";
}
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
......@@ -29,17 +29,13 @@ computeEigenvalues(const cpppetsc::Matrix<cpppetsc::ParallelComputePolicy> &);
template std::vector<std::complex<typename LinearEigenvalueProblemSolver<
cpppetsc::SequentialComputePolicy>::real_type>>
computeGeneralizedEigenvalues(
const cpppetsc::Matrix<cpppetsc::SequentialComputePolicy> &,
const cpppetsc::Matrix<cpppetsc::SequentialComputePolicy> &,
const std::size_t number_of_eigenvalues);
computeEigenvalues(const cpppetsc::Matrix<cpppetsc::SequentialComputePolicy> &,
const cpppetsc::Matrix<cpppetsc::SequentialComputePolicy> &);
template std::vector<std::complex<typename LinearEigenvalueProblemSolver<
cpppetsc::ParallelComputePolicy>::real_type>>
computeGeneralizedEigenvalues(
const cpppetsc::Matrix<cpppetsc::ParallelComputePolicy> &A,
const cpppetsc::Matrix<cpppetsc::ParallelComputePolicy> &B,
const std::size_t number_of_eigenvalues);
computeEigenvalues(const cpppetsc::Matrix<cpppetsc::ParallelComputePolicy> &A,
const cpppetsc::Matrix<cpppetsc::ParallelComputePolicy> &B);
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2021 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#pragma once
#include "ae108/cpppetsc/Vector.h"
namespace ae108 {
namespace cppslepc {
template <class Policy> struct EigenPair {
using vector_type = cpppetsc::distributed<cpppetsc::Vector<Policy>>;
using value_type = typename vector_type::value_type::value_type;
#ifdef AE108_PETSC_COMPLEX
value_type value;
vector_type vector;
#else
std::complex<value_type> value;
vector_type vector_real;
vector_type vector_imag;
#endif
};
extern template class EigenPair<cpppetsc::SequentialComputePolicy>;
extern template class EigenPair<cpppetsc::ParallelComputePolicy>;
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#pragma once
#include "ae108/cpppetsc/Exception.h"
namespace ae108 {
namespace cppslepc {
struct InvalidEigenvalueIndexException : cpppetsc::Exception {
const char *what() const noexcept override;
};
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
// © 2022 ETH Zurich, Mechanics and Materials Lab
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#pragma once
#include "ae108/cpppetsc/Exception.h"
namespace ae108 {
namespace cppslepc {
struct InvalidProblemTypeException : cpppetsc::Exception {
const char *what() const noexcept override;
};
} // namespace cppslepc
} // namespace ae108
\ No newline at end of file
......@@ -18,7 +18,10 @@
#include "ae108/cpppetsc/ParallelComputePolicy_fwd.h"
#include "ae108/cpppetsc/SequentialComputePolicy_fwd.h"
#include "ae108/cpppetsc/UniqueEntity.h"
#include <slepceps.h>
#include "ae108/cpppetsc/Vector.h"
#include "ae108/cppslepc/EigenPair.h"
#include "ae108/cppslepc/EigenvalueProblemSolverDivergedException.h"
#include <slepc/slepceps.h>
namespace ae108 {
namespace cppslepc {
......@@ -28,15 +31,80 @@ public:
using size_type = PetscInt;
using value_type = PetscScalar;
using real_type = PetscReal;
using complex_type = PetscComplex;
using vector_type = cpppetsc::distributed<cpppetsc::Vector<Policy>>;
using matrix_type = cpppetsc::Matrix<Policy>;
explicit LinearEigenvalueProblemSolver();
enum class Type {
hermitian = EPS_HEP,
nonhermitian = EPS_NHEP,
generalized_hermitian = EPS_GHEP,
generalized_nonhermitian = EPS_GNHEP,
generalized_nonhermitian_spd = EPS_PGNHEP,
generalized_indefinite = EPS_GHIEP
};
/**
* @brief Sets the matrices associated with a standard eigenvalue problem,
* i.e. Ax = lambda * x.
* @throw InvalidProblemTypeException if the type does not represent a
* nongeneralized problem
*/
void setOperators(const matrix_type *A, const Type type = Type::nonhermitian);
/**
* @brief Sets the matrices associated with a generalized eigenvalue
* problem, i.e. Ax = lambda * Bx.
* @throw InvalidProblemTypeException if the type does not represent a
* generalized problem
*/
void setOperators(const matrix_type *A, const matrix_type *B,
const Type type = Type::generalized_nonhermitian);
/**
* @brief Solve linear eigenvalue problem.
* @throw NoOperatorsSetException if the operators have not been set before
* the call to solve()
*/
void solve();
/**
* @brief Number of discovered eigenpairs.
*/
size_type numberOfEigenpairs() const;
/**
* @brief Get nth eigenpair.
* @throw InvalidEigenvalueIndexException if `n` is out of range
*/
void getEigenpair(const size_type n, EigenPair<Policy> *out) const;
/**
* @brief Get nth eigenvalue.
* @throw InvalidEigenvalueIndexException if `n` is out of range
*/
complex_type getEigenvalue(const size_type n) const;
/**
* @brief Returns the internal solver.
*/
EPS data() const;
private:
enum class GeneralizedProblem {
no = 0,
yes = 1,
};
/**
* @brief Sets the type of the eigenvalue problem.
* @throw InvalidProblemTypeException if the type does not match the state of
* `generalized`
*/
void setType(const Type type, const GeneralizedProblem generalized);
cpppetsc::UniqueEntity<EPS> eps_;
};
......@@ -48,6 +116,10 @@ extern template class LinearEigenvalueProblemSolver<
} // namespace cppslepc
} // namespace ae108
#include "ae108/cppslepc/InvalidEigenvalueIndexException.h"
#include "ae108/cppslepc/InvalidProblemTypeException.h"
#include "ae108/cppslepc/NoOperatorsSetException.h"
namespace ae108 {
namespace cppslepc {
......@@ -59,7 +131,112 @@ LinearEigenvalueProblemSolver<Policy>::LinearEigenvalueProblemSolver()
return cpppetsc::UniqueEntity<EPS>(
solver, [](EPS eps) { Policy::handleError(EPSDestroy(&eps)); });
}()) {
Policy::handleError(EPSSetFromOptions(eps_.get()));
Policy::handleError(EPSSetFromOptions(this->data()));
}
template <class Policy>
void LinearEigenvalueProblemSolver<Policy>::setOperators(const matrix_type *A,
const Type type) {
Policy::handleError(EPSSetOperators(this->data(), A->data(), NULL));
setType(type, GeneralizedProblem::no);
}
template <class Policy>
void LinearEigenvalueProblemSolver<Policy>::setOperators(const matrix_type *A,
const matrix_type *B,
const Type type) {
Policy::handleError(EPSSetOperators(this->data(), A->data(), B->data()));
setType(type, GeneralizedProblem::yes);
}
template <class Policy> void LinearEigenvalueProblemSolver<Policy>::solve() {
const auto numberOfMatrices = [&]() {
auto st = ST();
Policy::handleError(EPSGetST(data(), &st));
auto n = size_type();
Policy::handleError(STGetNumMatrices(st, &n));
return n;
};
if (numberOfMatrices() == 0) {
throw NoOperatorsSetException();
}
Policy::handleError(EPSSolve(this->data()));
const auto hasError = [&]() {
auto reason = EPSConvergedReason{};
Policy::handleError(EPSGetConvergedReason(this->data(), &reason));
return reason < 0;
}();
if (hasError) {
throw EigenvalueProblemSolverDivergedException{};
}
}
template <class Policy>
typename LinearEigenvalueProblemSolver<Policy>::size_type
LinearEigenvalueProblemSolver<Policy>::numberOfEigenpairs() const {
auto size = size_type{};
Policy::handleError(EPSGetConverged(this->data(), &size));
return size;
}
template <class Policy>
void LinearEigenvalueProblemSolver<Policy>::setType(
const Type type, const GeneralizedProblem generalized) {
const auto generalizedType =
bool{type == Type::generalized_hermitian ||
type == Type::generalized_nonhermitian ||
type == Type::generalized_nonhermitian_spd ||
type == Type::generalized_indefinite};
if (generalizedType != (generalized == GeneralizedProblem::yes)) {
throw InvalidProblemTypeException();
}
Policy::handleError(
EPSSetProblemType(this->data(), static_cast<EPSProblemType>(type)));
}
template <class Policy>
void LinearEigenvalueProblemSolver<Policy>::getEigenpair(
const size_type n, EigenPair<Policy> *out) const {
if (n < 0 || n >= numberOfEigenpairs()) {
throw InvalidEigenvalueIndexException();
}
#ifdef AE108_PETSC_COMPLEX
Policy::handleError(EPSGetEigenpair(this->data(), n, &out->value, NULL,
out->vector.unwrap().data(), NULL));
#else
value_type real, imag;
Policy::handleError(EPSGetEigenpair(this->data(), n, &real, &imag,
out->vector_real.unwrap().data(),
out->vector_imag.unwrap().data()));
out->value.real(real);
out->value.imag(imag);
#endif
}
template <class Policy>
typename LinearEigenvalueProblemSolver<Policy>::complex_type
LinearEigenvalueProblemSolver<Policy>::getEigenvalue(const size_type n) const {
if (n < 0 || n >= numberOfEigenpairs()) {
throw InvalidEigenvalueIndexException();
}
auto eigenvalue = complex_type();
#ifdef AE108_PETSC_COMPLEX