ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPInitFixedFreeStd.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include <sstream>
45 
46 #include "ConstrainedOptPack_QPInitFixedFreeStd.hpp"
47 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
48 
49 namespace ConstrainedOptPack {
50 namespace QPSchurPack {
51 
53  :n_(0)
54  ,n_R_(0)
55  ,m_(0)
56  ,G_(NULL)
57  ,A_(NULL)
58  ,Ko_(NULL)
59  ,constraints_(NULL)
60 {}
61 
63  const DVectorSlice &g
64  ,const MatrixSymOp &G
65  ,const MatrixOp *A
66  ,size_type n_R
67  ,const size_type i_x_free[]
68  ,const size_type i_x_fixed[]
69  ,const EBounds bnd_fixed[]
70  ,const DVectorSlice &b_X
71  ,const MatrixSymOpNonsing &Ko
72  ,const DVectorSlice &fo
73  ,Constraints *constraints
74  ,std::ostream *out
75  ,bool test_setup
76  ,value_type warning_tol
77  ,value_type error_tol
78  ,bool print_all_warnings
79  )
80 {
81  namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
82 
83  if(!constraints)
84  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
85  "constraints == NULL is not allowed." );
86 
87  // Validate the sizes of the input arguments
88  const size_type
89  n = constraints->n(),
90  n_X = n - n_R;
91  if( n_R > n )
92  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
93  "n_R > constraints->n() is not allowed." );
94  if(g.dim() !=n)
95  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
96  "g.dim() != constraints->n()." );
97  if(G.rows() != n || G.cols() != n)
98  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
99  "G.rows() != constraints->n() or G.cols() != constraints->n()." );
100  size_type
101  m = 0;
102  if(A) {
103  m = A->cols();
104  if( A->rows() != n )
105  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
106  "A->rows() != constraints->n()." );
107  }
108  if(b_X.dim() != n_X)
109  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
110  "b_X.dim() != constraints->n() - n_R." );
111  if(Ko.rows() != n_R+m || Ko.cols() != n_R+m)
112  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
113  "Ko.rows() != n_R+A->cols() or Ko.cols() != n_R+A->cols()." );
114  if(fo.dim() != n_R+m)
115  throw std::invalid_argument( "QPInitFixedFreeStd::initialize(...) : Error, "
116  "fo.dim() != n_R+A->cols()." );
117 
118  // Setup x_init, l_x_X_map, i_x_X_map
119 
120  const int NOT_SET_YET = -9999; // Can't be FREE, LOWER, UPPER, or EQUALITY
121  if(test_setup)
122  x_init_.assign( n, (EBounds)NOT_SET_YET );
123  else
124  x_init_.resize(n);
125  l_x_X_map_.assign(n,0);
126  i_x_X_map_.assign(n_X,0);
127 
128  // Set free portion of x_init
129  if( i_x_free == NULL ) {
130  for( size_type i = 0; i < n_R; ++i )
131  x_init_[i] = FREE;
132  }
133  else {
134  if(test_setup) {
135  for( const size_type *i_x_R = i_x_free; i_x_R != i_x_free + n_R; ++i_x_R ) {
136  if( *i_x_R < 1 || *i_x_R > n ) {
137  std::ostringstream omsg;
138  omsg
139  << "QPInitFixedFreeStd::initialize(...) : Error, "
140  << "i_x_free[" << i_x_R-i_x_free << "] = "
141  << (*i_x_R) << " is out of bounds";
142  throw std::invalid_argument( omsg.str() );
143  }
144  if( x_init_(*i_x_R) != NOT_SET_YET ) {
145  std::ostringstream omsg;
146  omsg
147  << "QPInitFixedFreeStd::initialize(...) : Error, "
148  << "Duplicate entries for i_x_free[i] = "
149  << (*i_x_R);
150  throw std::invalid_argument( omsg.str() );
151  }
152  x_init_(*i_x_R) = FREE;
153  }
154  }
155  else {
156  for( const size_type *i_x_R = i_x_free; i_x_R != i_x_free + n_R; ++i_x_R ) {
157  x_init_(*i_x_R) = FREE;
158  }
159  }
160  }
161 
162  // Setup the fixed portion of x_init and l_x_X_map and i_x_X_map
163  {
164  const size_type
165  *i_x_X = i_x_fixed;
166  const EBounds
167  *bnd = bnd_fixed;
168  i_x_X_map_t::iterator
169  i_x_X_map_itr = i_x_X_map_.begin();
170  if(test_setup) {
171  for( size_type l = 1; l <= n_X; ++l, ++i_x_X, ++bnd, ++i_x_X_map_itr ) {
172  if( *i_x_X < 1 || *i_x_X > n ) {
173  std::ostringstream omsg;
174  omsg
175  << "QPInitFixedFreeStd::initialize(...) : Error, "
176  << "i_x_fixed[" << i_x_X-i_x_fixed << "] = "
177  << (*i_x_X) << " is out of bounds";
178  throw std::invalid_argument( omsg.str() );
179  }
180  if( *bnd == FREE ) {
181  std::ostringstream omsg;
182  omsg
183  << "QPInitFixedFreeStd::initialize(...) : Error, "
184  << "bnd_fixed[" << l-1 << "] can not equal FREE";
185  throw std::invalid_argument( omsg.str() );
186  }
187  if( x_init_(*i_x_X) != NOT_SET_YET ) {
188  std::ostringstream omsg;
189  omsg
190  << "QPInitFixedFreeStd::initialize(...) : Error, "
191  << "Duplicate entries for i_x_fixed[i] = "
192  << (*i_x_X);
193  throw std::invalid_argument( omsg.str() );
194  }
195  x_init_(*i_x_X) = *bnd;
196  l_x_X_map_(*i_x_X) = l;
197  *i_x_X_map_itr = *i_x_X;
198  }
199  // Check that x_init was filled up properly
200  for( size_type i = 1; i <= n; ++i ) {
201  if( x_init_(i) == NOT_SET_YET ) {
202  std::ostringstream omsg;
203  omsg
204  << "QPInitFixedFreeStd::initialize(...) : Error, "
205  << "x_init(" << i << ") has not been set by"
206  " i_x_free[] or i_x_fixed[].";
207  throw std::invalid_argument( omsg.str() );
208  }
209  }
210  }
211  else {
212  for( size_type l = 1; l <= n_X; ++l, ++i_x_X, ++bnd, ++i_x_X_map_itr ) {
213  x_init_(*i_x_X) = *bnd;
214  l_x_X_map_(*i_x_X) = l;
215  *i_x_X_map_itr = *i_x_X;
216  }
217  }
218  }
219 
220  // Setup Q_R and Q_X
221  Q_R_row_i_.resize( i_x_free ? n_R : 0 );
222  Q_R_col_j_.resize( i_x_free ? n_R : 0 );
223  Q_X_row_i_.resize(n_X);
224  Q_X_col_j_.resize(n_X);
225  initialize_Q_R_Q_X(
226  n_R,n_X,i_x_free,i_x_fixed,test_setup
227  ,n_R && i_x_free ? &Q_R_row_i_[0] : NULL
228  ,n_R && i_x_free ? &Q_R_col_j_[0] : NULL
229  ,&Q_R_
230  ,n_X ? &Q_X_row_i_[0] : NULL
231  ,n_X ? &Q_X_col_j_[0] : NULL
232  ,&Q_X_
233  );
234 
235  // Setup other arguments
236  n_ = n;
237  n_R_ = n_R;
238  m_ = m;
239  g_.bind(const_cast<DVectorSlice&>(g));
240  G_ = &G;
241  A_ = A;
242  b_X_.bind(const_cast<DVectorSlice&>(b_X));
243  Ko_ = &Ko;
244  fo_.bind(const_cast<DVectorSlice&>(fo));
245  constraints_ = constraints;
246 }
247 
248 // Overridden from QP
249 
251 {
252  assert_initialized();
253  return n_;
254 }
255 
257 {
258  assert_initialized();
259  return m_;
260 }
261 
262 const DVectorSlice QPInitFixedFreeStd::g() const
263 {
264  assert_initialized();
265  return g_;
266 }
267 
268 const MatrixSymOp& QPInitFixedFreeStd::G() const
269 {
270  assert_initialized();
271  return *G_;
272 }
273 
274 const MatrixOp& QPInitFixedFreeStd::A() const
275 {
276  TEUCHOS_TEST_FOR_EXCEPT( !( A_ ) ); // ToDo: make this throw an exception
277  return *A_;
278 }
279 
281 {
282  assert_initialized();
283  return n_R_;
284 }
285 
287 {
288  assert_initialized();
289  return x_init_;
290 }
291 
293 {
294  assert_initialized();
295  return l_x_X_map_;
296 }
297 
299 {
300  assert_initialized();
301  return i_x_X_map_;
302 }
303 
304 const DVectorSlice QPInitFixedFreeStd::b_X() const
305 {
306  assert_initialized();
307  return b_X_;
308 }
309 
310 const GenPermMatrixSlice& QPInitFixedFreeStd::Q_R() const
311 {
312  assert_initialized();
313  return Q_R_;
314 }
315 
316 const GenPermMatrixSlice& QPInitFixedFreeStd::Q_X() const
317 {
318  assert_initialized();
319  return Q_X_;
320 }
321 
322 const MatrixSymOpNonsing& QPInitFixedFreeStd::Ko() const
323 {
324  assert_initialized();
325  return *Ko_;
326 }
327 
328 const DVectorSlice QPInitFixedFreeStd::fo() const
329 {
330  assert_initialized();
331  return fo_;
332 }
333 
335 {
336  assert_initialized();
337  return *constraints_;
338 }
339 
341 {
342  assert_initialized();
343  return *constraints_;
344 }
345 
346 // private member functions
347 
348 void QPInitFixedFreeStd::assert_initialized() const
349 {
350  if( !n_ )
351  throw std::logic_error( "QPInitFixedFreeStd::assert_initialized(), Error "
352  "object not initialized\n" );
353 }
354 
355 } // end namespace QPSchurPack
356 } // end namespace ConstrainedOptPack
void initialize(const DVectorSlice &g, const MatrixSymOp &G, const MatrixOp *A, size_type n_R, const size_type i_x_free[], const size_type i_x_fixed[], const EBounds bnd_fixed[], const DVectorSlice &b_X, const MatrixSymOpNonsing &Ko, const DVectorSlice &fo, Constraints *constraints, std::ostream *out=NULL, bool test_setup=false, value_type warning_tol=1e-10, value_type error_tol=1e-5, bool print_all_warnings=false)
Initialize.
size_t size_type
Represents the extra constraints in the QP to be satisfied by the schur complement QP solver QPSchur ...
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)