ROL
poisson-control/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_Bounds.hpp"
53#include "ROL_Algorithm.hpp"
56#include "ROL_StatusTest.hpp"
57#include "ROL_Types.hpp"
58
59#include "ROL_Stream.hpp"
60#include "Teuchos_GlobalMPISession.hpp"
61#include "Teuchos_XMLParameterListHelpers.hpp"
62
63#include <iostream>
64#include <algorithm>
65
66
67
68template <class Real>
69class StatusTest_PDAS : public ROL::StatusTest<Real> {
70private:
71
72 Real gtol_;
73 Real stol_;
75
76public:
77
78 virtual ~StatusTest_PDAS() {}
79
80 StatusTest_PDAS( Real gtol = 1.e-6, Real stol = 1.e-12, int max_iter = 100 ) :
81 gtol_(gtol), stol_(stol), max_iter_(max_iter) {}
82
85 virtual bool check( ROL::AlgorithmState<Real> &state ) {
86 if ( (state.gnorm > this->gtol_) &&
87 (state.snorm > this->stol_) &&
88 (state.iter < this->max_iter_) ) {
89 return true;
90 }
91 else {
92 if ( state.iter < 2 ) {
93 return true;
94 }
95 return false;
96 }
97 }
98
99};
100
101typedef double RealT;
102
103int main(int argc, char *argv[]) {
104
105 typedef std::vector<RealT> vector;
106 typedef ROL::Vector<RealT> V;
107 typedef ROL::StdVector<RealT> SV;
108
109 typedef typename vector::size_type uint;
110
111
112
113 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
114
115 // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
116 int iprint = argc - 1;
117 ROL::Ptr<std::ostream> outStream;
118 ROL::nullstream bhs; // outputs nothing
119 if (iprint > 0)
120 outStream = ROL::makePtrFromRef(std::cout);
121 else
122 outStream = ROL::makePtrFromRef(bhs);
123
124 int errorFlag = 0;
125
126 // *** Example body.
127
128 try {
129 uint dim = 256; // Set problem dimension.
130 RealT alpha = 1.e-4;
132
133 ROL::Ptr<vector> l_ptr = ROL::makePtr<vector>(dim);
134 ROL::Ptr<vector> u_ptr = ROL::makePtr<vector>(dim);
135
136 ROL::Ptr<V> lo = ROL::makePtr<SV>(l_ptr);
137 ROL::Ptr<V> up = ROL::makePtr<SV>(u_ptr);
138
139 for ( uint i = 0; i < dim; i++ ) {
140 if ( i < dim/3.0 || i > 2*dim/3.0 ) {
141 (*l_ptr)[i] = 0.0;
142 (*u_ptr)[i] = 0.25;
143 }
144 else {
145 (*l_ptr)[i] = 0.75;
146 (*u_ptr)[i] = 1.0;
147 }
148 }
149 ROL::Bounds<RealT> icon(lo,up);
150
151 // Primal dual active set.
152 std::string filename = "input.xml";
153 auto parlist = ROL::getParametersFromXmlFile( filename );
154
155 // Krylov parameters.
156 parlist->sublist("General").sublist("Krylov").set("Absolute Tolerance",1.e-4);
157 parlist->sublist("General").sublist("Krylov").set("Relative Tolerance",1.e-2);
158 parlist->sublist("General").sublist("Krylov").set("Iteration Limit",50);
159
160 // PDAS parameters.
161 parlist->sublist("Step").sublist("Primal Dual Active Set").set("Relative Step Tolerance",1.e-8);
162 parlist->sublist("Step").sublist("Primal Dual Active Set").set("Relative Gradient Tolerance",1.e-6);
163 parlist->sublist("Step").sublist("Primal Dual Active Set").set("Iteration Limit", 1);
164 parlist->sublist("Step").sublist("Primal Dual Active Set").set("Dual Scaling",(alpha>0.0)?alpha:1.e-4);
165
166 // Status test parameters.
167 parlist->sublist("Status Test").set("Gradient Tolerance",1.e-12);
168 parlist->sublist("Status Test").set("Step Tolerance",1.e-14);
169 parlist->sublist("Status Test").set("Iteration Limit",100);
170
171 // Define algorithm.
172 ROL::Ptr<ROL::Step<RealT>> step = ROL::makePtr<ROL::PrimalDualActiveSetStep<RealT>>(*parlist);
173 ROL::Ptr<ROL::StatusTest<RealT>> status = ROL::makePtr<ROL::StatusTest<RealT>>(*parlist);
174 ROL::Ptr<ROL::Algorithm<RealT>> algo = ROL::makePtr<ROL::Algorithm<RealT>>(step,status,false);
175
176 // Iteration vector.
177 ROL::Ptr<vector> x_ptr = ROL::makePtr<vector>(dim, 0.0);
178 SV x(x_ptr);
179
180 // Run algorithm.
181 x.zero();
182 algo->run(x, obj, icon, true, *outStream);
183 std::ofstream file;
184 file.open("control_PDAS.txt");
185
186 for ( uint i = 0; i < dim; i++ ) {
187 file << (*x_ptr)[i] << "\n";
188 }
189 file.close();
190
191 // Projected Newton.
192 // re-load parameters
193 Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
194 // Set algorithm.
195 step = ROL::makePtr<ROL::TrustRegionStep<RealT>>(*parlist);
196 status = ROL::makePtr<ROL::StatusTest<RealT>>(*parlist);
197 algo = ROL::makePtr<ROL::Algorithm<RealT>>(step,status,false);
198 // Iteration vector.
199 ROL::Ptr<vector> y_ptr = ROL::makePtr<vector>(dim, 0.0);
200 SV y(y_ptr);
201
202 // Run Algorithm
203 y.zero();
204 algo->run(y, obj, icon, true, *outStream);
205
206 std::ofstream file_tr;
207 file_tr.open("control_TR.txt");
208 for ( uint i = 0; i < dim; i++ ) {
209 file_tr << (*y_ptr)[i] << "\n";
210 }
211 file_tr.close();
212
213 ROL::Ptr<V> error = x.clone();
214 error->set(x);
215 error->axpy(-1.0,y);
216 *outStream << "\nError between PDAS solution and TR solution is " << error->norm() << "\n";
217 errorFlag = ((error->norm() > 1e2*std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 1 : 0);
218 }
219 catch (std::logic_error& err) {
220 *outStream << err.what() << "\n";
221 errorFlag = -1000;
222 }; // end try
223
224 if (errorFlag != 0)
225 std::cout << "End Result: TEST FAILED\n";
226 else
227 std::cout << "End Result: TEST PASSED\n";
228
229 return 0;
230
231}
232
Vector< Real > V
Contains definitions for Poisson optimal control.
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 the elementwise interface to apply upper and lower bound constraints.
Provides an interface to check status of optimization algorithms.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual bool check(ROL::AlgorithmState< Real > &state)
Check algorithm status.
StatusTest_PDAS(Real gtol=1.e-6, Real stol=1.e-12, int max_iter=100)
int main(int argc, char *argv[])
State for algorithm class. Will be used for restarts.
constexpr auto dim