From b56a1d1d57c977e645e6cd65ce7bd21c02bc0920 Mon Sep 17 00:00:00 2001
From: Manuel Weberndorfer <manuel.weberndorfer@id.ethz.ch>
Date: Tue, 26 Oct 2021 11:01:14 +0000
Subject: [PATCH] add cpppetsc::LeastSquaresSolver
---
cpppetsc/src/CMakeLists.txt | 1 +
cpppetsc/src/LeastSquaresSolver.cc | 32 +++
.../ae108/cpppetsc/LeastSquaresSolver.h | 186 ++++++++++++++++++
cpppetsc/test/CMakeLists.txt | 1 +
cpppetsc/test/LeastSquaresSolver_Test.cc | 88 +++++++++
5 files changed, 308 insertions(+)
create mode 100644 cpppetsc/src/LeastSquaresSolver.cc
create mode 100644 cpppetsc/src/include/ae108/cpppetsc/LeastSquaresSolver.h
create mode 100644 cpppetsc/test/LeastSquaresSolver_Test.cc
diff --git a/cpppetsc/src/CMakeLists.txt b/cpppetsc/src/CMakeLists.txt
index 8249b2b2..b0761fa4 100644
--- a/cpppetsc/src/CMakeLists.txt
+++ b/cpppetsc/src/CMakeLists.txt
@@ -20,6 +20,7 @@ add_library(${PROJECT_NAME}-cpppetsc
./Exception.cc
./GeneralizedMeshBoundaryConditition.cc
./InvalidParametersException.cc
+ ./LeastSquaresSolver.cc
./LinearSolver.cc
./LinearSolverDivergedException.cc
./LocalVertexView.cc
diff --git a/cpppetsc/src/LeastSquaresSolver.cc b/cpppetsc/src/LeastSquaresSolver.cc
new file mode 100644
index 00000000..cc9f275a
--- /dev/null
+++ b/cpppetsc/src/LeastSquaresSolver.cc
@@ -0,0 +1,32 @@
+// © 2020, 2021 ETH Zurich, Mechanics and Materials Lab
+// © 2020 California Institute of Technology
+//
+// 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/LeastSquaresSolver.h"
+
+#ifndef AE108_PETSC_COMPLEX
+
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/SequentialComputePolicy.h"
+
+namespace ae108 {
+namespace cpppetsc {
+
+template class LeastSquaresSolver<SequentialComputePolicy>;
+template class LeastSquaresSolver<ParallelComputePolicy>;
+
+} // namespace cpppetsc
+} // namespace ae108
+
+#endif
\ No newline at end of file
diff --git a/cpppetsc/src/include/ae108/cpppetsc/LeastSquaresSolver.h b/cpppetsc/src/include/ae108/cpppetsc/LeastSquaresSolver.h
new file mode 100644
index 00000000..1be860cb
--- /dev/null
+++ b/cpppetsc/src/include/ae108/cpppetsc/LeastSquaresSolver.h
@@ -0,0 +1,186 @@
+// © 2020, 2021 ETH Zurich, Mechanics and Materials Lab
+// © 2020 California Institute of Technology
+//
+// 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
+
+#ifndef AE108_PETSC_COMPLEX
+
+#include "ae108/cpppetsc/Matrix.h"
+#include "ae108/cpppetsc/ParallelComputePolicy_fwd.h"
+#include "ae108/cpppetsc/SequentialComputePolicy_fwd.h"
+#include "ae108/cpppetsc/TAOSolverDivergedException.h"
+#include "ae108/cpppetsc/TaggedVector.h"
+#include "ae108/cpppetsc/UniqueEntity.h"
+#include "ae108/cpppetsc/Vector.h"
+#include <functional>
+#include <petscmat.h>
+#include <petscmath.h>
+#include <petscsys.h>
+#include <petsctao.h>
+#include <petscvec.h>
+
+namespace ae108 {
+namespace cpppetsc {
+
+template <class Policy> class LeastSquaresSolver {
+public:
+ using vector_type = Vector<Policy>;
+ using matrix_type = Matrix<Policy>;
+
+ using size_type = typename Vector<Policy>::size_type;
+ using value_type = typename Vector<Policy>::value_type;
+ using real_type = typename Vector<Policy>::real_type;
+
+ /**
+ * @brief Initializes the nonlinear solver with the problem dimension.
+ *
+ * @param bufferMatrix The matrix to use as an internal buffer.
+ * @param bufferVector The matrix to use as an internal buffer.
+ */
+ explicit LeastSquaresSolver(matrix_type bufferMatrix,
+ distributed<vector_type> bufferVector);
+
+ /**
+ * @remark The first parameter is the input, the second parameter is the
+ * output.
+ */
+ using ObjectiveFunctor = std::function<void(const distributed<vector_type> &,
+ distributed<vector_type> *)>;
+
+ /**
+ * @remark The first parameter is the input, the second parameter is the
+ * output.
+ */
+ using GradientFunctor =
+ std::function<void(const distributed<vector_type> &, matrix_type *)>;
+
+ /**
+ * @brief Minimizes the sum of the squared elements of the objective function.
+ *
+ * @param objective The function to minimize.
+ * @param gradient The gradient of the objective.
+ * @param guess The starting point for the optimization.
+ *
+ * @throw TAOSolverDivergedException if solve did not converge.
+ *
+ * @return The minimizer.
+ */
+ distributed<vector_type> solve(ObjectiveFunctor objective,
+ GradientFunctor gradient,
+ distributed<vector_type> guess) const;
+
+private:
+ static PetscErrorCode objectiveAdapter(Tao, Vec, Vec, void *);
+ static PetscErrorCode gradientAdapter(Tao, Vec, Mat, Mat, void *);
+
+ static Tao createTao();
+
+ /**
+ * @throw TAOSolverDivergedException if solve did not converge.
+ */
+ void checkConvergence() const;
+
+ UniqueEntity<Tao> _tao;
+
+ matrix_type _buffer_matrix;
+ distributed<vector_type> _buffer_vector;
+};
+
+extern template class LeastSquaresSolver<SequentialComputePolicy>;
+extern template class LeastSquaresSolver<ParallelComputePolicy>;
+} // namespace cpppetsc
+} // namespace ae108
+
+/********************************************************************
+ * implementations
+ *******************************************************************/
+
+namespace ae108 {
+namespace cpppetsc {
+template <class Policy> Tao LeastSquaresSolver<Policy>::createTao() {
+ Tao tao;
+ Policy::handleError(TaoCreate(Policy::communicator(), &tao));
+ return tao;
+}
+
+template <class Policy>
+LeastSquaresSolver<Policy>::LeastSquaresSolver(
+ matrix_type bufferMatrix, distributed<vector_type> bufferVector)
+ : _tao(makeUniqueEntity<Policy>(createTao())),
+ _buffer_matrix(std::move(bufferMatrix)),
+ _buffer_vector(std::move(bufferVector)) {
+ Policy::handleError(TaoSetFromOptions(_tao.get()));
+ Policy::handleError(TaoSetType(_tao.get(), TAOBRGN));
+ Policy::handleError(TaoBRGNSetRegularizerWeight(_tao.get(), 0.));
+}
+
+template <class Policy>
+void LeastSquaresSolver<Policy>::checkConvergence() const {
+ const auto errorCode = [this]() {
+ auto reason = TaoConvergedReason{};
+ Policy::handleError(TaoGetConvergedReason(_tao.get(), &reason));
+ return reason;
+ }();
+
+ if (errorCode < 0) {
+ throw TAOSolverDivergedException{};
+ }
+}
+
+template <class Policy>
+distributed<typename LeastSquaresSolver<Policy>::vector_type>
+LeastSquaresSolver<Policy>::solve(ObjectiveFunctor objective,
+ GradientFunctor gradient,
+ distributed<vector_type> guess) const {
+ Policy::handleError(
+ TaoSetResidualRoutine(_tao.get(), _buffer_vector.unwrap().data(),
+ &LeastSquaresSolver::objectiveAdapter, &objective));
+ Policy::handleError(TaoSetJacobianResidualRoutine(
+ _tao.get(), _buffer_matrix.data(), _buffer_matrix.data(),
+ &LeastSquaresSolver::gradientAdapter, &gradient));
+ Policy::handleError(TaoSetInitialVector(_tao.get(), guess.unwrap().data()));
+ Policy::handleError(TaoSolve(_tao.get()));
+
+ checkConvergence();
+
+ return guess;
+}
+
+template <class Policy>
+PetscErrorCode LeastSquaresSolver<Policy>::objectiveAdapter(Tao, Vec input,
+ Vec output,
+ void *context) {
+ const auto functionPtr = static_cast<ObjectiveFunctor *>(context);
+ distributed<vector_type> wrappedOutput(UniqueEntity<Vec>(output, [](Vec) {}));
+ (*functionPtr)(distributed<vector_type>(UniqueEntity<Vec>(input, [](Vec) {})),
+ &wrappedOutput);
+ return 0;
+}
+
+template <class Policy>
+PetscErrorCode LeastSquaresSolver<Policy>::gradientAdapter(Tao, Vec input,
+ Mat output, Mat,
+ void *context) {
+ const auto functionPtr = static_cast<GradientFunctor *>(context);
+ matrix_type wrappedOutput(UniqueEntity<Mat>(output, [](Mat) {}));
+ (*functionPtr)(distributed<vector_type>(UniqueEntity<Vec>(input, [](Vec) {})),
+ &wrappedOutput);
+ return 0;
+}
+
+} // namespace cpppetsc
+} // namespace ae108
+
+#endif
\ No newline at end of file
diff --git a/cpppetsc/test/CMakeLists.txt b/cpppetsc/test/CMakeLists.txt
index b7cf6dc4..2167eeab 100644
--- a/cpppetsc/test/CMakeLists.txt
+++ b/cpppetsc/test/CMakeLists.txt
@@ -17,6 +17,7 @@ cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
add_executable(${PROJECT_NAME}-CppPetscTests
./InvalidParametersException_Test.cc
+ ./LeastSquaresSolver_Test.cc
./LinearSolverDivergedException_Test.cc
./LinearSolver_Test.cc
./Matrix_Test.cc
diff --git a/cpppetsc/test/LeastSquaresSolver_Test.cc b/cpppetsc/test/LeastSquaresSolver_Test.cc
new file mode 100644
index 00000000..a5a7d4e1
--- /dev/null
+++ b/cpppetsc/test/LeastSquaresSolver_Test.cc
@@ -0,0 +1,88 @@
+// © 2020, 2021 ETH Zurich, Mechanics and Materials Lab
+// © 2020 California Institute of Technology
+//
+// 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/LeastSquaresSolver.h"
+#include "ae108/cpppetsc/ParallelComputePolicy.h"
+#include "ae108/cpppetsc/SequentialComputePolicy.h"
+#include "ae108/cpppetsc/TaggedVector.h"
+#include "ae108/cppptest/Matchers.h"
+#include "ae108/cppptest/isLocal.h"
+#include <gmock/gmock.h>
+
+using ae108::cppptest::isLocal;
+using ae108::cppptest::ScalarNear;
+using testing::Test;
+using testing::Types;
+
+namespace ae108 {
+namespace cpppetsc {
+namespace {
+
+template <class Policy> struct LeastSquaresSolver_Test : Test {
+ using solver_type = LeastSquaresSolver<Policy>;
+ using vector_type = typename solver_type::vector_type;
+ using matrix_type = typename solver_type::matrix_type;
+
+ /**
+ * @brief Solves {x == 1, 2x == 4} by minimizing the squared distance.
+ *
+ * @return The solution.
+ */
+ distributed<vector_type>
+ solveConflictingEquations(const solver_type &solver) const {
+ return solver.solve(
+ [](const distributed<vector_type> &input,
+ distributed<vector_type> *const output) {
+ const auto full = vector_type::fromDistributed(input);
+ const auto replacer = output->unwrap().replace();
+ if (isLocal(output->unwrap(), 0)) {
+ replacer(0) = full(0) - 1.;
+ }
+ if (isLocal(output->unwrap(), 1)) {
+ replacer(1) = 2. * full(0) - 4.;
+ }
+ },
+ [](const distributed<vector_type> &, matrix_type *const output) {
+ const auto replacer = output->assemblyView().replace();
+ if (isLocal(*output, 0)) {
+ replacer(0, 0) = 1.;
+ }
+ if (isLocal(*output, 1)) {
+ replacer(1, 0) = 2.;
+ }
+ },
+ tag<DistributedTag>(vector_type::fromList({7.})));
+ }
+};
+
+using Policies = Types<SequentialComputePolicy, ParallelComputePolicy>;
+TYPED_TEST_CASE(LeastSquaresSolver_Test, Policies);
+
+TYPED_TEST(LeastSquaresSolver_Test,
+ finds_optimal_solution_for_conflicting_equations) {
+ using matrix_type = typename TestFixture::solver_type::matrix_type;
+ using vector_type = typename TestFixture::solver_type::vector_type;
+
+ typename TestFixture::solver_type solver(matrix_type(2, 1),
+ distributed<vector_type>(2));
+ const auto solution =
+ vector_type::fromDistributed(this->solveConflictingEquations(solver));
+
+ EXPECT_THAT(solution(0), ScalarNear(18. / 10., 1e-6));
+}
+
+} // namespace
+} // namespace cpppetsc
+} // namespace ae108
--
GitLab