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_PoissonControl.hpp"
52 #include "ROL_Algorithm.hpp"
54 #include "ROL_TrustRegionStep.hpp"
55 #include "ROL_StatusTest.hpp"
56 #include "ROL_Types.hpp"
57 #include "Teuchos_oblackholestream.hpp"
58 #include "Teuchos_GlobalMPISession.hpp"
59 #include "Teuchos_XMLParameterListHelpers.hpp"
60 
61 #include <iostream>
62 #include <algorithm>
63 
64 #include "ROL_BoundConstraint.hpp"
65 
66 
67 template <class Real>
68 class StatusTest_PDAS : public ROL::StatusTest<Real> {
69 private:
70 
71  Real gtol_;
72  Real stol_;
73  int max_iter_;
74 
75 public:
76 
77  virtual ~StatusTest_PDAS() {}
78 
79  StatusTest_PDAS( Real gtol = 1.e-6, Real stol = 1.e-12, int max_iter = 100 ) :
80  gtol_(gtol), stol_(stol), max_iter_(max_iter) {}
81 
84  virtual bool check( ROL::AlgorithmState<Real> &state ) {
85  if ( (state.gnorm > this->gtol_) &&
86  (state.snorm > this->stol_) &&
87  (state.iter < this->max_iter_) ) {
88  return true;
89  }
90  else {
91  if ( state.iter < 2 ) {
92  return true;
93  }
94  return false;
95  }
96  }
97 
98 };
99 
100 typedef double RealT;
101 
102 int main(int argc, char *argv[]) {
103 
104  typedef std::vector<RealT> vector;
105  typedef ROL::Vector<RealT> V;
106  typedef ROL::StdVector<RealT> SV;
107 
108  typedef typename vector::size_type uint;
109 
110  using Teuchos::RCP; using Teuchos::rcp;
111 
112  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
113 
114  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
115  int iprint = argc - 1;
116  Teuchos::RCP<std::ostream> outStream;
117  Teuchos::oblackholestream bhs; // outputs nothing
118  if (iprint > 0)
119  outStream = Teuchos::rcp(&std::cout, false);
120  else
121  outStream = Teuchos::rcp(&bhs, false);
122 
123  int errorFlag = 0;
124 
125  // *** Example body.
126 
127  try {
128  uint dim = 256; // Set problem dimension.
129  RealT alpha = 1.e-4;
131 
132  RCP<vector> l_rcp = rcp( new vector(dim) );
133  RCP<vector> u_rcp = rcp( new vector(dim) );
134 
135  RCP<V> lo = rcp( new SV(l_rcp) );
136  RCP<V> up = rcp( new SV(u_rcp) );
137 
138  for ( uint i = 0; i < dim; i++ ) {
139  if ( i < dim/3.0 || i > 2*dim/3.0 ) {
140  (*l_rcp)[i] = 0.0;
141  (*u_rcp)[i] = 0.25;
142  }
143  else {
144  (*l_rcp)[i] = 0.75;
145  (*u_rcp)[i] = 1.0;
146  }
147  }
148  ROL::BoundConstraint<RealT> icon(lo,up);
149 
150  // Primal dual active set.
151  std::string filename = "input.xml";
152  RCP<Teuchos::ParameterList> parlist = rcp( new Teuchos::ParameterList() );
153  Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
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  RCP<ROL::Algorithm<RealT> > algo = rcp(new ROL::Algorithm<RealT>("Primal Dual Active Set",*parlist,false));
173 
174  // Iteration vector.
175  RCP<vector> x_rcp = rcp( new vector(dim, 0.0) );
176  SV x(x_rcp);
177 
178  // Run algorithm.
179  x.zero();
180  algo->run(x, obj, icon, true, *outStream);
181  std::ofstream file;
182  file.open("control_PDAS.txt");
183 
184  for ( uint i = 0; i < dim; i++ ) {
185  file << (*x_rcp)[i] << "\n";
186  }
187  file.close();
188 
189  // Projected Newton.
190  // re-load parameters
191  Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );
192  // Set algorithm.
193  algo = Teuchos::rcp(new ROL::Algorithm<RealT>("Trust Region",*parlist,false));
194  // Iteration vector.
195  RCP<vector> y_rcp = rcp( new vector(dim, 0.0) );
196  SV y(y_rcp);
197 
198  // Run Algorithm
199  y.zero();
200  algo->run(y, obj, icon, true, *outStream);
201 
202  std::ofstream file_tr;
203  file_tr.open("control_TR.txt");
204  for ( uint i = 0; i < dim; i++ ) {
205  file_tr << (*y_rcp)[i] << "\n";
206  }
207  file_tr.close();
208 
209  RCP<V> error = x.clone();
210  error->set(x);
211  error->axpy(-1.0,y);
212  *outStream << "\nError between PDAS solution and TR solution is " << error->norm() << "\n";
213  errorFlag = ((error->norm() > 1e2*std::sqrt(ROL::ROL_EPSILON<RealT>())) ? 1 : 0);
214  }
215  catch (std::logic_error err) {
216  *outStream << err.what() << "\n";
217  errorFlag = -1000;
218  }; // end try
219 
220  if (errorFlag != 0)
221  std::cout << "End Result: TEST FAILED\n";
222  else
223  std::cout << "End Result: TEST PASSED\n";
224 
225  return 0;
226 
227 }
228 
int main(int argc, char *argv[])
StatusTest_PDAS(Real gtol=1.e-6, Real stol=1.e-12, int max_iter=100)
Contains definitions of custom data types in ROL.
Contains definitions for Poisson optimal control.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:91
Provides the std::vector implementation of the ROL::Vector interface.
Poisson distributed control.
Provides an interface to run optimization algorithms.
Provides an interface to check status of optimization algorithms.
Provides the interface to apply upper and lower bound constraints.
double RealT
virtual bool check(ROL::AlgorithmState< Real > &state)
Check algorithm status.