TAOSolver_Test.cc 10.3 KB
Newer Older
1
// © 2020, 2022 ETH Zurich, Mechanics and Materials Lab
webmanue's avatar
webmanue committed
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// © 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.

webmanue's avatar
webmanue committed
16
17
#ifndef AE108_PETSC_COMPLEX

webmanue's avatar
webmanue committed
18
19
20
21
22
23
24
25
26
27
#include "ae108/cpppetsc/ParallelComputePolicy.h"
#include "ae108/cpppetsc/SequentialComputePolicy.h"
#include "ae108/cpppetsc/TAOSolver.h"
#include "ae108/cpppetsc/TaggedVector.h"
#include "ae108/cppptest/Matchers.h"
#include "ae108/cppptest/isLocal.h"
#include <cmath>
#include <gmock/gmock.h>

using ae108::cppptest::isLocal;
28
using ae108::cppptest::ScalarEqIfLocal;
webmanue's avatar
webmanue committed
29
using ae108::cppptest::ScalarNearIfLocal;
webmanue's avatar
webmanue committed
30
31
32
33
34
35
36
37
38
39
using testing::Test;
using testing::Types;

namespace ae108 {
namespace cpppetsc {
namespace {

template <class Policy> struct TAOSolver_Test : Test {
  using solver_type = TAOSolver<Policy>;
  using value_type = typename solver_type::value_type;
webmanue's avatar
webmanue committed
40
  using real_type = typename solver_type::real_type;
webmanue's avatar
webmanue committed
41
42
43
44
45
46
47
48
  using vector_type = typename solver_type::vector_type;
  using matrix_type = typename solver_type::matrix_type;
};

using Policies = Types<SequentialComputePolicy, ParallelComputePolicy>;
TYPED_TEST_CASE(TAOSolver_Test, Policies);

TYPED_TEST(TAOSolver_Test, minimizes_x_minus_1_squared) {
webmanue's avatar
webmanue committed
49
  using real_type = typename TestFixture::real_type;
webmanue's avatar
webmanue committed
50
51
52
53
54
  using vector_type = typename TestFixture::vector_type;
  using matrix_type = typename TestFixture::matrix_type;

  typename TestFixture::solver_type solver(matrix_type(1, 1));
  const auto solution = solver.solve(
webmanue's avatar
webmanue committed
55
      [](const distributed<vector_type> &input, real_type *const output) {
webmanue's avatar
webmanue committed
56
        const auto full = vector_type::fromDistributed(input);
webmanue's avatar
webmanue committed
57
        *output = std::norm(full(0) - 1.);
webmanue's avatar
webmanue committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
      },
      [](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) = 2. * (full(0) - 1.);
        }
      },
      [](const distributed<vector_type> &, matrix_type *const output) {
        const auto replacer = output->assemblyView().replace();
        if (isLocal(*output, 0)) {
          replacer(0, 0) = 2.;
        }
      },
      tag<DistributedTag>(vector_type::fromList({7.})));

75
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, 1.));
webmanue's avatar
webmanue committed
76
77
78
}

TYPED_TEST(TAOSolver_Test, minimizes_two_dimensional_problem) {
webmanue's avatar
webmanue committed
79
  using real_type = typename TestFixture::real_type;
webmanue's avatar
webmanue committed
80
81
82
83
84
  using vector_type = typename TestFixture::vector_type;
  using matrix_type = typename TestFixture::matrix_type;

  typename TestFixture::solver_type solver(matrix_type(2, 2));
  const auto solution = solver.solve(
webmanue's avatar
webmanue committed
85
      [](const distributed<vector_type> &input, real_type *const output) {
webmanue's avatar
webmanue committed
86
        const auto full = vector_type::fromDistributed(input);
webmanue's avatar
webmanue committed
87
88
        *output =
            std::norm(full(0) - 1) + std::norm(2. * full(0) - full(1)) + 7;
webmanue's avatar
webmanue committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
      },
      [](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) = 2. * (full(0) - 1.) + 4. * (2. * full(0) - full(1));
        }
        if (isLocal(output->unwrap(), 1)) {
          replacer(1) = -2. * (2. * full(0) - full(1));
        }
      },
      [](const distributed<vector_type> &, matrix_type *const output) {
        const auto replacer = output->assemblyView().replace();
        if (isLocal(*output, 0)) {
          replacer(0, 0) = 10.;
          replacer(0, 1) = -4.;
        }
        if (isLocal(*output, 1)) {
          replacer(1, 0) = -4.;
          replacer(1, 1) = 2.;
        }
      },
      tag<DistributedTag>(vector_type::fromList({3., 4.})));

114
115
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, 1.));
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(1, 2.));
webmanue's avatar
webmanue committed
116
117
118
}

TYPED_TEST(TAOSolver_Test, raises_exception_on_nonconvergence) {
webmanue's avatar
webmanue committed
119
  using real_type = typename TestFixture::real_type;
webmanue's avatar
webmanue committed
120
121
122
123
124
125
  using vector_type = typename TestFixture::vector_type;
  using matrix_type = typename TestFixture::matrix_type;

  typename TestFixture::solver_type solver(matrix_type(1, 1));
  const auto callSolve = [&solver]() {
    solver.solve(
webmanue's avatar
webmanue committed
126
        [](const distributed<vector_type> &, real_type *const output) {
webmanue's avatar
webmanue committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
          *output = 1.;
        },
        [](const distributed<vector_type> &,
           distributed<vector_type> *const output) {
          const auto replacer = output->unwrap().replace();
          replacer(0) = 1.;
        },
        [](const distributed<vector_type> &, matrix_type *const output) {
          const auto replacer = output->assemblyView().replace();
          replacer(0, 0) = 0.;
        },
        tag<DistributedTag>(vector_type::fromList({1.})));
  };

  EXPECT_THROW(callSolve(), TAOSolverDivergedException);
}

144
145
146
147
148
149
TYPED_TEST(TAOSolver_Test, minimizes_with_equality_constraints) {
  using real_type = typename TestFixture::real_type;
  using vector_type = typename TestFixture::vector_type;
  using matrix_type = typename TestFixture::matrix_type;

  typename TestFixture::solver_type solver(
webmanue's avatar
webmanue committed
150
      matrix_type(2, 2), TestFixture::solver_type::Type::almm);
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187

  const auto value = 5.;

  solver.setConstraints({
      [&](const distributed<vector_type> &input,
          distributed<vector_type> *const output) {
        const auto full = vector_type::fromDistributed(input);
        output->unwrap().setZero();
        output->unwrap().replace().element(0, full(0) - value);
      },
      [](const distributed<vector_type> &, matrix_type *const output) {
        output->setZero();
        const auto replacer =
            output->assemblyView().replace().element(0, 0, 1.);
      },
      {tag<DistributedTag>(vector_type(1)), matrix_type(1, 2)},
  });

  const auto solution = solver.solve(
      [](const distributed<vector_type> &input, real_type *const output) {
        const auto full = vector_type::fromDistributed(input);
        *output = std::pow(full(0) - 1., 2.) + std::pow(full(1) - 2., 2.);
      },
      [](const distributed<vector_type> &input,
         distributed<vector_type> *const output) {
        const auto full = vector_type::fromDistributed(input);
        const auto replacer = output->unwrap().replace();
        replacer(0) = 2. * (full(0) - 1.);
        replacer(1) = 2. * (full(1) - 2.);
      },
      [](const distributed<vector_type> &, matrix_type *const output) {
        const auto replacer = output->assemblyView().replace();
        replacer(0, 0) = 2.;
        replacer(1, 1) = 2.;
      },
      tag<DistributedTag>(vector_type::fromList({3., 4.})));

webmanue's avatar
webmanue committed
188
  EXPECT_THAT(solution.unwrap(), ScalarNearIfLocal(0, value, 1e-7));
189
190
191
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(1, 2.));
}

webmanue's avatar
webmanue committed
192
193
194
template <class Policy> struct TAOSolver_BoundsTest : Test {
  using solver_type = TAOSolver<Policy>;
  using value_type = typename solver_type::value_type;
webmanue's avatar
webmanue committed
195
  using real_type = typename solver_type::real_type;
webmanue's avatar
webmanue committed
196
197
198
199
200
201
202
  using vector_type = typename solver_type::vector_type;
  using matrix_type = typename solver_type::matrix_type;

  solver_type solver{matrix_type(1, 1), solver_type::Type::blmvm};

  distributed<vector_type> minimizeQuadratic(const value_type guess) {
    return solver.solve(
webmanue's avatar
webmanue committed
203
        [](const distributed<vector_type> &input, real_type *const output) {
webmanue's avatar
webmanue committed
204
          const auto full = vector_type::fromDistributed(input);
webmanue's avatar
webmanue committed
205
          *output = std::norm(full(0) - 1.);
webmanue's avatar
webmanue committed
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
        },
        [](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) = 2. * (full(0) - 1.);
          }
        },
        [](const distributed<vector_type> &, matrix_type *const output) {
          const auto replacer = output->assemblyView().replace();
          if (isLocal(*output, 0)) {
            replacer(0, 0) = 2.;
          }
        },
        tag<DistributedTag>(vector_type::fromList({guess})));
  }
};

using Policies = Types<SequentialComputePolicy, ParallelComputePolicy>;
TYPED_TEST_CASE(TAOSolver_BoundsTest, Policies);

TYPED_TEST(TAOSolver_BoundsTest, minimizes_with_lower_bounds) {
  using vector_type = typename TestFixture::vector_type;
  this->solver.setBounds(tag<DistributedTag>(vector_type::fromList({2.})),
                         tag<DistributedTag>(vector_type::fromList(
                             {TestFixture::solver_type::no_upper_bound})));

  const auto solution = this->minimizeQuadratic(7.);
235
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, 2.));
webmanue's avatar
webmanue committed
236
237
238
239
240
241
242
243
244
}

TYPED_TEST(TAOSolver_BoundsTest, minimizes_with_upper_bounds) {
  using vector_type = typename TestFixture::vector_type;
  this->solver.setBounds(tag<DistributedTag>(vector_type::fromList(
                             {TestFixture::solver_type::no_lower_bound})),
                         tag<DistributedTag>(vector_type::fromList({-2.})));

  const auto solution = this->minimizeQuadratic(-7.);
245
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, -2.));
webmanue's avatar
webmanue committed
246
247
248
249
250
251
252
253
}

TYPED_TEST(TAOSolver_BoundsTest, minimizes_with_both_bounds) {
  using vector_type = typename TestFixture::vector_type;
  this->solver.setBounds(tag<DistributedTag>(vector_type::fromList({2.})),
                         tag<DistributedTag>(vector_type::fromList({2.})));

  const auto solution = this->minimizeQuadratic(-7.);
254
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, 2.));
webmanue's avatar
webmanue committed
255
256
257
258
259
260
261
262
263
}

TYPED_TEST(TAOSolver_BoundsTest, lmvm_does_not_apply_bounds) {
  this->solver.setType(TestFixture::solver_type::Type::lmvm);
  using vector_type = typename TestFixture::vector_type;
  this->solver.setBounds(tag<DistributedTag>(vector_type::fromList({2.})),
                         tag<DistributedTag>(vector_type::fromList({2.})));

  const auto solution = this->minimizeQuadratic(-7.);
264
  EXPECT_THAT(solution.unwrap(), ScalarEqIfLocal(0, 1.));
webmanue's avatar
webmanue committed
265
266
267
268
}
} // namespace
} // namespace cpppetsc
} // namespace ae108
webmanue's avatar
webmanue committed
269
270

#endif