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.
Definition: ROL_Bounds.hpp:59
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.
Definition: ROL_Vector.hpp:84
Poisson distributed control.
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.
Definition: ROL_Types.hpp:143
constexpr auto dim