ROL
poisson-inversion/example_01.cpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49#define USE_HESSVEC 1
50
51#include "ROL_Types.hpp"
53#include "ROL_Algorithm.hpp"
56#include "ROL_StatusTest.hpp"
57#include "ROL_Stream.hpp"
58#include "Teuchos_GlobalMPISession.hpp"
59
60#include <iostream>
61#include <algorithm>
62
63typedef double RealT;
64
65int main(int argc, char *argv[]) {
66
67 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
68
69 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
70 int iprint = argc - 1;
71 ROL::Ptr<std::ostream> outStream;
72 ROL::nullstream bhs; // outputs nothing
73 if (iprint > 0)
74 outStream = ROL::makePtrFromRef(std::cout);
75 else
76 outStream = ROL::makePtrFromRef(bhs);
77
78 int errorFlag = 0;
79
80 // *** Example body.
81
82 try {
83
84 int dim = 128; // Set problem dimension.
86
87 // Define algorithm.
88 ROL::ParameterList parlist;
89 std::string stepname = "Trust Region";
90 parlist.sublist("Step").sublist(stepname).set("Subproblem Solver", "Truncated CG");
91 parlist.sublist("General").sublist("Krylov").set("Iteration Limit",50);
92 parlist.sublist("General").sublist("Krylov").set("Relative Tolerance",1e-2);
93 parlist.sublist("General").sublist("Krylov").set("Absolute Tolerance",1e-4);
94 parlist.sublist("Status Test").set("Gradient Tolerance",1.e-12);
95 parlist.sublist("Status Test").set("Step Tolerance",1.e-14);
96 parlist.sublist("Status Test").set("Iteration Limit",100);
97 ROL::Ptr<ROL::Step<RealT>>
98 step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(parlist);
99 ROL::Ptr<ROL::StatusTest<RealT>>
100 status = ROL::makePtr<ROL::StatusTest<RealT>>(parlist);
101 ROL::Algorithm<RealT> algo(step,status,false);
102
103 // Iteration vector.
104 ROL::Ptr<std::vector<RealT> > x_ptr = ROL::makePtr<std::vector<RealT>>(dim, 0.0);
105 // Set initial guess.
106 for (int i=0; i<dim; i++) {
107 (*x_ptr)[i] = 0.1;
108 }
109 ROL::StdVector<RealT> x(x_ptr);
110
111 // Run algorithm.
112 algo.run(x, obj, true, *outStream);
113
114 // Compute dense Hessian matrix.
115 Teuchos::SerialDenseMatrix<int, RealT> H(x.dimension(), x.dimension());
116 H = ROL::computeDenseHessian<RealT>(obj, x);
117 //H.print(*outStream);
118
119 // Compute and print eigenvalues.
120 std::vector<std::vector<RealT> > eigenvals = ROL::computeEigenvalues<RealT>(H);
121
122 *outStream << "\nEigenvalues:\n";
123 for (unsigned i=0; i<(eigenvals[0]).size(); i++) {
124 if (i==0) {
125 *outStream << std::right
126 << std::setw(28) << "Real"
127 << std::setw(28) << "Imag"
128 << "\n";
129 }
130 *outStream << std::scientific << std::setprecision(16) << std::right
131 << std::setw(28) << (eigenvals[0])[i]
132 << std::setw(28) << (eigenvals[1])[i]
133 << "\n";
134 }
135
136 // Compute and print generalized eigenvalues.
137 Teuchos::SerialDenseMatrix<int, RealT> M = computeDotMatrix(x);
138 //M.print(*outStream);
139 std::vector<std::vector<RealT> > genEigenvals = ROL::computeGenEigenvalues<RealT>(H, M);
140
141 *outStream << "\nGeneralized eigenvalues:\n";
142 for (unsigned i=0; i<(genEigenvals[0]).size(); i++) {
143 if (i==0) {
144 *outStream << std::right
145 << std::setw(28) << "Real"
146 << std::setw(28) << "Imag"
147 << "\n";
148 }
149 *outStream << std::scientific << std::setprecision(16) << std::right
150 << std::setw(28) << (genEigenvals[0])[i]
151 << std::setw(28) << (genEigenvals[1])[i]
152 << "\n";
153 }
154
155 // Sort and compare eigenvalues and generalized eigenvalues - should be close.
156 std::sort((eigenvals[0]).begin(), (eigenvals[0]).end());
157 std::sort((eigenvals[1]).begin(), (eigenvals[1]).end());
158 std::sort((genEigenvals[0]).begin(), (genEigenvals[0]).end());
159 std::sort((genEigenvals[1]).begin(), (genEigenvals[1]).end());
160
161 RealT errtol = std::sqrt(ROL::ROL_EPSILON<RealT>());
162 for (unsigned i=0; i<(eigenvals[0]).size(); i++) {
163 if ( std::abs( (genEigenvals[0])[i] - (eigenvals[0])[i] ) > errtol*((eigenvals[0])[i]+ROL::ROL_THRESHOLD<RealT>()) ) {
164 errorFlag++;
165 *outStream << std::scientific << std::setprecision(20) << "Real genEigenvals - eigenvals (" << i << ") = " << std::abs( (genEigenvals[0])[i] - (eigenvals[0])[i] ) << " > " << errtol*((eigenvals[0])[i]+1e4*ROL::ROL_THRESHOLD<RealT>()) << "\n";
166 }
167 if ( std::abs( (genEigenvals[1])[i] - (eigenvals[1])[i] ) > errtol*((eigenvals[1])[i]+ROL::ROL_THRESHOLD<RealT>()) ) {
168 errorFlag++;
169 *outStream << std::scientific << std::setprecision(20) << "Imag genEigenvals - eigenvals (" << i << ") = " << std::abs( (genEigenvals[1])[i] - (eigenvals[1])[i] ) << " > " << errtol*((eigenvals[1])[i]+ROL::ROL_THRESHOLD<RealT>()) << "\n";
170 }
171 }
172
173 // Compute inverse of Hessian.
174 Teuchos::SerialDenseMatrix<int, RealT> invH = ROL::computeInverse<RealT>(H);
175 Teuchos::SerialDenseMatrix<int, RealT> HinvH(H);
176
177 // Multiply with Hessian and verify that it gives the identity (l2 dot matrix M from above).
178 HinvH.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, H, invH, 0.0);
179 //*outStream << std::scientific << std::setprecision(6); HinvH.print(*outStream);
180 HinvH -= M;
181 if (HinvH.normOne() > errtol) {
182 errorFlag++;
183 *outStream << std::scientific << std::setprecision(20) << "1-norm of H*inv(H) - I = " << HinvH.normOne() << " > " << errtol << "\n";
184 }
185
186 // Use Newton algorithm with line search.
187 stepname = "Line Search";
188 parlist.sublist("Step").sublist(stepname).sublist("Descent Method").set("Type", "Newton's Method");
189 ROL::Ptr<ROL::Step<RealT>>
190 newton_step = ROL::makePtr<ROL::LineSearchStep<RealT>>(parlist);
191 ROL::Ptr<ROL::StatusTest<RealT>>
192 newton_status = ROL::makePtr<ROL::StatusTest<RealT>>(parlist);
193 ROL::Algorithm<RealT> newton_algo(newton_step,newton_status,false);
194
195 // Reset initial guess.
196 for (int i=0; i<dim; i++) {
197 (*x_ptr)[i] = 0.1;
198 }
199
200 // Run Newton algorithm.
201 newton_algo.run(x, obj, true, *outStream);
202
203 ROL::Ptr<const ROL::AlgorithmState<RealT> > new_state = newton_algo.getState();
204 ROL::Ptr<const ROL::AlgorithmState<RealT> > old_state = algo.getState();
205 *outStream << "old_optimal_value = " << old_state->value << std::endl;
206 *outStream << "new_optimal_value = " << new_state->value << std::endl;
207 if ( std::abs(new_state->value - old_state->value) / std::abs(old_state->value) > errtol ) {
208 errorFlag++;
209 *outStream << std::scientific << std::setprecision(20) << "\nabs(new_optimal_value - old_optimal_value) / abs(old_optimal_value) = " << std::abs(new_state->value - old_state->value) / std::abs(old_state->value) << " > " << errtol << "\n";
210 }
211
212 }
213 catch (std::logic_error& err) {
214 *outStream << err.what() << "\n";
215 errorFlag = -1000;
216 }; // end try
217
218 if (errorFlag != 0)
219 std::cout << "End Result: TEST FAILED\n";
220 else
221 std::cout << "End Result: TEST PASSED\n";
222
223 return 0;
224
225}
226
Contains definitions for Poisson material inversion.
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
Contains definitions of custom data types in ROL.
Provides an interface to run optimization algorithms.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
int dimension() const
Return dimension of the vector space.
int main(int argc, char *argv[])
constexpr auto dim