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_QPSchur.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 //
43 // 7/4/2002: RAB: I was able to update this class using the
44 // functions in AbstractLinAlgPack_LinAlgOpPackHack.hpp. This is wastefull in that
45 // I am creating temporaries every time any operation is performed
46 // but this was the easiest way to get things going.
47 //
48 // 7/4/2002: RAB : ToDo: In the future it would be good to create
49 // some type of temporary vector server so that I could avoid
50 // creating all of these temporaries. This will take some thought
51 // and may not be worth it for now.
52 //
53 
54 #include <assert.h>
55 
56 #include <ostream>
57 #include <iomanip>
58 #include <limits>
59 
60 #include "ConstrainedOptPack_QPSchur.hpp"
61 #include "ConstrainedOptPack_ComputeMinMult.hpp"
62 #include "AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
63 #include "AbstractLinAlgPack_SpVectorOp.hpp"
64 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
65 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
66 #include "AbstractLinAlgPack_GenPermMatrixSliceOut.hpp"
67 #include "AbstractLinAlgPack_SpVectorOut.hpp"
68 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
69 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
70 #include "AbstractLinAlgPack_VectorStdOps.hpp"
71 #include "AbstractLinAlgPack_EtaVector.hpp"
72 #include "DenseLinAlgPack_AssertOp.hpp"
73 #include "DenseLinAlgPack_DVectorClass.hpp"
74 #include "DenseLinAlgPack_DVectorClassExt.hpp"
75 #include "DenseLinAlgPack_DVectorOp.hpp"
76 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
77 #include "DenseLinAlgPack_DVectorOut.hpp"
78 #include "DenseLinAlgPack_DMatrixOut.hpp"
79 #include "Teuchos_Workspace.hpp"
80 #include "Teuchos_Assert.hpp"
81 
82 namespace LinAlgOpPack {
87 }
88 
89 namespace {
90 
91 // Some local helper functions.
92 
93 template< class T >
94 inline
95 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
96 
97 //
98 // Print a bnd as a string
99 //
100 inline
101 const char* bnd_str( ConstrainedOptPack::EBounds bnd )
102 {
103  switch(bnd) {
104  case ConstrainedOptPack::FREE:
105  return "FREE";
106  case ConstrainedOptPack::UPPER:
107  return "UPPER";
108  case ConstrainedOptPack::LOWER:
109  return "LOWER";
110  case ConstrainedOptPack::EQUALITY:
111  return "EQUALITY";
112  }
113  TEUCHOS_TEST_FOR_EXCEPT(true); // should never be executed
114  return 0;
115 }
116 
117 //
118 // print a bool
119 //
120 inline
121 const char* bool_str( bool b )
122 {
123  return b ? "true" : "false";
124 }
125 
126 //
127 // Return a std::string that has a file name, line number and
128 // error message.
129 //
130 std::string error_msg(
131  const char file_name[], const int line_num, const char err_msg[]
132  )
133 {
134  std::ostringstream omsg;
135  omsg
136  << file_name << ":" << line_num << ":" << err_msg;
137  return omsg.str();
138 }
139 
140 //
141 // Deincrement all indices less that k_remove
142 //
143 void deincrement_indices(
145  ,std::vector<DenseLinAlgPack::size_type> *indice_vector
146  ,size_t len_vector
147  )
148 {
150  typedef std::vector<DenseLinAlgPack::size_type> vec_t;
151  TEUCHOS_TEST_FOR_EXCEPT( !( len_vector <= indice_vector->size() ) );
152  for( vec_t::iterator itr = indice_vector->begin(); itr != indice_vector->begin() + len_vector; ++itr ) {
153  if( *itr > k_remove )
154  --(*itr);
155  }
156 }
157 
158 //
159 // Insert the element (r_v,c_v) into r[] and c[] sorted by r[]
160 //
161 void insert_pair_sorted(
164  ,size_t len_vector // length of the new vector
165  ,std::vector<DenseLinAlgPack::size_type> *r
166  ,std::vector<DenseLinAlgPack::size_type> *c
167  )
168 {
169  typedef std::vector<DenseLinAlgPack::size_type> rc_t;
170  TEUCHOS_TEST_FOR_EXCEPT( !( r->size() >= len_vector && c->size() >= len_vector ) );
171  // find the insertion point in r[]
172  rc_t::iterator
173  itr = std::lower_bound( r->begin(), r->begin() + len_vector-1, r_v );
174  const DenseLinAlgPack::size_type p = itr - r->begin();
175  // Shift all of the stuff out of the way to make room for the insert
176  {for( rc_t::iterator itr_last = r->begin() + len_vector-1;
177  itr_last > r->begin() + p; --itr_last )
178  {
179  *itr_last = *(itr_last-1);
180  }}
181  {for( rc_t::iterator itr_last = c->begin() + len_vector-1;
182  itr_last > c->begin() + p; --itr_last )
183  {
184  *itr_last = *(itr_last-1);
185  }}
186  // Insert the new elements
187  (*r)[p] = r_v;
188  (*c)[p] = c_v;
189 }
190 
191 //
192 // z_hat = inv(S_hat) * ( d_hat - U_hat'*vo )
193 //
194 void calc_z(
196  ,const DenseLinAlgPack::DVectorSlice &d_hat
197  ,const AbstractLinAlgPack::MatrixOp &U_hat
198  ,const DenseLinAlgPack::DVectorSlice *vo // If NULL then assumed zero
199  ,DenseLinAlgPack::DVectorSlice *z_hat
200  )
201 {
204  using Teuchos::Workspace;
206 
207  Workspace<DenseLinAlgPack::value_type> t_ws(wss,d_hat.dim());
208  DenseLinAlgPack::DVectorSlice t(&t_ws[0],t_ws.size());
209  t = d_hat;
210  if(vo)
211  Vp_StMtV( &t, -1.0, U_hat, BLAS_Cpp::trans, *vo );
212  V_InvMtV( z_hat, S_hat, BLAS_Cpp::no_trans, t );
213 }
214 
215 //
216 // v = inv(Ko) * ( fo - U_hat * z_hat )
217 //
218 void calc_v(
220  ,const DenseLinAlgPack::DVectorSlice *fo // If NULL then assumed to be zero
221  ,const AbstractLinAlgPack::MatrixOp &U_hat
222  ,const DenseLinAlgPack::DVectorSlice &z_hat // Only accessed if U_hat.cols() > 0
223  ,DenseLinAlgPack::DVectorSlice *v
224  )
225 {
226  using DenseLinAlgPack::norm_inf;
229  using Teuchos::Workspace;
231 
232  Workspace<DenseLinAlgPack::value_type> t_ws(wss,v->dim());
233  DenseLinAlgPack::DVectorSlice t(&t_ws[0],t_ws.size());
234  if(fo) {
235  t = *fo;
236  }
237  else {
238  t = 0.0;
239  }
240  if( U_hat.cols() )
241  Vp_StMtV( &t, -1.0, U_hat, BLAS_Cpp::no_trans, z_hat );
242  if( norm_inf(t) > 0.0 )
243  V_InvMtV( v, Ko, BLAS_Cpp::no_trans, t );
244  else
245  *v = 0.0;
246 }
247 
248 //
249 // mu_D_hat =
250 // - Q_XD_hat' * g
251 // - Q_XD_hat' * G * x
252 // - Q_XD_hat' * A * v(n_R+1:n_R+m)
253 // - Q_XD_hat' * A_bar * P_plus_hat * z_hat
254 //
255 void calc_mu_D(
257  ,const DenseLinAlgPack::DVectorSlice &x
258  ,const DenseLinAlgPack::DVectorSlice &v
259  ,DenseLinAlgPack::DVectorSlice *mu_D
260  )
261 {
262  using BLAS_Cpp::no_trans;
263  using BLAS_Cpp::trans;
264  using LinAlgOpPack::V_MtV;
265  using LinAlgOpPack::V_StMtV;
267  using AbstractLinAlgPack::V_MtV;
268  using AbstractLinAlgPack::Vp_MtV;
269  using Teuchos::Workspace;
271 
273  &qp = act_set.qp();
275  n = qp.n(),
276  n_R = qp.n_R(),
277  m = qp.m();
278 
279  const AbstractLinAlgPack::GenPermMatrixSlice &Q_XD_hat = act_set.Q_XD_hat();
280  const DenseLinAlgPack::DVectorSlice g = qp.g();
281  const AbstractLinAlgPack::MatrixSymOp &G = qp.G();
282  // mu_D_hat = - Q_XD_hat' * g
283  V_StMtV( mu_D, -1.0, Q_XD_hat, trans, g );
284  // mu_D_hat += - Q_XD_hat' * G * x
285  Vp_StPtMtV( mu_D, -1.0, Q_XD_hat, trans, G, no_trans, x );
286  // mu_D_hat += - Q_XD_hat' * A * v(n_R+1:n_R+m)
287  if( m ) {
288  Vp_StPtMtV( mu_D, -1.0, Q_XD_hat, trans, qp.A(), no_trans, v(n_R+1,n_R+m) );
289  }
290  // p_mu_D_hat += - Q_XD_hat' * A_bar * P_plus_hat * z_hat
291  if( act_set.q_plus_hat() && act_set.q_hat() ) {
292  const DenseLinAlgPack::DVectorSlice z_hat = act_set.z_hat();
293  AbstractLinAlgPack::SpVector P_plus_hat_z_hat;
294  V_MtV( &P_plus_hat_z_hat, act_set.P_plus_hat(), no_trans, z_hat );
295  Vp_StPtMtV( mu_D, -1.0, Q_XD_hat, trans
296  , qp.constraints().A_bar(), no_trans, P_plus_hat_z_hat() );
297  }
298 }
299 
300 //
301 // p_mu_D_hat =
302 // - Q_XD_hat' * G * Q_R * p_v(1:n_R)
303 // - Q_XD_hat' * G * P_XF_hat * p_z_hat
304 // - Q_XD_hat' * A * p_v(n_R+1:n_R+m)
305 // - Q_XD_hat' * A_bar * (P_plus_hat * p_z_hat + e(ja))
306 //
307 void calc_p_mu_D(
309  ,const DenseLinAlgPack::DVectorSlice &p_v
310  ,const DenseLinAlgPack::DVectorSlice &p_z_hat
311  ,const DenseLinAlgPack::size_type *ja // If != NULL then we will include the term e(ja)
312  ,DenseLinAlgPack::DVectorSlice *p_mu_D
313  )
314 {
315  using BLAS_Cpp::no_trans;
316  using BLAS_Cpp::trans;
317  using AbstractLinAlgPack::V_MtV;
318  using AbstractLinAlgPack::Vp_MtV;
320  using Teuchos::Workspace;
322 
324  &qp = act_set.qp();
326  &constraints = qp.constraints();
328  n = qp.n(),
329  n_R = qp.n_R(),
330  m = qp.m();
331 
332  const AbstractLinAlgPack::GenPermMatrixSlice &Q_XD_hat = act_set.Q_XD_hat();
333  const AbstractLinAlgPack::MatrixSymOp &G = qp.G();
334  // p_mu_D_hat = - Q_XD_hat' * G * Q_R * p_v(1:n_R)
335  {
336  AbstractLinAlgPack::SpVector Q_R_p_v1;
337  V_MtV( &Q_R_p_v1, qp.Q_R(), no_trans, p_v(1,n_R) );
338  Vp_StPtMtV( p_mu_D, -1.0, Q_XD_hat, trans, G, no_trans, Q_R_p_v1(), 0.0 );
339  }
340  // p_mu_D_hat += - Q_XD_hat' * G * P_XF_hat * p_z_hat
341  if( act_set.q_F_hat() ) {
342  AbstractLinAlgPack::SpVector P_XF_hat_p_z_hat;
343  V_MtV( &P_XF_hat_p_z_hat, act_set.P_XF_hat(), no_trans, p_z_hat );
344  Vp_StPtMtV( p_mu_D, -1.0, Q_XD_hat, trans, G, no_trans, P_XF_hat_p_z_hat() );
345  }
346  // p_mu_D_hat += - Q_XD_hat' * A * p_v(n_R+1:n_R+m)
347  if( m ) {
348  Vp_StPtMtV( p_mu_D, -1.0, Q_XD_hat, trans, qp.A(), no_trans, p_v(n_R+1,n_R+m) );
349  }
350  // p_mu_D_hat += - Q_XD_hat' * A_bar * ( P_plus_hat * p_z_hat + e(ja) )
351  if( act_set.q_plus_hat() || ja ) {
352  AbstractLinAlgPack::SpVector p_lambda_bar(
353  n+constraints.m_breve(), act_set.q_plus_hat() + (ja ? 1 : 0) );
354  if( act_set.q_plus_hat() ) // p_lambda_bar = P_plus_hat * p_z_hat
355  Vp_MtV( &p_lambda_bar, act_set.P_plus_hat(), no_trans, p_z_hat );
356  if( ja ) // p_lambda_bar += e(ja) (non-duplicate indices?)
357  p_lambda_bar.insert_element(AbstractLinAlgPack::SpVector::element_type(*ja,1.0));
358  Vp_StPtMtV( p_mu_D, -1.0, Q_XD_hat, trans
359  , constraints.A_bar(), no_trans, p_lambda_bar() );
360  }
361 }
362 
363 //
364 // Calculate the residual of the augmented KKT system:
365 //
366 // [ ro ] = [ Ko U_hat ] [ v ] + [ ao * bo ]
367 // [ ra ] [ U_hat' V_hat ] [ z_hat ] [ aa * ba ]
368 //
369 // Expanding this out we have:
370 //
371 // ro = Ko * v + U_hat * z_hat + ao * bo
372 //
373 // = [ Q_R'*G*Q_R Q_R'*A ] * [ x_R ] + [ Q_R'*G*P_XF_hat + Q_R'*A_bar*P_plus_hat ] * z_hat + [ ao*boR ]
374 // [ A'*Q_R ] [ lambda ] [ A'*P_XF_hat ] [ ao*bom ]
375 //
376 // = [ Q_R'*G*(Q_R*x_R + P_XF_hat*z_hat) + Q_R'*A*lambda + Q_R'*A_bar*P_plus_hat*z_hat + ao*boR ]
377 // [ A'*(Q_R*x_R + P_XF_hat*z_hat) + ao*bom ]
378 //
379 // ra = [ U_hat' * v + V_hat * z_hat + aa*ba
380 //
381 // = [ P_XF_hat'*G*Q_R + P_plus_hat'*A_bar'*Q_R , P_XF_hat'*A ] * [ x_R ; lambda ]
382 // + [ P_XF_hat'*G*P_XF_hat + P_XF_hat'*A_bar*P_plus_hat + P_plus_hat'*A_bar'*P_XF_hat
383 // + P_F_tilde'*P_C_hat + P_C_hat'*P_F_tilde ] * z_hat + aa*ba
384 //
385 // = P_XF_hat'*G*(Q_R*x_R + P_XF_hat*z_hat) + P_plus_hat'*A_bar'*(Q_R*x_R + P_XF_hat*z_hat)
386 // + P_XF_hat'*A*lambda + P_XF_hat'*A_bar*P_plus_hat*z_hat
387 // + (P_F_tilde'*P_C_hat + P_C_hat'*P_F_tilde)*z_hat + aa*ba
388 //
389 // Given the QP matrices G, A, and A_bar and the active set mapping matrices Q_R, P_XF_hat,
390 // P_plus_hat and P_FC_hat = P_F_tilde'*P_C_hat, we can compute the residual efficiently
391 // as follows:
392 //
393 // x_free = Q_R*x_R + P_XF_hat*z_hat (sparse)
394 //
395 // lambda_bar = P_plus_hat*z_hat (sparse)
396 //
397 // t1 = G*x_free
398 //
399 // t2 = A*lambda
400 //
401 // t3 = A_bar*lambda_bar
402 //
403 // roR = Q_R'*t1 + Q_R'*t2 + Q_R'*t3 + ao*boR
404 //
405 // rom = A'*t1 + ao*bom
406 //
407 // ra = P_XF_hat'*t1 + P_plus_hat'*A_bar'*x_free + P_XF_hat'*t2 + P_XF_hat'*t3
408 // + (P_FC_hat + P_FC_hat')*z_hat + aa*ba
409 //
410 // On output we will have set:
411 //
412 // roR_scaling = ||Q_R'*t1||inf + ||Q_R'*t2||inf + ||Q_R'*t3||inf + ||ao*boR||inf
413 //
414 // rom_scaling = ||A'*t1||inf + ||ao*bom||inf
415 //
416 // ra_scaling = ||P_XF_hat'*t1||inf + ||P_plus_hat'*A_bar'*t1||inf + ||P_XF_hat'*t2||inf
417 // + ||P_XF_hat'*t3||inf + ||(P_FC_hat + P_FC_hat')*z_hat|| + ||aa*ba||inf
418 //
419 // Note: In the future we could make this a little more efficent by combining (Q_R + P_XF_hat) into
420 // an single permulation matrix and then we could leave out the terms for the variables initially
421 // fixed and still fixed rather than computing the terms then throwing them away.
422 //
423 // Also, in the future, this could really be implemented for extended precision data types but
424 // it would require some real work!
425 //
426 template<class val_type>
427 void calc_resid(
429  ,const DenseLinAlgPack::DVectorSlice &v
430  ,const DenseLinAlgPack::DVectorSlice &z_hat // Only accessed if q_hat > 0
431  ,const DenseLinAlgPack::value_type ao // Only accessed if bo != NULL
432  ,const DenseLinAlgPack::DVectorSlice *bo // If NULL then considered 0
434  ,DenseLinAlgPack::value_type *roR_scaling
435  ,DenseLinAlgPack::value_type *rom_scaling // Only set if m > 0
436  ,const DenseLinAlgPack::value_type aa // Only accessed if q_hat > 0
437  ,const DenseLinAlgPack::DVectorSlice *ba // If NULL then considered 0, Only accessed if q_hat > 0
438  ,DenseLinAlgPack::VectorSliceTmpl<val_type> *ra // Only set if q_hat > 0
439  ,DenseLinAlgPack::value_type *ra_scaling // Only set if q_hat > 0
440  )
441 {
442  using BLAS_Cpp::no_trans;
443  using BLAS_Cpp::trans;
444  using DenseLinAlgPack::norm_inf;
445  using DenseLinAlgPack::DVectorSlice;
447  using AbstractLinAlgPack::V_MtV;
448  using AbstractLinAlgPack::SpVector;
450  using LinAlgOpPack::Vp_V;
451  using LinAlgOpPack::V_StV;
452  using LinAlgOpPack::V_MtV;
453  using LinAlgOpPack::Vp_MtV;
456  namespace COP = ConstrainedOptPack;
457  using Teuchos::Workspace;
459 
460  const COP::QPSchurPack::QP
461  &qp = act_set.qp();
462  const COP::QPSchurPack::Constraints
463  &constraints = qp.constraints();
464  const GenPermMatrixSlice
465  &Q_R = qp.Q_R(),
466  &P_XF_hat = act_set.P_XF_hat(),
467  &P_plus_hat = act_set.P_plus_hat();
469  n = qp.n(),
470  n_R = qp.n_R(),
471  m = qp.m(),
472  m_breve = constraints.m_breve(),
473  q_hat = act_set.q_hat(),
474  q_F_hat = act_set.q_F_hat(),
475  q_C_hat = act_set.q_C_hat(),
476  q_plus_hat = act_set.q_plus_hat();
477 
478  const DVectorSlice
479  x_R = v(1,n_R),
480  lambda = ( m ? v(n_R+1,n_R+m): DVectorSlice() ),
481  boR = ( bo ? (*bo)(1,n_R) : DVectorSlice() ),
482  bom = ( bo && m ? (*bo)(n_R+1,n_R+m) : DVectorSlice() );
483 
485  roR = (*ro)(1,n_R),
486  rom = ( m ? (*ro)(n_R+1,n_R+m) : DenseLinAlgPack::VectorSliceTmpl<val_type>() );
487 
488  Workspace<DenseLinAlgPack::value_type>
489  x_free_ws(wss,n);
490  DenseLinAlgPack::DVectorSlice
491  x_free(&x_free_ws[0],x_free_ws.size());
492  Workspace<val_type>
493  t1_ws(wss,n),
494  t2_ws(wss,n),
495  t3_ws(wss,n),
496  tR_ws(wss,n_R),
497  tm_ws(wss,m),
498  ta_ws(wss,q_hat);
500  t1(&t1_ws[0],t1_ws.size()),
501  t2(&t2_ws[0],t2_ws.size()),
502  t3(&t3_ws[0],t3_ws.size()),
503  tR(&tR_ws[0],tR_ws.size()),
504  tm(tm_ws.size()?&tm_ws[0]:NULL,tm_ws.size()),
505  ta(ta_ws.size()?&ta_ws[0]:NULL,ta_ws.size());
506 
507  *roR_scaling = 0.0;
508  if( m )
509  *rom_scaling = 0.0;
510  if( q_hat )
511  *ra_scaling = 0.0;
512 
513  // x_free = Q_R*x_R + P_XF_hat*z_hat (dense for now)
514  LinAlgOpPack::V_MtV( &x_free, Q_R, no_trans, x_R );
515  if( q_F_hat )
516  LinAlgOpPack::Vp_MtV( &x_free, P_XF_hat, no_trans, z_hat );
517  // lambda_bar = P_plus_hat*z_hat
518  SpVector lambda_bar;
519  if( q_plus_hat )
520  V_MtV( &lambda_bar, P_plus_hat, no_trans, z_hat );
521  // t1 = G*x_free
522  Vp_StMtV( &t1, 1.0, qp.G(), no_trans, x_free, 0.0 );
523  // t2 = A*lambda
524  if( m )
525  Vp_StMtV( &t2, 1.0, qp.A(), no_trans, lambda, 0.0 );
526  // t3 = A_bar*lambda_bar
527  if( q_plus_hat )
528  Vp_StMtV( &t3, 1.0, constraints.A_bar(), no_trans, lambda_bar(), 0.0 );
529  // roR = Q_R'*t1 + Q_R'*t2 + Q_R'*t3 + ao*boR
530  LinAlgOpPack::V_MtV( &tR, Q_R, trans, t1 ); // roR = Q_R'*t1
531  *roR_scaling += norm_inf(tR);
532  roR = tR;
533  if( m ) { // roR += Q_R'*t2
534  LinAlgOpPack::V_MtV( &tR, Q_R, trans, t2 );
535  *roR_scaling += norm_inf(tR);
536  LinAlgOpPack::Vp_V(&roR,tR);
537  }
538  if( q_plus_hat ) { // roR += Q_R'*t3
539  LinAlgOpPack::V_MtV( &tR, Q_R, trans, t3 );
540  *roR_scaling += norm_inf(tR);
541  LinAlgOpPack::Vp_V(&roR,tR);
542  }
543  if( bo ) {
544  LinAlgOpPack::Vp_StV( &roR, ao, boR ); // roR += ao*boR
545  *roR_scaling += std::fabs(ao) * norm_inf(boR);
546  }
547  // rom = A'*t1 + ao*bom
548  if( m ) {
549  Vp_StMtV( &tm, 1.0, qp.A(), trans, t1, 0.0 ); // A'*t1
550  *rom_scaling += norm_inf(tm);
551  LinAlgOpPack::Vp_V(&rom,tm);
552  if(bo) {
553  LinAlgOpPack::Vp_StV( &rom, ao, bom ); // rom += ao*bom
554  *rom_scaling += std::fabs(ao)*norm_inf(bom);
555  }
556  }
557  // ra = P_XF_hat'*t1 + P_plus_hat'*A_bar'*x_free + P_XF_hat'*t2 + P_XF_hat'*t3
558  // +(P_FC_hat + P_FC_hat')*z_hat + aa*ba
559  if( q_hat ) {
560  if(ba) { // ra = aa*ba
561  V_StV( ra, aa, *ba );
562  *ra_scaling += std::fabs(aa) * norm_inf(*ba);
563  }
564  else {
565  *ra = 0.0;
566  }
567  if( q_F_hat ) { // ra += P_XF_hat'*t1
568  LinAlgOpPack::V_MtV( &ta, P_XF_hat, trans, t1 );
569  *ra_scaling += norm_inf(ta);
570  LinAlgOpPack::Vp_V(ra,ta);
571  }
572  if( q_plus_hat ) { // ra += P_plus_hat'*A_bar'*x_free
573  Vp_StPtMtV( &ta, 1.0, P_plus_hat, trans, constraints.A_bar(), trans, x_free, 0.0 );
574  *ra_scaling += norm_inf(ta);
575  LinAlgOpPack::Vp_V(ra,ta);
576  }
577  if( q_F_hat && m ) { // ra += P_XF_hat'*t2
578  LinAlgOpPack::V_MtV( &ta, P_XF_hat, trans, t2 );
579  *ra_scaling += norm_inf(ta);
580  LinAlgOpPack::Vp_V(ra,ta);
581  }
582  if( q_F_hat && q_plus_hat ) { // ra += P_XF_hat'*t3
583  LinAlgOpPack::V_MtV( &ta, P_XF_hat, trans, t3 );
584  *ra_scaling += norm_inf(ta);
585  LinAlgOpPack::Vp_V(ra,ta);
586  }
587  if( q_C_hat ) { // ra += (P_FC_hat + P_FC_hat')*z_hat
588  const GenPermMatrixSlice
589  &P_FC_hat = act_set.P_FC_hat();
590  ta = 0.0;
591  for( GenPermMatrixSlice::const_iterator itr = P_FC_hat.begin(); itr != P_FC_hat.end(); ++itr ) {
592  ta(itr->row_i()) = z_hat(itr->col_j());
593  ta(itr->col_j()) = z_hat(itr->row_i());
594  }
595  *ra_scaling += norm_inf(ta);
596  LinAlgOpPack::Vp_V(ra,ta);
597  }
598  }
599 }
600 
601 //
602 // Correct a nearly degenerate Lagrange multiplier
603 //
604 // If < 0 is returned it means that the multiplier could not
605 // be corrected and this should be considered to be dual infeasible
606 // In this case the error output is sent to *out if print_dual_infeas
607 // == true.
608 //
609 // If 0 is returned then the multiplier was near degenerate and
610 // was corrected. In this case a warning message is printed to
611 // *out.
612 //
613 // If > 0 is returned then the multiplier's sign
614 // is just fine and no corrections were needed (no output).
615 //
616 int correct_dual_infeas(
617  const DenseLinAlgPack::size_type j // for output info only
618  ,const ConstrainedOptPack::EBounds bnd_j
619  ,const DenseLinAlgPack::value_type t_P // (> 0) full step length
620  ,const DenseLinAlgPack::value_type scale // (> 0) scaling value
621  ,const DenseLinAlgPack::value_type dual_infeas_tol
622  ,const DenseLinAlgPack::value_type degen_mult_val
623  ,std::ostream *out // Can be NULL
625  ,const bool print_dual_infeas
626  ,const char nu_j_n[] // Name of nu_j
627  ,DenseLinAlgPack::value_type *nu_j // required
628  ,DenseLinAlgPack::value_type *scaled_viol // = scale*nu_j*(bnd_j==UPPER ? 1.0: -1.0 ) (after output)
629  ,const char p_nu_j_n[] = NULL // Name of p_nu_j (can be NULL if p_nu_j==NULL)
630  ,DenseLinAlgPack::value_type *p_nu_j = NULL // optional (can be NULL)
631  ,const char nu_j_plus_n[] = NULL // Name of nu_j_plus (can be NULL if p_nu_j==NULL)
632  ,DenseLinAlgPack::value_type *nu_j_plus = NULL // optional (can be NULL)
633  )
634 {
635  typedef DenseLinAlgPack::value_type value_type;
636  namespace COP = ConstrainedOptPack;
637 
638  value_type nu_j_max = (*scaled_viol) = scale * (*nu_j) * (bnd_j == COP::UPPER ? +1.0 : -1.0);
639  if( nu_j_max > 0.0 || bnd_j == COP::EQUALITY ) // Leave any multiplier value with the correct sign alone!
640  return +1; // No correction needed
641  // See if we need to correct the multiplier
642  nu_j_max = std::fabs(nu_j_max);
643  if( nu_j_max < dual_infeas_tol ) {
644  // This is a near degenerate multiplier so adjust it
645  value_type degen_val = degen_mult_val * ( bnd_j == COP::UPPER ? +1.0 : -1.0 );
646  if( out && (int)olevel >= (int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
647  *out
648  << "\nWarning, the constriant a(" << j << ") currently at its "
649  << (bnd_j == COP::UPPER ? "UPPER" : "LOWER") << " bound"
650  << " has the wrong Lagrange multiplier value but\n"
651  << "scale*|"<<nu_j_n<<"| = " << scale << " * |" << (*nu_j)
652  << "| = " << nu_j_max << " < dual_infeas_tol = " << dual_infeas_tol
653  << "\nTherefore, this is considered a degenerate constraint and this "
654  << "multiplier is set to " << degen_val << std::endl;
655  }
656  if(p_nu_j) {
657  nu_j_max += std::fabs( t_P * (*p_nu_j) ) * scale;
658  if( nu_j_max < dual_infeas_tol ) {
659  // The full step is also degenerate so adjust it also
660  if( out && (int)olevel >= (int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
661  *out
662  << "Also, the maximum full step scale*(|"<<nu_j_n<<"|+|t_P*"<<p_nu_j_n<<"|) = "
663  << scale << " * (|" << (*nu_j) << "| + |" << t_P << " * " << (*p_nu_j) << "|) = "
664  << nu_j_max << " < dual_infeas_tol = " << dual_infeas_tol
665  << "\nTherefore, this is considered degenerate and therefore "
666  << "seting "<<p_nu_j_n<<" = 0";
667  if(nu_j_plus)
668  *out << " and "<< nu_j_plus_n <<" = " << degen_val;
669  *out << std::endl;
670  }
671  *p_nu_j = 0.0; // Don't let it limit the step length
672  if(nu_j_plus) {
673  *nu_j_plus = degen_val;
674  }
675  }
676  }
677  *nu_j = degen_val; // Now set it
678  *scaled_viol = scale * (*nu_j) * (bnd_j == COP::UPPER ? +1.0 : -1.0); // Not violated!
679  return 0;
680  }
681  else {
682  if( print_dual_infeas && out && (int)olevel >= (int)COP::QPSchur::OUTPUT_BASIC_INFO ) {
683  *out
684  << "\nError, the constriant a(" << j << ") currently at its "
685  << (bnd_j == COP::UPPER ? "UPPER" : "LOWER") << " bound"
686  << " has the wrong Lagrange multiplier value and\n"
687  << "scale*|"<<nu_j_n<<"| = " << scale << " * |" << (*nu_j)
688  << "| = " << nu_j_max << " > dual_infeas_tol = " << dual_infeas_tol
689  << "\nThis is an indication of instability in the calculations.\n"
690  << "The QP algorithm is terminated!\n";
691  }
692  return -1;
693  }
694  return 0; // Will never be executed!
695 }
696 
697 //
698 // Calculate the QP objective if it has not been already
699 //
700 // qp_grad_norm_inf = ||g + G*x||inf
701 //
702 // On input, if qp_grad_norm_inf != 0, then it will be assumed
703 // that this value has already been computed and the computation will
704 // be skipped.
705 //
706 void calc_obj_grad_norm_inf(
708  ,const DenseLinAlgPack::DVectorSlice &x
709  ,DenseLinAlgPack::value_type *qp_grad_norm_inf
710  )
711 {
712  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this?
713 }
714 
715 } // end namespace
716 
717 namespace ConstrainedOptPack {
718 
719 // public member functions for QPSchur::U_hat_t
720 
722  :G_(NULL)
723  ,A_(NULL)
724  ,A_bar_(NULL)
725  ,Q_R_(NULL)
726  ,P_XF_hat_(NULL)
727  ,P_plus_hat_(NULL)
728 {}
729 
731  const MatrixSymOp *G
732  ,const MatrixOp *A
733  ,const MatrixOp *A_bar
734  ,const GenPermMatrixSlice *Q_R
735  ,const GenPermMatrixSlice *P_XF_hat
736  ,const GenPermMatrixSlice *P_plus_hat
737  )
738 {
739  G_ = G;
740  A_ = A;
741  A_bar_ = A_bar;
742  Q_R_ = Q_R;
743  P_XF_hat_ = P_XF_hat;
744  P_plus_hat_ = P_plus_hat;
745 }
746 
747 size_type QPSchur::U_hat_t::rows() const
748 {
749  return Q_R_->cols() + ( A_ ? A_->cols() : 0 );
750 }
751 
752 size_type QPSchur::U_hat_t::cols() const
753 {
754  return P_plus_hat_->cols();
755 }
756 
757 /* 10/25/00: I don't think we need this function!
758 void QPSchur::U_hat_t::Mp_StM(DMatrixSlice* C, value_type a
759  , BLAS_Cpp::Transp M_trans ) const
760 {
761  using BLAS_Cpp::no_trans;
762  using BLAS_Cpp::trans;
763  using LinAlgOpPack::Mp_StMtP;
764  using LinAlgOpPack::Mp_StPtMtP;
765 
766  // C += a * op(U_hat)
767 
768  LinAlgOpPack::Mp_M_assert_sizes( C->rows(), C->cols(), no_trans
769  , rows(), cols(), M_trans );
770 
771  const size_type
772  n_R = Q_R_->cols(),
773  m = A() ? A()->cols() : 0;
774 
775  if( M_trans == no_trans ) {
776  //
777  // C += a * op(U_hat)
778  //
779  // = a * [ Q_R' * G * P_XF_hat + Q_R' * A_bar * P_plus_hat ]
780  // [ A' * P_XF_hat ]
781  //
782  // C1 += a * Q_R' * G * P_XF_hat + a * Q_R' * A_bar * P_plus_hat
783  //
784  // C2 += a * A' * P_XF_hat
785  //
786  DMatrixSlice
787  C1 = (*C)(1,n_R,1,C->cols()),
788  C2 = m ? (*C)(n_R+1,n_R+m,1,C->cols()) : DMatrixSlice();
789  // C1 += a * Q_R' * G * P_XF_hat
790  if( P_XF_hat().nz() )
791  Mp_StPtMtP( &C1, a, Q_R(), trans, G(), no_trans, P_XF_hat(), no_trans );
792  // C1 += a * Q_R' * A_bar * P_plus_hat
793  if( P_plus_hat().nz() )
794  Mp_StPtMtP( &C1, a, Q_R(), trans, A_bar(), no_trans, P_plus_hat(), no_trans );
795  // C2 += a * A' * P_XF_hat
796  if(m && P_XF_hat().nz())
797  Mp_StMtP( &C2, a, *A(), trans, P_plus_hat(), no_trans );
798  }
799  else {
800  TEUCHOS_TEST_FOR_EXCEPT(true); // Implement this!
801  }
802 }
803 */
804 
806  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
807  ,const DVectorSlice& x, value_type b
808  ) const
809 {
810  using BLAS_Cpp::no_trans;
811  using BLAS_Cpp::trans;
812  using DenseLinAlgPack::Vt_S;
813  using AbstractLinAlgPack::V_MtV;
815  using LinAlgOpPack::V_MtV;
818  using Teuchos::Workspace;
820 
821  LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),rows(),cols(),M_trans,x.dim());
822 
823  //
824  // U_hat = [ Q_R' * G * P_XF_hat + Q_R' * A_bar * P_plus_hat ]
825  // [ A' * P_XF_hat ]
826  //
827 
828  const size_type
829  n_R = Q_R_->cols(),
830  m = A() ? A()->cols() : 0;
831 
832  if( M_trans == BLAS_Cpp::no_trans ) {
833  //
834  // y = b*y + a * U_hat * x
835  //
836  // = b*y + a * [ Q_R' * G * P_XF_hat + Q_R' * A_bar * P_plus_hat ] * x
837  // [ A' * P_XF_hat ]
838  //
839  // =>
840  //
841  // y1 = b * y1 + a * Q_R' * G * P_XF_hat * x + a * Q_R' * A_bar * P_plus_hat * x
842  // y2 = b * y2 + a * A' * P_XF_hat * x
843  //
844  DVectorSlice
845  y1 = (*y)(1,n_R),
846  y2 = m ? (*y)(n_R+1,n_R+m) : DVectorSlice();
847  SpVector
848  P_XF_hat_x,
849  P_plus_hat_x;
850  // P_XF_hat_x = P_XF_hat * x
851  if( P_XF_hat().nz() )
852  V_MtV( &P_XF_hat_x, P_XF_hat(), no_trans, x );
853  // P_plus_hat_x = P_plus_hat * x
854  if(P_plus_hat().nz())
855  V_MtV( &P_plus_hat_x, P_plus_hat(), no_trans, x );
856  // y1 = b * y1
857  if(b==0.0) y1=0.0;
858  else if(b!=1.0) Vt_S(&y1,b);
859  // y1 += a * Q_R' * G * P_XF_hat_x
860  if(P_XF_hat().nz())
861  Vp_StPtMtV( &y1, a, Q_R(), trans, G(), no_trans, P_XF_hat_x() );
862  // y1 += a * Q_R' * A_bar * P_plus_hat_x
863  if(P_plus_hat().nz())
864  Vp_StPtMtV( &y1, a, Q_R(), trans, A_bar(), no_trans, P_plus_hat_x() );
865  if(m) {
866  // y2 = b * y2
867  if(b==0.0) y2=0.0;
868  else if(b!=1.0) Vt_S(&y2,b);
869  // y2 += a * A' * P_XF_hat_x
870  if( P_XF_hat().nz() )
871  Vp_StMtV( &y2, a, *A(), trans, P_XF_hat_x() );
872  }
873  }
874  else if( M_trans == BLAS_Cpp::trans ) {
875  //
876  // y = b*y + a * U_hat' * x
877  //
878  // = b*y + a * [ P_XF_hat' * G * Q_R + P_plus_hat' * A_bar' * Q_R, P_XF_hat' * A ] * [ x1 ]
879  // [ x2 ]
880  // =>
881  //
882  // y = b * y + a * P_XF_hat' * G * Q_R * x1 + a * P_plus_hat' * A_bar' * Q_R * x1
883  // + a * P_XF_hat' * A * x2
884  //
885  const DVectorSlice
886  x1 = x(1,n_R),
887  x2 = m ? x(n_R+1,n_R+m) : DVectorSlice();
888  SpVector
889  Q_R_x1;
890  // Q_R_x1 = Q_R * x1
891  V_MtV( &Q_R_x1, Q_R(), no_trans, x1 );
892  // y = b*y
893  if(b==0.0) *y = 0.0;
894  else if(b!=1.0) Vt_S( y, b );
895  // y += a * P_XF_hat' * G * Q_R_x1
896  if(P_XF_hat().nz())
897  Vp_StPtMtV( y, a, P_XF_hat(), trans, G(), no_trans, Q_R_x1() );
898  // y += a * P_plus_hat' * A_bar' * Q_R_x1
899  if(P_plus_hat().nz())
900  Vp_StPtMtV( y, a, P_plus_hat(), trans, A_bar(), trans, Q_R_x1() );
901  // y += a * P_XF_hat' * A * x2
902  if( m && P_XF_hat().nz() )
903  Vp_StPtMtV( y, a, P_XF_hat(), trans, *A(), no_trans, x2 );
904  }
905  else {
906  TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid value for M_trans
907  }
908 }
909 
911  DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
912  ,const SpVectorSlice& x, value_type b
913  ) const
914 {
915 // // Uncomment to use the default version
916 // MatrixOp::Vp_StMtV(y,a,M_trans,x,b); return;
917 
918  using BLAS_Cpp::no_trans;
919  using BLAS_Cpp::trans;
920  using DenseLinAlgPack::Vt_S;
921  using LinAlgOpPack::V_MtV;
922  using AbstractLinAlgPack::V_MtV;
924  using LinAlgOpPack::V_MtV;
927  using Teuchos::Workspace;
929 
930  LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),rows(),cols(),M_trans,x.dim());
931 
932  //
933  // U_hat = [ Q_R' * G * P_XF_hat + Q_R' * A_bar * P_plus_hat ]
934  // [ A' * P_XF_hat ]
935  //
936 
937  const size_type
938  n_R = Q_R_->cols(),
939  m = A() ? A()->cols() : 0;
940 
941  if( M_trans == BLAS_Cpp::no_trans ) {
942  //
943  // y = b*y + a * U_hat * x
944  //
945  // = b*y + a * [ Q_R' * G * P_XF_hat + Q_R' * A_bar * P_plus_hat ] * x
946  // [ A' * P_XF_hat ]
947  //
948  // =>
949  //
950  // y1 = b * y1 + a * Q_R' * G * P_XF_hat * x + a * Q_R' * A_bar * P_plus_hat * x
951  // y2 = b * y2 + a * A' * P_XF_hat * x
952  //
953  DVectorSlice
954  y1 = (*y)(1,n_R),
955  y2 = m ? (*y)(n_R+1,n_R+m) : DVectorSlice();
956  SpVector
957  P_XF_hat_x,
958  P_plus_hat_x;
959  // P_XF_hat_x = P_XF_hat * x
960  if( P_XF_hat().nz() )
961  V_MtV( &P_XF_hat_x, P_XF_hat(), no_trans, x );
962  // P_plus_hat_x = P_plus_hat * x
963  if(P_plus_hat().nz())
964  V_MtV( &P_plus_hat_x, P_plus_hat(), no_trans, x );
965  // y1 = b * y1
966  if(b==0.0) y1=0.0;
967  else if(b!=1.0) Vt_S(&y1,b);
968  // y1 += a * Q_R' * G * P_XF_hat_x
969  if(P_XF_hat().nz())
970  Vp_StPtMtV( &y1, a, Q_R(), trans, G(), no_trans, P_XF_hat_x() );
971  // y1 += a * Q_R' * A_bar * P_plus_hat_x
972  if(P_plus_hat().nz())
973  Vp_StPtMtV( &y1, a, Q_R(), trans, A_bar(), no_trans, P_plus_hat_x() );
974  if(m) {
975  // y2 = b * y2
976  if(b==0.0) y2=0.0;
977  else if(b!=1.0) Vt_S(&y2,b);
978  // y2 += a * A' * P_XF_hat_x
979  if(P_XF_hat().nz())
980  Vp_StMtV( &y2, a, *A(), trans, P_XF_hat_x() );
981  }
982  }
983  else if( M_trans == BLAS_Cpp::trans ) {
984  //
985  // y = b*y + a * U_hat' * x
986  //
987  // = b*y + a * [ P_XF_hat' * G * Q_R + P_plus_hat' * A_bar' * Q_R, P_XF_hat' * A ] * [ x1 ]
988  // [ x2 ]
989  // =>
990  //
991  // y = b * y + a * P_XF_hat' * G * Q_R * x1 + a * P_plus_hat' * A_bar' * Q_R * x1
992  // + a * P_XF_hat' * A * x2
993  //
994  const SpVectorSlice
995  x1 = x(1,n_R),
996  x2 = m ? x(n_R+1,n_R+m) : SpVectorSlice(NULL,0,0,0);
997  SpVector
998  Q_R_x1;
999  // Q_R_x1 = Q_R * x1
1000  V_MtV( &Q_R_x1, Q_R(), no_trans, x1 );
1001  // y = b*y
1002  if(b ==0.0) *y = 0.0;
1003  else if(b!=1.0) Vt_S( y, b );
1004  // y += a * P_XF_hat' * G * Q_R_x1
1005  if(P_XF_hat().nz())
1006  Vp_StPtMtV( y, a, P_XF_hat(), trans, G(), no_trans, Q_R_x1() );
1007  // y += a * P_plus_hat' * A_bar' * Q_R_x1
1008  if(P_plus_hat().nz())
1009  Vp_StPtMtV( y, a, P_plus_hat(), trans, A_bar(), trans, Q_R_x1() );
1010  // y += a * P_XF_hat' * A * x2
1011  if( m && P_XF_hat().nz() )
1012  Vp_StPtMtV( y, a, P_XF_hat(), trans, *A(), no_trans, x2 );
1013  }
1014  else {
1015  TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid value for M_trans
1016  }
1017 }
1018 
1019 // public member functions for QPSchur::ActiveSet
1020 
1021 QPSchur::ActiveSet::ActiveSet(
1022  const schur_comp_ptr_t& schur_comp
1023  ,MSADU::PivotTolerances pivot_tols
1024  )
1025  :schur_comp_(schur_comp)
1026  ,pivot_tols_(pivot_tols)
1027  ,initialized_(false)
1028  ,test_(false)
1029  ,qp_(NULL)
1030  ,x_init_(NULL)
1031  ,n_(0)
1032  ,n_R_(0)
1033  ,m_(0)
1034  ,q_plus_hat_(0)
1035  ,q_F_hat_(0)
1036  ,q_C_hat_(0)
1037 {}
1038 
1040  QP& qp, size_type num_act_change, const int ij_act_change[]
1041  ,const EBounds bnds[], bool test, bool salvage_init_schur_comp
1042  ,std::ostream *out, EOutputLevel output_level )
1043 {
1044  using LinAlgOpPack::V_mV;
1045  using LinAlgOpPack::V_MtV;
1046  using LinAlgOpPack::V_InvMtV;
1048  using AbstractLinAlgPack::V_MtV;
1051  using DenseLinAlgPack::sym;
1052  typedef MatrixSymAddDelUpdateable MSADU;
1053  namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
1054  using Teuchos::Workspace;
1056 
1057  const size_type
1058  n = qp.n(),
1059  n_R = qp.n_R(),
1060  n_X = n - n_R,
1061  m = qp.m();
1062  const QP::x_init_t
1063  &x_init = qp.x_init();
1064  const QP::l_x_X_map_t
1065  &l_x_X_map = qp.l_x_X_map();
1066  const QP::i_x_X_map_t
1067  &i_x_X_map = qp.i_x_X_map();
1068  const DVectorSlice
1069  b_X = qp.b_X();
1070  const DVectorSlice
1071  g = qp.g();
1072  const MatrixSymOp
1073  &G = qp.G();
1074  const QP::Constraints
1075  &constraints = qp.constraints();
1076  const size_type
1077  m_breve = constraints.m_breve();
1078 
1079  try {
1080 
1081  // Count the number of each type of change
1082  size_type
1083  q_plus_hat = 0,
1084  q_F_hat = 0,
1085  q_C_hat = 0;
1086  if( num_act_change ) {
1087  for( size_type k = 1; k <= num_act_change; ++k ) {
1088  const int ij = ij_act_change[k-1];
1089  const EBounds bnd = bnds[k-1];
1090  if( ij < 0 ) {
1091  // Initially fixed variable being freed.
1092  if( x_init(-ij) == FREE ) {
1093  std::ostringstream omsg;
1094  omsg
1095  << "QPSchur::ActiveSet::initialize(...) : Error, "
1096  << "The variable x(" << -ij << ") is not initially fixed and can not "
1097  << "be freed by ij_act_change[" << k-1 << "]\n";
1098  throw std::invalid_argument( omsg.str() );
1099  }
1100  if( x_init(-ij) == EQUALITY ) {
1101  std::ostringstream omsg;
1102  omsg
1103  << "QPSchur::ActiveSet::initialize(...) : Error, "
1104  << "The variable x(" << -ij << ") is equality fixed and therefore can not "
1105  << "be freed by ij_act_change[" << k-1 << "]\n";
1106  throw std::invalid_argument( omsg.str() );
1107  }
1108  ++q_F_hat;
1109  }
1110  else {
1111  // Constraint being added to the active-set
1112  if( ij <= n ) {
1113  // Fixing a variable to a bound
1114  EBounds x_init_bnd = x_init(ij);
1115  if( x_init_bnd == FREE ) {
1116  // initially free variable being fixed
1117  ++q_plus_hat;
1118  }
1119  else if ( x_init_bnd == EQUALITY ) {
1120  // ToDo: Throw exception
1122  }
1123  else if( x_init_bnd == bnd ) {
1124  // ToDo: Throw exception
1126  }
1127  else {
1128  // Initially fixed variable being fixed to another bound
1129  ++q_F_hat; // We must free the variable first
1130  ++q_C_hat; // Now we fix it to a different bound.
1131  }
1132  }
1133  else {
1134  // Adding a general inequality (or equality) constraint
1135  if( ij > n + m_breve ) {
1136  // ToDo: Throw exception
1138  }
1139  ++q_plus_hat;
1140  }
1141  }
1142  }
1143  }
1144 
1145  const size_type
1146  q_D_hat = (n - n_R) - q_F_hat,
1147  q_D_hat_max = n_X;
1148 
1149  // Now let's set stuff up: ij_map, constr_norm, bnds and part of d_hat
1150  const size_type
1151  q_hat = q_plus_hat + q_F_hat + q_C_hat,
1152  q_hat_max = n_X + n, // If all the initially fixed variables where freed
1153  // Then all the degrees of freedom where used up with other constraints.
1154  q_F_hat_max = n_X,
1155  q_plus_hat_max = n;
1156 
1157  ij_map_.resize(q_hat_max);
1158  constr_norm_.resize(q_hat_max);
1159  bnds_.resize(q_hat_max);
1160  d_hat_.resize(q_hat_max); // set the terms involving the bounds first.
1161 
1162  if( num_act_change ) {
1163  size_type s = 0;
1164  for( size_type k = 1; k <= num_act_change; ++k ) {
1165  const int ij = ij_act_change[k-1];
1166  const EBounds bnd = bnds[k-1];
1167  if( ij < 0 ) {
1168  // Initially fixed variable being freed.
1169  ij_map_[s] = ij;
1170  constr_norm_[s] = 1.0;
1171  bnds_[s] = FREE;
1172  d_hat_[s] = - g(-ij); // - g^X_{l^{(-)}}
1173  ++s;
1174  }
1175  else {
1176  // Constraint being added to the active-set
1177  if( ij <= n ) {
1178  // Fixing a variable to a bound
1179  EBounds x_init_bnd = x_init(ij);
1180  if( x_init_bnd == FREE ) {
1181  // initially free variable being fixed
1182  ij_map_[s] = ij;
1183  constr_norm_[s] = 1.0;
1184  bnds_[s] = bnd;
1185  d_hat_[s] = constraints.get_bnd(ij,bnd);
1186  ++s;
1187  }
1188  else {
1189  // Initially fixed variable being fixed to another bound
1190  // Free the variable first
1191  ij_map_[s] = ij;
1192  constr_norm_[s] = 1.0;
1193  bnds_[s] = FREE;
1194  d_hat_[s] = - g(ij); // - g^X_{l^{(-)}}
1195  ++s;
1196  // Now fix to a different bound
1197  ij_map_[s] = ij;
1198  constr_norm_[s] = 1.0;
1199  bnds_[s] = bnd;
1200  d_hat_[s] = constraints.get_bnd(ij,bnd) - b_X(l_x_X_map(ij));
1201  ++s;
1202  }
1203  }
1204  else {
1205  // Adding a general inequality (or equality) constraint
1206  ij_map_[s] = ij;
1207  constr_norm_[s] = 1.0; // ToDo: We need to compute this in an efficient way!
1208  bnds_[s] = bnd;
1209  d_hat_[s] = constraints.get_bnd(ij,bnd); // \bar{b}_{j^{(+)}}
1210  ++s;
1211  }
1212  }
1213  }
1214  TEUCHOS_TEST_FOR_EXCEPT( !( s == q_hat ) );
1215  }
1216 
1217  // Setup P_XF_hat_ and P_plus_hat_
1218  P_XF_hat_row_.resize(q_F_hat_max);
1219  P_XF_hat_col_.resize(q_F_hat_max);
1220  P_FC_hat_row_.resize(q_F_hat_max);
1221  P_FC_hat_col_.resize(q_F_hat_max);
1222  P_plus_hat_row_.resize(q_plus_hat_max);
1223  P_plus_hat_col_.resize(q_plus_hat_max);
1224  if(q_hat) {
1225  // See ConstrainedOptPack_QPSchur.hpp for description of P_XF_hat and P_plus_hat
1226  size_type
1227  k_XF_hat = 0, // zero based
1228  k_plus_hat = 0; // zero based
1229  ij_map_t::const_iterator
1230  ij_itr = ij_map_.begin(),
1231  ij_itr_end = ij_itr + q_hat;
1232  for( size_type s = 1; ij_itr != ij_itr_end; ++ij_itr, ++s ) {
1233  const int ij = *ij_itr;
1234  EBounds x_init_ij;
1235  if( ij < 0 ) {
1236  const size_type i = -ij;
1237  TEUCHOS_TEST_FOR_EXCEPT( !( i <= n ) );
1238  // [P_XF_hat](:,s) = e(i)
1239  P_XF_hat_row_[k_XF_hat] = i;
1240  P_XF_hat_col_[k_XF_hat] = s;
1241  ++k_XF_hat;
1242  }
1243  else if( !(ij <= n && (x_init_ij = x_init(ij)) != FREE ) ) {
1244  const size_type j = ij;
1245  TEUCHOS_TEST_FOR_EXCEPT( !( 0 < j && j <= n + m_breve ) );
1246  // [P_plus_hat](:,s) = e(j)
1247  P_plus_hat_row_[k_plus_hat] = j;
1248  P_plus_hat_col_[k_plus_hat] = s;
1249  ++k_plus_hat;
1250  }
1251  }
1252  TEUCHOS_TEST_FOR_EXCEPT( !( k_XF_hat == q_F_hat ) );
1253  TEUCHOS_TEST_FOR_EXCEPT( !( k_plus_hat == q_plus_hat ) );
1254  }
1255  P_XF_hat_.initialize_and_sort(
1256  n,q_hat,q_F_hat,0,0,GPMSTP::BY_ROW
1257  ,q_F_hat ? &P_XF_hat_row_[0] : NULL
1258  ,q_F_hat ? &P_XF_hat_col_[0] : NULL
1259  ,test
1260  );
1261  P_plus_hat_.initialize_and_sort(
1262  n+m_breve,q_hat,q_plus_hat,0,0,GPMSTP::BY_ROW
1263  ,q_plus_hat ? &P_plus_hat_row_[0] : NULL
1264  ,q_plus_hat ? &P_plus_hat_col_[0] : NULL
1265  ,test
1266  );
1267 
1268  // Setup P_FC_hat_
1269  if( q_C_hat ) {
1270  throw std::logic_error(
1271  error_msg(__FILE__,__LINE__,"QPSchur::ActiveSet::initialize(...) : "
1272  "Error, q_C_hat != 0, now supported yet!"));
1273  // ToDo: We should implement this but it is unlikely to be needed
1274  }
1275  P_FC_hat_.initialize_and_sort(
1276  q_hat,q_hat,q_C_hat,0,0,GPMSTP::BY_ROW
1277  ,q_C_hat ? &P_FC_hat_row_[0] : NULL
1278  ,q_C_hat ? &P_FC_hat_col_[0] : NULL
1279  ,test
1280  );
1281 
1282  // Setup Q_XD_hat_
1283  Q_XD_hat_row_.resize(q_D_hat_max);
1284  Q_XD_hat_col_.resize(q_D_hat_max);
1285  if(q_D_hat) {
1286  // See ConstrainedOptPack_QPSchur.hpp for description of Q_XD_hat
1287  size_type
1288  k_XD_hat = 0; // zero based
1289  GenPermMatrixSlice::const_iterator
1290  Q_X_itr = qp.Q_X().begin(); // This is sorted by row already!
1291  P_row_t::const_iterator
1292  XF_search = P_XF_hat_row_.begin(), // These are already sorted by row!
1293  XF_search_end = XF_search + q_F_hat;
1294  for( size_type l = 1; l <= n_X; ++l, ++Q_X_itr ) {
1295  const size_type i = Q_X_itr->row_i(); // Already sorted by row
1296  // Look for i in XF
1297  for( ; XF_search != XF_search_end && *XF_search < i; ++XF_search ) ;
1298  if( XF_search == XF_search_end || (XF_search < XF_search_end && *XF_search > i) ) {
1299  // We went right past i and did not find it so
1300  // this variable has not been freed so lets add it!
1301  Q_XD_hat_row_[k_XD_hat] = i;
1302  Q_XD_hat_col_[k_XD_hat] = k_XD_hat + 1;
1303  ++k_XD_hat;
1304  }
1305  }
1306  TEUCHOS_TEST_FOR_EXCEPT( !( k_XD_hat == q_D_hat ) );
1307  }
1308  Q_XD_hat_.initialize(
1309  n,q_D_hat,q_D_hat,0,0,GPMSTP::BY_ROW // Should already be sorted by row!
1310  , q_D_hat ? &Q_XD_hat_row_[0] : NULL
1311  , q_D_hat ? &Q_XD_hat_col_[0] : NULL
1312  ,test
1313  );
1314 
1315  // Setup l_fxfx
1316  l_fxfx_.resize(q_D_hat_max);
1317  if(q_D_hat) {
1318  for( size_type k = 0; k < q_D_hat; ++k ) {
1319  l_fxfx_[k] = l_x_X_map(Q_XD_hat_row_[k]);
1320  TEUCHOS_TEST_FOR_EXCEPT( !( l_fxfx_[k] != 0 ) );
1321  }
1322  }
1323 
1324  // Set the rest of the terms in d_hat involving matrices
1325  //
1326  // d_hat += - P_XF_hat' * G * b_XX - P_plus_hat' * A_bar' * b_XX
1327  //
1328  // where: b_XX = Q_X * b_X
1329  //
1330  if( q_hat ) {
1331  SpVector b_XX;
1332  V_MtV( &b_XX, qp.Q_X(), BLAS_Cpp::no_trans, b_X );
1333  Vp_StPtMtV( &d_hat_(1,q_hat), -1.0, P_XF_hat_, BLAS_Cpp::trans
1334  , G, BLAS_Cpp::no_trans, b_XX() );
1335  Vp_StPtMtV( &d_hat_(1,q_hat), -1.0, P_plus_hat_, BLAS_Cpp::trans
1336  , constraints.A_bar(), BLAS_Cpp::trans, b_XX() );
1337  }
1338 
1339  // Setup U_hat
1340  U_hat_.initialize( &G, m ? &qp.A() : NULL, &constraints.A_bar()
1341  , &qp.Q_R(), &P_XF_hat_, &P_plus_hat_ );
1342 
1343  // Set the rest of the members
1344  test_ = test;
1345  qp_ = &qp;
1346  x_init_ = &x_init;
1347  n_ = n;
1348  n_R_ = n_R;
1349  m_ = m;
1350  m_breve_ = m_breve;
1351  q_plus_hat_ = q_plus_hat;
1352  q_F_hat_ = q_F_hat;
1353  q_C_hat_ = q_C_hat;
1354 
1355  // Resize storage for z_hat, p_z_hat, mu_D_hat, and p_mu_D_hat and set to zero by default
1356  z_hat_.resize(q_hat_max);
1357  p_z_hat_.resize(q_hat_max);
1358  mu_D_hat_.resize(n_X);
1359  p_mu_D_hat_.resize(n_X);
1360 
1361  initialized_ = true; // set to true tenatively so that we can
1362  // print this stuff.
1363 
1364  if( out && (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
1365  *out
1366  << "\nPrint definition of Active-Set before the Schur complement is formed...\n";
1367  dump_act_set_quantities( *this, *out, false );
1368  }
1369 
1370  // Initialize and factorize the schur complement
1371  if( q_hat ) {
1372  // Temporary storage for S (dense)
1373  DMatrix S_store(q_hat+1,q_hat+1);
1374  DMatrixSliceSym S_sym( S_store(2,q_hat+1,1,q_hat), BLAS_Cpp::lower );
1375  MatrixSymPosDefCholFactor S(&S_store());
1376  // S = -1.0 * U_hat' * inv(Ko) * U_hat
1377  M_StMtInvMtM( &S, -1.0, U_hat_, BLAS_Cpp::trans, qp.Ko()
1378  , MatrixSymNonsing::DUMMY_ARG );
1379  // Now add parts of V_hat
1380  if( q_F_hat ) {
1381  // S += P_XF_hat' * G * P_XF_hat
1382  Mp_StPtMtP( &S, 1.0, MatrixSymOp::DUMMY_ARG, qp_->G(), P_XF_hat_, BLAS_Cpp::no_trans );
1383  }
1384  if( q_F_hat && q_plus_hat ) {
1385  // S += P_XF_hat' * A_bar * P_plus_hat + P_plus_hat' * A_bar' * P_XF_hat
1387  qp_->constraints().A_bar()
1388  ,BLAS_Cpp::no_trans, 1.0
1389  ,P_XF_hat_, BLAS_Cpp::no_trans
1390  ,P_plus_hat_, BLAS_Cpp::no_trans
1391  ,1.0, &S );
1392  }
1393  if( q_F_hat && q_C_hat ) {
1394  // S += P_FC_hat + P_FC_hat'
1395  throw std::logic_error(
1396  error_msg(__FILE__,__LINE__,"QPSchur::ActiveSet::initialize(...) : "
1397  "Error, q_C_hat != 0, now supported yet!"));
1398  // ToDo: We should implement this but it is unlikely to be needed
1399  }
1400 
1401  if( out && (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
1402  *out
1403  << "\nIninitial Schur Complement before it is nonsingular:\n"
1404  << "\nS_hat =\nLower triangular part (ignore nonzeros above diagonal)\n"
1405  << S_store();
1406  }
1407  // Initialize and factorize the schur complement!
1408  try {
1409  schur_comp().update_interface().initialize(
1410  S_sym,q_hat_max,true,MSADU::Inertia( q_plus_hat + q_C_hat, 0, q_F_hat )
1411  ,pivot_tols() );
1412  }
1413  catch(const MSADU::WarnNearSingularUpdateException& excpt) {
1414  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1415  *out
1416  << "\nActiveSet::initialize(...) : " << excpt.what()
1417  << std::endl;
1418  }
1419  }
1420  catch(const MSADU::SingularUpdateException& excpt) {
1421  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1422  *out
1423  << "\nActiveSet::initialize(...) : " << excpt.what()
1424  << std::endl;
1425  }
1426  if(salvage_init_schur_comp) {
1427  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1428  *out
1429  << "\nsalvage_init_schur_comp == true\n"
1430  << "We will attempt to add as many rows/cols of the "
1431  << "initial Schur complement as possible ...\n";
1432  }
1433  // ToDo: We will build the schur complement one row/col at a time
1434  // skipping those updates that cause it to become singular. For each
1435  // update that causes the schur compplement to become singular we
1436  // will remove the corresponding change.
1437 
1438  throw; // For now just rethrow the exception!
1439  }
1440  else {
1441  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1442  *out
1443  << "\nsalvage_init_schur_comp == false\n"
1444  << "We will just throw this singularity exception out of here ...\n";
1445  }
1446  throw;
1447  }
1448  }
1449  if( out && (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
1450  *out
1451  << "\nSchur Complement after factorization:\n"
1452  << "\nS_hat =\n"
1453  << schur_comp().op_interface();
1454  }
1455  }
1456  else {
1457  schur_comp().update_interface().set_uninitialized();
1458  }
1459 
1460  // Success, we are initialized!
1461  initialized_ = true;
1462  return;
1463 
1464  } // end try
1465  catch(...) {
1466  initialized_ = false;
1467  throw;
1468  }
1469 }
1470 
1472 {
1473  // ToDo: Finish Me
1475 }
1476 
1478  size_type ja, EBounds bnd_ja, bool update_steps
1479  ,std::ostream *out, EOutputLevel output_level
1480  ,bool force_refactorization
1481  ,bool allow_any_cond
1482  )
1483 {
1484  using BLAS_Cpp::no_trans;
1485  using BLAS_Cpp::trans;
1486  using DenseLinAlgPack::dot;
1488  using LinAlgOpPack::V_StMtV;
1490  using LinAlgOpPack::V_InvMtV;
1491  using Teuchos::Workspace;
1493 
1494  typedef AbstractLinAlgPack::EtaVector eta_t;
1495 
1496  assert_initialized();
1497 
1498  bool wrote_output = false;
1499 
1501  &constraints = qp_->constraints();
1502 
1503  if( is_init_fixed(ja) && (*x_init_)(ja) == bnd_ja ) {
1504  //
1505  // This is a variable that was initially fixed, then freed and now
1506  // is being fixed back to its original bound.
1507  //
1508  // Here we will shrink the augmented KKT system.
1509  //
1510  const size_type q_hat = this->q_hat();
1511  const size_type sd = s_map(-int(ja));
1512  const size_type la = qp_->l_x_X_map()(ja);
1513  TEUCHOS_TEST_FOR_EXCEPT( !( sd ) );
1514  TEUCHOS_TEST_FOR_EXCEPT( !( la ) );
1515  wrote_output = remove_augmented_element(
1516  sd,force_refactorization
1517  ,MatrixSymAddDelUpdateable::EIGEN_VAL_POS
1518  ,out,output_level,allow_any_cond);
1519  // We must remove (ja,sd) from P_XF_hat
1520  P_row_t::iterator
1521  itr = std::lower_bound( P_XF_hat_row_.begin(), P_XF_hat_row_.begin()+q_F_hat_, ja );
1522  TEUCHOS_TEST_FOR_EXCEPT( !( itr != P_XF_hat_row_.end() ) );
1523  const size_type p = itr - P_XF_hat_row_.begin();
1524  std::copy( P_XF_hat_row_.begin() + p + 1, P_XF_hat_row_.begin()+q_F_hat_,
1525  P_XF_hat_row_.begin() + p );
1526  std::copy( P_XF_hat_col_.begin() + p + 1, P_XF_hat_col_.begin()+q_F_hat_,
1527  P_XF_hat_col_.begin() + p );
1528  // Deincrement all counters in permutation matrices for removed element
1529  if( q_F_hat_ > 1 )
1530  deincrement_indices( sd, &P_XF_hat_col_, q_F_hat_-1 );
1531  if( q_C_hat_ > 0 )
1532  deincrement_indices( sd, &P_FC_hat_col_, q_C_hat_ );
1533  if( q_plus_hat_ > 0 )
1534  deincrement_indices( sd, &P_plus_hat_col_, q_plus_hat_ );
1535  //
1536  // Add the multiplier for mu_D_hat(...)
1537  //
1538  const size_type q_D_hat = this->q_D_hat();
1539  // add l_fxfx(q_D_hat+1) = la
1540  l_fxfx_[q_D_hat] = la;
1541  // add mu_D_hat(q_D_hat+1) = 0.0
1542  mu_D_hat_[q_D_hat] = 0.0;
1543  // add p_mu_D_hat(q_D_hat+1) = 1.0
1544  if(update_steps)
1545  p_mu_D_hat_[q_D_hat] = 1.0;
1546  // Insert the pair (ja,q_D_hat+1) into Q_XD_hat(...) (sorted by row)
1547  insert_pair_sorted(ja,q_D_hat+1,q_D_hat+1,&Q_XD_hat_row_,&Q_XD_hat_col_);
1548  //
1549  // Update the counts
1550  //
1551  --q_F_hat_;
1552  }
1553  else {
1554  //
1555  // Expand the size of the schur complement to add the new constraint
1556  //
1557 
1558  // Compute the terms for the update
1559 
1560  value_type d_p = 0.0;
1561  const size_type q_hat = this->q_hat();
1562  Workspace<value_type> t_hat_ws(wss,q_hat);
1563  DVectorSlice t_hat(t_hat_ws.size()?&t_hat_ws[0]:NULL,t_hat_ws.size());
1564  value_type alpha_hat = 0.0;
1565  bool changed_bounds = false;
1566  size_type sd = 0; // Only used if changed_bounds == true
1567 
1568  if( ja <= n_ && !is_init_fixed(ja) ) {
1569  //
1570  // Fix an initially free variable is being fixed
1571  //
1572  // u_p = [ Q_R' * e(ja) ] <: R^(n_R+m)
1573  // [ 0 ]
1574  //
1575  const size_type
1576  la = qp_->Q_R().lookup_col_j(ja);
1577  TEUCHOS_TEST_FOR_EXCEPT( !( la ) );
1578  const eta_t u_p = eta_t(la,n_R_+m_);
1579  // r = inv(Ko)*u_p
1580  DVector r; // ToDo: Make this sparse!
1581  V_InvMtV( &r, qp_->Ko(), no_trans, u_p() );
1582  // t_hat = - U_hat' * r
1583  if(q_hat)
1584  V_StMtV( &t_hat(), -1.0, U_hat_, trans, r() );
1585  // alpha_hat = - u_p ' * r
1586  alpha_hat = - dot( u_p(), r() );
1587  // d_p = \bar{b}_{j^{(+)}}
1588  d_p = constraints.get_bnd( ja, bnd_ja );
1589 
1590  changed_bounds = false;
1591  }
1592  else if ( is_init_fixed(ja) ) {
1593  //
1594  // An intially fixed variable was freed and
1595  // is now being fixed to the other bound.
1596  //
1597  // Here we must expand the augmented KKT system for this
1598  // simple change.
1599  //
1600  // u_p = 0
1601  //
1602  // v_p = e(sd) <: R^(q_hat), where sd = s_map(-ja)
1603  //
1604  sd = s_map(-int(ja));
1605  TEUCHOS_TEST_FOR_EXCEPT( !( sd ) );
1606  const size_type la = qp_->l_x_X_map()(ja);
1607  TEUCHOS_TEST_FOR_EXCEPT( !( la ) );
1608  // t_hat = e(sd)
1609  t_hat = 0.0;
1610  t_hat(sd) = 1.0;
1611  // alpha_hat = 0.0
1612  alpha_hat = 0.0;
1613  // d_p = \bar{b}_{j^{(+)}} - b_X(la)
1614  d_p = constraints.get_bnd( ja, bnd_ja ) - qp_->b_X()(la);
1615 
1616  changed_bounds = true;
1617  }
1618  else { // ja > n
1619  //
1620  // Add an extra equality or inequality constraint.
1621  //
1622  // u_p = [ Q_R' * A_bar * e(ja) ] n_R
1623  // [ 0 ] m
1624  const eta_t e_ja = eta_t(ja,n_+m_breve_);
1625  const MatrixOp &A_bar = constraints.A_bar();
1626  DVector u_p( n_R_ + m_ ); // ToDo: make this sparse
1627  Vp_StPtMtV( &u_p(1,n_R_), 1.0, qp_->Q_R(), trans, A_bar, no_trans, e_ja(), 0.0 );
1628  if( m_ )
1629  u_p(n_R_+1,n_R_+m_) = 0.0;
1630  // r = inv(Ko) * u_p
1631  DVector r; // ToDo: Make this sparse!
1632  V_InvMtV( &r, qp_->Ko(), no_trans, u_p() );
1633  if(q_hat) {
1634  // t_hat = v_p - U_hat' * r
1635  // where: v_p = P_XF_hat' * A_bar * e(ja)
1636  V_StMtV( &t_hat(), -1.0, U_hat_, trans, r() );
1637  Vp_StPtMtV( &t_hat(), 1.0, P_XF_hat_, trans, A_bar, no_trans, e_ja() );
1638  }
1639  // alpha_hat = - u_p ' * r
1640  alpha_hat = - dot( u_p(), r() );
1641  // d_p = \bar{b}_{j^{(+)}} - b_X' * Q_X' * A_bar * e(ja)
1642  //
1643  // d_p = \bar{b}_{j^{(+)}}
1644  d_p = constraints.get_bnd( ja, bnd_ja );
1645  if(n_ > n_R_) {
1646  // d_p += - b_X' * Q_X' * A_bar * e(ja)
1647  r.resize( n_ - n_R_ ); // reuse storage
1648  Vp_StPtMtV( &r(), 1.0, qp_->Q_X(), trans, A_bar, no_trans, e_ja(), 0.0 );
1649  d_p += - dot( qp_->b_X(), r() );
1650  }
1651 
1652  changed_bounds = false;
1653  }
1654 
1655  // Update the schur complement if nonsingular. These
1656  // with throw exceptions if the matrix is singular
1657  // or has the wrong inertia
1658  if(q_hat) {
1659  try {
1660  schur_comp().update_interface().augment_update(
1661  &t_hat(), alpha_hat, force_refactorization
1662  ,MatrixSymAddDelUpdateable::EIGEN_VAL_NEG
1663  ,MSADU::PivotTolerances(
1664  pivot_tols().warning_tol
1665  ,allow_any_cond ? 0.0 : pivot_tols().singular_tol
1666  ,allow_any_cond ? 0.0 : pivot_tols().wrong_inertia_tol ) );
1667  }
1668  catch(const MSADU::WarnNearSingularUpdateException& excpt) {
1669  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1670  *out
1671  << "\nActiveSet::add_constraint(...) : " << excpt.what()
1672  << std::endl;
1673  wrote_output = true;
1674  }
1675  }
1676  }
1677  else {
1678  schur_comp().update_interface().initialize(
1679  alpha_hat, (n_-n_R_) + n_-m_ );
1680  }
1681  // Update the rest of the augmented KKT system
1682  if( changed_bounds )
1683  ++q_C_hat_;
1684  else
1685  ++q_plus_hat_;
1686  const size_type q_hat_new = q_F_hat_ + q_C_hat_ + q_plus_hat_;
1687  // Add ij_map(q_hat) = ja to ij_map(...)
1688  ij_map_[q_hat_new - 1] = ja;
1689  // Add constr_norm(q_hat) to constr_norm(...)
1690  constr_norm_[q_hat_new - 1] = 1.0; // ToDo: Compute this for real!
1691  // Add bnds(q_hat)_ = bnd_ja to bnds(...)
1692  bnds_[q_hat_new - 1] = bnd_ja;
1693  // Augment d_hat = [ d_hat; d_p ]
1694  d_hat_(q_hat_new) = d_p;
1695  // Augment z_hat with new (zeroed) multiplier value, z_hat = [ z_hat; 0 ]
1696  z_hat_(q_hat_new) = 0.0;
1697  if( update_steps ) {
1698  // Set the step for this multiplier to 1.0, p_z_hat = [ p_z_hat; 1 ]
1699  p_z_hat_(q_hat_new) = 1.0;
1700  }
1701  if( !changed_bounds ) {
1702  // Insert (ja, q_hat_new) into P_plus_hat, sorted by row
1703  insert_pair_sorted(ja,q_hat_new,q_plus_hat_,&P_plus_hat_row_,&P_plus_hat_col_);
1704  }
1705  else {
1706  TEUCHOS_TEST_FOR_EXCEPT( !( sd ) );
1707  // Insert (sd,q_hat_new) into P_FC_hat, sorted by row)
1708  insert_pair_sorted(sd,q_hat_new,q_C_hat_,&P_FC_hat_row_,&P_FC_hat_col_);
1709  }
1710  }
1711  // Update the permutation matrices and U_hat
1712  reinitialize_matrices(test_);
1713  return wrote_output;
1714 }
1715 
1717  int jd, std::ostream *out, EOutputLevel output_level
1718  ,bool force_refactorization, bool allow_any_cond
1719  )
1720 {
1721  using BLAS_Cpp::no_trans;
1722  using BLAS_Cpp::trans;
1723  using DenseLinAlgPack::dot;
1725  using LinAlgOpPack::V_StMtV;
1726  using LinAlgOpPack::V_MtV;
1727  using LinAlgOpPack::Vp_MtV;
1729  using LinAlgOpPack::V_InvMtV;
1731  using Teuchos::Workspace;
1733 
1734  typedef AbstractLinAlgPack::EtaVector eta_t;
1735 
1736  assert_initialized();
1737 
1738  bool wrote_output = false;
1739 
1740  if( jd < 0 ) {
1741  //
1742  // A variable initially fixed is being freed.
1743  // Increase the dimension of the augmented the KKT system!
1744  //
1745  size_type
1746  q_hat = this->q_hat(),
1747  q_F_hat = this->q_F_hat(),
1748  q_plus_hat = this->q_plus_hat(),
1749  q_D_hat = this->q_D_hat();
1750  // Get indexes
1751  const size_type id = -jd;
1752  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= id && id <= n_ ) );
1753  const size_type ld = qp_->l_x_X_map()(-jd);
1754  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= ld && ld <= n_ - n_R_ ) );
1755  size_type kd; // Find kd (this is unsorted)
1756  {for( kd = 1; kd <= q_D_hat; ++kd ) {
1757  if( l_fxfx_[kd-1] == ld ) break;
1758  }}
1759  TEUCHOS_TEST_FOR_EXCEPT( !( kd <= q_D_hat ) );
1760  // Get references
1761  const MatrixSymOp
1762  &G = qp_->G();
1763  const DVectorSlice
1764  g = qp_->g();
1765  const MatrixOp
1766  &A_bar = qp_->constraints().A_bar();
1767  const MatrixSymOpNonsing
1768  &Ko = qp_->Ko();
1769  const MatrixOp
1770  &U_hat = this->U_hat();
1771  const GenPermMatrixSlice
1772  &Q_R = qp_->Q_R(),
1773  &Q_X = qp_->Q_X(),
1774  &P_XF_hat = this->P_XF_hat(),
1775  &P_plus_hat = this->P_plus_hat();
1776  const DVectorSlice
1777  b_X = qp_->b_X();
1778  //
1779  // Compute the update quantities to augmented KKT system
1780  //
1781  // e_id
1782  eta_t e_id(id,n_);
1783  // u_p = [ Q_R'*G*e_id ; A'*e_id ] <: R^(n_R+m)
1784  DVector u_p(n_R_+m_);
1785  Vp_StPtMtV( &u_p(1,n_R_), 1.0, Q_R, trans, G, no_trans, e_id(), 0.0 );
1786  if( m_ )
1787  V_MtV( &u_p(n_R_+1,n_R_+m_), qp_->A(), trans, e_id() );
1788  const value_type
1789  nrm_u_p = DenseLinAlgPack::norm_inf( u_p() );
1790  // sigma = e_id'*G*e_id <: R
1791  const value_type
1792  sigma = transVtMtV( e_id(), G, no_trans, e_id() );
1793  // d_p = - g(id) - b_X'*(Q_X'*G*e_id) <: R
1794  DVector Q_X_G_e_id(Q_X.cols());
1795  Vp_StPtMtV( &Q_X_G_e_id(), 1.0, Q_X, trans, G, no_trans, e_id(), 0.0 );
1796  const value_type
1797  d_p = -g(id) - dot( b_X, Q_X_G_e_id() );
1798  // r = inv(Ko)*u_p <: R^(n_R+m)
1799  DVector r;
1800  if( nrm_u_p > 0.0 )
1801  V_InvMtV( &r, Ko, no_trans, u_p() );
1802  // t_hat = v_p - U_hat'*r
1803  // where: v_p = P_XF_hat'*G*e_id + P_plus_hat'*A_bar'*e_id <: R^(q_hat)
1804  Workspace<value_type>
1805  t_hat_ws(wss,q_hat);
1806  DVectorSlice
1807  t_hat(&t_hat_ws[0],q_hat);
1808  if(q_hat) {
1809  t_hat = 0.0;
1810  // t_hat += v_p
1811  if(q_F_hat_)
1812  Vp_StPtMtV( &t_hat(), 1.0, P_XF_hat, trans, G, no_trans, e_id() );
1813  if(q_plus_hat_)
1814  Vp_StPtMtV( &t_hat(), 1.0, P_plus_hat, trans, A_bar, trans, e_id() );
1815  // t_hat += U_hat'*r
1816  if( nrm_u_p > 0.0 )
1817  Vp_MtV( &t_hat(), U_hat, trans, r() );
1818  }
1819  // alpha_hat = sigma - u_p'*r
1820  const value_type
1821  alpha_hat = sigma - ( nrm_u_p > 0.0 ? dot(u_p(),r()) : 0.0 );
1822  //
1823  // Update the schur complement if nonsingular. These
1824  // with throw exceptions if the matrix is singular
1825  // or has the wrong inertia
1826  //
1827  if(q_hat) {
1828  try {
1829  schur_comp().update_interface().augment_update(
1830  &t_hat(), alpha_hat, force_refactorization
1831  ,MatrixSymAddDelUpdateable::EIGEN_VAL_POS );
1832  }
1833  catch(const MSADU::WarnNearSingularUpdateException& excpt) {
1834  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
1835  *out
1836  << "\nActiveSet::drop_constraint(...) : " << excpt.what()
1837  << std::endl;
1838  wrote_output = true;
1839  }
1840  }
1841  }
1842  else {
1843  schur_comp().update_interface().initialize(
1844  alpha_hat, (n_-n_R_) + n_-m_ );
1845  }
1846  //
1847  // Remove multiplier from mu_D_hat(...)
1848  //
1849  // remove l_fxfx(kd) == ld from l_fxfx(...)
1850  std::copy( l_fxfx_.begin() + kd, l_fxfx_.begin() + q_D_hat
1851  , l_fxfx_.begin() + (kd-1) );
1852  // remove mu_D_hat(kd) from mu_D_hat(...)
1853  std::copy( mu_D_hat_.begin() + kd, mu_D_hat_.begin() + q_D_hat
1854  , mu_D_hat_.begin() + (kd-1) );
1855  // remove Q_XD_hat(id,ld) from Q_XD_hat(...)
1856  P_row_t::iterator
1857  itr = std::lower_bound( Q_XD_hat_row_.begin(), Q_XD_hat_row_.begin()+q_D_hat, id );
1858  TEUCHOS_TEST_FOR_EXCEPT( !( itr != Q_XD_hat_row_.end() ) );
1859  const size_type p = itr - Q_XD_hat_row_.begin();
1860  std::copy( Q_XD_hat_row_.begin() + p + 1, Q_XD_hat_row_.begin()+q_D_hat,
1861  Q_XD_hat_row_.begin() + p );
1862  std::copy( Q_XD_hat_col_.begin() + p + 1, Q_XD_hat_col_.begin()+q_D_hat,
1863  Q_XD_hat_col_.begin() + p );
1864  if( q_D_hat > 1 )
1865  deincrement_indices( kd, &Q_XD_hat_col_, q_D_hat-1 );
1866  //
1867  // Update the counts
1868  //
1869  ++q_F_hat_;
1870  q_hat = this->q_hat();
1871  //
1872  // Add the elements for newly freed variable
1873  //
1874  // add ij_map(q_hat) == -id to ij_map(...)
1875  ij_map_[q_hat-1] = -id;
1876  // add s_map(-id) == q_hat to s_map(...)
1877  // ToDo: implement s_map(...)
1878  // add bnds(q_hat) == FREE to bnds(...)
1879  bnds_[q_hat-1] = FREE;
1880  // add d_hat(q_hat) == d_p to d_hat(...)
1881  d_hat_[q_hat-1] = d_p;
1882  // add p_X(ld) == 0 to the end of z_hat(...)
1883  z_hat_[q_hat-1] = 0.0; // This is needed so that (z_hat + beta*t_D*p_z_hat)(q_hat) == 0
1884  // Insert (id,q_hat) into P_XF_hat sorted by row
1885  insert_pair_sorted(id,q_hat,q_F_hat_,&P_XF_hat_row_,&P_XF_hat_col_);
1886  }
1887  else {
1888  //
1889  // Shrink the dimension of the augmented KKT system to remove the constraint!
1890  //
1891  const size_type q_hat = this->q_hat();
1892  const size_type sd = s_map(jd);
1893  TEUCHOS_TEST_FOR_EXCEPT( !( sd ) );
1894  wrote_output = remove_augmented_element(
1895  sd,force_refactorization
1896  ,MatrixSymAddDelUpdateable::EIGEN_VAL_NEG
1897  ,out,output_level,allow_any_cond
1898  );
1899  if( is_init_fixed(jd) ) {
1900  // This must be an intially fixed variable, currently fixed at a different bound.
1901  // We must remove this element from P_FC_hat(...)
1902  const size_type sd1 = s_map(-jd); // This is the position in the schur complement where first freed
1903  TEUCHOS_TEST_FOR_EXCEPT( !( sd1 ) );
1904  // Remove P_FC_hat(sd1,sd) from P_FC_hat(...)
1905  P_row_t::iterator
1906  itr = std::lower_bound( P_FC_hat_row_.begin(), P_FC_hat_row_.begin()+q_C_hat_, sd1 );
1907  TEUCHOS_TEST_FOR_EXCEPT( !( itr != P_FC_hat_row_.end() ) );
1908  const size_type p = itr - P_FC_hat_row_.begin();
1909  std::copy( P_FC_hat_row_.begin() + p + 1, P_FC_hat_row_.begin()+q_C_hat_,
1910  P_FC_hat_row_.begin() + p );
1911  std::copy( P_FC_hat_col_.begin() + p + 1, P_FC_hat_col_.begin()+q_C_hat_,
1912  P_FC_hat_col_.begin() + p );
1913  --q_C_hat_;
1914  }
1915  else {
1916  // We must remove P_plus_hat(jd,sd) from P_plus_hat(...)
1917  P_row_t::iterator
1918  itr = std::lower_bound( P_plus_hat_row_.begin(), P_plus_hat_row_.begin()+q_plus_hat_, jd );
1919  TEUCHOS_TEST_FOR_EXCEPT( !( itr != P_plus_hat_row_.end() ) );
1920  const size_type p = itr - P_plus_hat_row_.begin();
1921  std::copy( P_plus_hat_row_.begin() + p + 1, P_plus_hat_row_.begin()+q_plus_hat_,
1922  P_plus_hat_row_.begin() + p );
1923  std::copy( P_plus_hat_col_.begin() + p + 1, P_plus_hat_col_.begin()+q_plus_hat_,
1924  P_plus_hat_col_.begin() + p );
1925  --q_plus_hat_;
1926  }
1927  // Deincrement all counters in permutation matrices for removed element
1928  if( q_F_hat_ > 0 )
1929  deincrement_indices( sd, &P_XF_hat_col_, q_F_hat_ );
1930  if( q_C_hat_ > 0 )
1931  deincrement_indices( sd, &P_FC_hat_col_, q_C_hat_ );
1932  if( q_plus_hat_ > 0 )
1933  deincrement_indices( sd, &P_plus_hat_col_, q_plus_hat_ );
1934  }
1935  // Update the permutation matrices and U_hat
1936  reinitialize_matrices(test_);
1937  return wrote_output;
1938 }
1939 
1941  int jd, size_type ja, EBounds bnd_ja, bool update_steps
1942  ,std::ostream *out, EOutputLevel output_level
1943  )
1944 {
1945  bool wrote_output = false;
1946  if( drop_constraint( jd, out, output_level, false, true ) )
1947  wrote_output = true;
1948  if( add_constraint( ja, bnd_ja, update_steps, out, output_level, true, true ) )
1949  wrote_output = true;
1950  return wrote_output;
1951 }
1952 
1955 {
1956  assert_initialized();
1957  return *qp_;
1958 }
1959 
1962 {
1963  assert_initialized();
1964  return *qp_;
1965 }
1966 
1967 size_type QPSchur::ActiveSet::q_hat() const
1968 {
1969  assert_initialized();
1970  return q_plus_hat_ + q_F_hat_ + q_C_hat_;
1971 }
1972 
1974 {
1975  assert_initialized();
1976  return q_plus_hat_;
1977 }
1978 
1980 {
1981  assert_initialized();
1982  return q_F_hat_;
1983 }
1984 
1986 {
1987  assert_initialized();
1988  return q_C_hat_;
1989 }
1990 
1992 {
1993  assert_initialized();
1994  return (n_ - n_R_) - q_F_hat_; // n^{X} - \hat{q}^{F}
1995 }
1996 
1997 int QPSchur::ActiveSet::ij_map( size_type s ) const
1998 {
1999  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= s && s <= this->q_hat() ) );
2000  return ij_map_[s-1];
2001 }
2002 
2003 size_type QPSchur::ActiveSet::s_map( int ij ) const
2004 {
2005  ij_map_t::const_iterator
2006  begin = ij_map_.begin(),
2007  end = begin + q_hat(),
2008  itr = std::find( begin, end, ij );
2009  return ( itr != end ? (itr - begin) + 1 : 0 );
2010 }
2011 
2012 value_type QPSchur::ActiveSet::constr_norm( size_type s ) const
2013 {
2014  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= s && s <= this->q_hat() ) );
2015  return constr_norm_(s);
2016 }
2017 
2018 EBounds QPSchur::ActiveSet::bnd( size_type s ) const
2019 {
2020  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= s && s <= this->q_hat() ) );
2021  return bnds_[s-1];
2022 }
2023 
2024 size_type QPSchur::ActiveSet::l_fxfx( size_type k ) const
2025 {
2026  TEUCHOS_TEST_FOR_EXCEPT( !( 1 <= k && k <= this->q_D_hat() ) );
2027  return l_fxfx_[k-1];
2028 }
2029 
2031 {
2032  assert_initialized();
2033  return U_hat_;
2034 }
2035 
2036 const MatrixSymOpNonsing& QPSchur::ActiveSet::S_hat() const
2037 {
2038  assert_initialized();
2039  return schur_comp().op_interface();
2040 }
2041 
2042 const GenPermMatrixSlice& QPSchur::ActiveSet::P_XF_hat() const
2043 {
2044  assert_initialized();
2045  return P_XF_hat_;
2046 }
2047 
2048 const GenPermMatrixSlice& QPSchur::ActiveSet::P_FC_hat() const
2049 {
2050  assert_initialized();
2051  return P_FC_hat_;
2052 }
2053 
2054 const GenPermMatrixSlice& QPSchur::ActiveSet::P_plus_hat() const
2055 {
2056  assert_initialized();
2057  return P_plus_hat_;
2058 }
2059 
2060 const GenPermMatrixSlice& QPSchur::ActiveSet::Q_XD_hat() const
2061 {
2062  assert_initialized();
2063  return Q_XD_hat_;
2064 }
2065 
2066 const DVectorSlice QPSchur::ActiveSet::d_hat() const
2067 {
2068  assert_initialized();
2069  return d_hat_(1,q_hat());
2070 }
2071 
2073 {
2074  assert_initialized();
2075  return z_hat_(1,q_hat());
2076 }
2077 
2078 const DVectorSlice QPSchur::ActiveSet::z_hat() const
2079 {
2080  assert_initialized();
2081  return z_hat_(1,q_hat());
2082 }
2083 
2085 {
2086  assert_initialized();
2087  return p_z_hat_(1,q_hat());
2088 }
2089 
2090 const DVectorSlice QPSchur::ActiveSet::p_z_hat() const
2091 {
2092  assert_initialized();
2093  return p_z_hat_(1,q_hat());
2094 }
2095 
2097 {
2098  assert_initialized();
2099  return mu_D_hat_(1,q_D_hat());
2100 }
2101 
2102 const DVectorSlice QPSchur::ActiveSet::mu_D_hat() const
2103 {
2104  assert_initialized();
2105  return mu_D_hat_(1,q_D_hat());
2106 }
2107 
2109 {
2110  assert_initialized();
2111  return p_mu_D_hat_(1,q_D_hat());
2112 }
2113 
2114 const DVectorSlice QPSchur::ActiveSet::p_mu_D_hat() const
2115 {
2116  assert_initialized();
2117  return p_mu_D_hat_(1,q_D_hat());
2118 }
2119 
2120 bool QPSchur::ActiveSet::is_init_fixed( size_type j ) const
2121 {
2122  assert_initialized();
2123  return j <= n_ && (*x_init_)(j) != FREE;
2124 }
2125 
2127 {
2128  return n_ - m_ == (n_ - n_R_) - q_F_hat_ + q_C_hat_ + q_plus_hat_;
2129 }
2130 
2131 // private member functions for QPSchur::ActiveSet
2132 
2133 void QPSchur::ActiveSet::assert_initialized() const
2134 {
2135  if( !initialized_ )
2136  throw std::logic_error(
2137  error_msg(__FILE__,__LINE__,"QPSchur::ActiveSet::assert_initialized() : Error, "
2138  "The active set has not been initialized yet!") );
2139 }
2140 
2141 void QPSchur::ActiveSet::assert_s( size_type s) const
2142 {
2143  TEUCHOS_TEST_FOR_EXCEPT( !( s <= q_hat() ) ); // ToDo: Throw an exception
2144 }
2145 
2146 void QPSchur::ActiveSet::reinitialize_matrices(bool test)
2147 {
2148  namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
2149 
2150  const size_type q_hat = this->q_hat();
2151  const size_type q_D_hat = this->q_D_hat();
2152 
2153  P_XF_hat_.initialize(
2154  n_,q_hat,q_F_hat_,0,0,GPMSTP::BY_ROW
2155  ,q_F_hat_ ? &P_XF_hat_row_[0] : NULL
2156  ,q_F_hat_ ? &P_XF_hat_col_[0] : NULL
2157  ,test
2158  );
2159  P_FC_hat_.initialize(
2160  q_hat,q_hat,q_C_hat_,0,0,GPMSTP::BY_ROW
2161  ,q_C_hat_ ? &P_FC_hat_row_[0] : NULL
2162  ,q_C_hat_ ? &P_FC_hat_col_[0] : NULL
2163  ,test
2164  );
2165  P_plus_hat_.initialize(
2166  n_+m_breve_,q_hat,q_plus_hat_,0,0,GPMSTP::BY_ROW
2167  ,q_plus_hat_ ? &P_plus_hat_row_[0] : NULL
2168  ,q_plus_hat_ ? &P_plus_hat_col_[0] : NULL
2169  ,test
2170  );
2171  Q_XD_hat_.initialize(
2172  n_,q_D_hat,q_D_hat,0,0,GPMSTP::BY_ROW
2173  ,q_D_hat ? &Q_XD_hat_row_[0] : NULL
2174  ,q_D_hat ? &Q_XD_hat_col_[0] : NULL
2175  ,test
2176  );
2177  U_hat_.initialize(
2178  &qp_->G(), m_ ? &qp_->A() : NULL, &qp_->constraints().A_bar()
2179  ,&qp_->Q_R(), &P_XF_hat_, &P_plus_hat_);
2180 }
2181 
2182 bool QPSchur::ActiveSet::remove_augmented_element(
2183  size_type sd, bool force_refactorization
2184  ,MatrixSymAddDelUpdateable::EEigenValType eigen_val_drop
2185  ,std::ostream *out, EOutputLevel output_level
2186  ,bool allow_any_cond
2187  )
2188 {
2189  bool wrote_output = false;
2190  const size_type q_hat = this->q_hat();
2191  // Delete the sd row and column for S_hat
2192  try {
2193  schur_comp().update_interface().delete_update(
2194  sd,force_refactorization,eigen_val_drop
2195  ,MSADU::PivotTolerances(
2196  pivot_tols().warning_tol
2197  ,allow_any_cond ? 0.0 : pivot_tols().singular_tol
2198  ,allow_any_cond ? 0.0 : pivot_tols().wrong_inertia_tol ));
2199  }
2200  catch(const MSADU::WarnNearSingularUpdateException& excpt) {
2201  if( out && (int)output_level >= QPSchur::OUTPUT_BASIC_INFO ) {
2202  *out
2203  << "\nActiveSet::drop_constraint(...) : " << excpt.what()
2204  << std::endl;
2205  wrote_output = true;
2206  }
2207  }
2208  // Remove the ij_map(s) = jd element from ij_map(...)
2209  std::copy( ij_map_.begin() + sd, ij_map_.begin() + q_hat
2210  , ij_map_.begin() + (sd-1) );
2211  // Remove the constr_norm(s) elment from constr_norm(...)
2212  std::copy( constr_norm_.begin() + sd, constr_norm_.begin() + q_hat
2213  , constr_norm_.begin() + (sd-1) );
2214  // Remove the bnds(s) element from bnds(...)
2215  std::copy( bnds_.begin() + sd, bnds_.begin() + q_hat
2216  , bnds_.begin() + (sd-1) );
2217  // Remove the d_hat(s) element from d_hat(...)
2218  std::copy( d_hat_.begin() + sd, d_hat_.begin() + q_hat
2219  , d_hat_.begin() + (sd-1) );
2220  // Remove the z_hat(s) element from z_hat(...)
2221  std::copy( z_hat_.begin() + sd, z_hat_.begin() + q_hat
2222  , z_hat_.begin() + (sd-1) );
2223  // Remove the p_z_hat(s) element from p_z_hat(...)
2224  std::copy( p_z_hat_.begin() + sd, p_z_hat_.begin() + q_hat
2225  , p_z_hat_.begin() + (sd-1) );
2226  return wrote_output;
2227 }
2228 
2229 // public member functions for QPSchur
2230 
2231 value_type QPSchur::DEGENERATE_MULT = std::numeric_limits<value_type>::min();
2232 
2233 void QPSchur::pivot_tols( MSADU::PivotTolerances pivot_tols )
2234 {
2235  act_set_.pivot_tols(pivot_tols);
2236 }
2237 
2238 QPSchur::MSADU::PivotTolerances QPSchur::pivot_tols() const
2239 {
2240  return act_set_.pivot_tols();
2241 }
2242 
2244  const schur_comp_ptr_t& schur_comp
2245  ,size_type max_iter
2246  ,value_type max_real_runtime
2247  ,value_type feas_tol
2248  ,value_type loose_feas_tol
2249  ,value_type dual_infeas_tol
2250  ,value_type huge_primal_step
2251  ,value_type huge_dual_step
2252  ,value_type warning_tol
2253  ,value_type error_tol
2254  ,size_type iter_refine_min_iter
2255  ,size_type iter_refine_max_iter
2256  ,value_type iter_refine_opt_tol
2257  ,value_type iter_refine_feas_tol
2258  ,bool iter_refine_at_solution
2259  ,bool salvage_init_schur_comp
2260  ,MSADU::PivotTolerances pivot_tols
2261  )
2262  :schur_comp_(schur_comp)
2263  ,max_iter_(max_iter)
2264  ,max_real_runtime_(max_real_runtime)
2265  ,feas_tol_(feas_tol)
2266  ,loose_feas_tol_(loose_feas_tol)
2267  ,dual_infeas_tol_(dual_infeas_tol)
2268  ,huge_primal_step_(huge_primal_step)
2269  ,huge_dual_step_(huge_dual_step)
2270  ,warning_tol_(warning_tol)
2271  ,error_tol_(error_tol)
2272  ,iter_refine_min_iter_(iter_refine_min_iter)
2273  ,iter_refine_max_iter_(iter_refine_max_iter)
2274  ,iter_refine_opt_tol_(iter_refine_opt_tol)
2275  ,iter_refine_feas_tol_(iter_refine_feas_tol)
2276  ,iter_refine_at_solution_(iter_refine_at_solution)
2277  ,salvage_init_schur_comp_(salvage_init_schur_comp)
2278  ,act_set_(schur_comp,pivot_tols)
2279 {}
2280 
2282  QP& qp
2283  ,size_type num_act_change, const int ij_act_change[], const EBounds bnds[]
2284  ,std::ostream *out, EOutputLevel output_level, ERunTests test_what
2285  ,DVectorSlice* x, SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
2286  ,size_type* iter, size_type* num_adds, size_type* num_drops
2287  )
2288 {
2289  using std::setw;
2290  using std::endl;
2291  using std::right;
2292  using DenseLinAlgPack::norm_inf;
2293  using AbstractLinAlgPack::norm_inf;
2294  using LinAlgOpPack::V_InvMtV;
2295  using Teuchos::Workspace;
2298 
2299  const value_type inf = std::numeric_limits<value_type>::max();
2300 
2301  if( !out )
2302  output_level = NO_OUTPUT;
2303 
2304  const int dbl_min_w = 20;
2305  const int dbl_w = (out ? my_max(dbl_min_w,int(out->precision()+8)) : 20 );
2306 
2307  // Set the schur complement
2308  act_set_.set_schur_comp( schur_comp_ );
2309 
2310  ESolveReturn
2311  solve_return = SUBOPTIMAL_POINT;
2312 
2313  // Print QPSchur output header
2314  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2315  *out
2316  << "\n*** Entering QPSchur::solve_qp(...)\n";
2317  }
2318 
2319  // Start the timer!
2320  stopwatch timer;
2321  timer.reset();
2322  timer.start();
2323 
2324  // Print the definition of the QP to be solved.
2325  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2326  *out
2327  << "\n*** Dump the definition of the QP to be solved ***\n";
2328  qp.dump_qp(*out);
2329  }
2330 
2331  // Print warm start info.
2332  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2333  *out
2334  << "\n*** Warm start info\n"
2335  << "\nNumber of variables = " << qp.n()
2336  << "\nNumber of initially fixed variables (not in Ko) = " << (qp.n() - qp.n_R())
2337  << "\nNumber of changes to the initial KKT system (num_act_change) = " << num_act_change << endl;
2338  const size_type
2339  n = qp.n();
2340  size_type
2341  num_var_fixed = 0, num_var_freed = 0, num_gen_equ = 0, num_gen_inequ = 0;
2342  for( size_type s = 1; s <= num_act_change; ++s ) {
2343  const int ij = ij_act_change[s-1];
2344  const EBounds bnd = bnds[s-1];
2345  if( ij < 0 )
2346  ++num_var_freed;
2347  else if( ij < n )
2348  ++num_var_fixed;
2349  else if( bnd == EQUALITY )
2350  ++num_gen_equ;
2351  else
2352  ++num_gen_inequ;
2353  }
2354  *out
2355  << "\n Number of initially fixed variables freed from a bound = " << num_var_freed
2356  << "\n Number of initially free variables fixed to a bound = " << num_var_fixed
2357  << "\n Number of general equality constraints added = " << num_gen_equ
2358  << "\n Number of general inequality constraints added = " << num_gen_inequ << endl;
2359  }
2360 
2361  if( num_act_change > 0 && (int)output_level >= (int)OUTPUT_ACT_SET ) {
2362  *out
2363  << endl
2364  << right << setw(5) << "s"
2365  << right << setw(20) << "ij_act_change"
2366  << right << setw(10) << "bnds" << endl
2367  << right << setw(5) << "---"
2368  << right << setw(20) << "-------------"
2369  << right << setw(10) << "--------" << endl;
2370  for( size_type s = 1; s <= num_act_change; ++s )
2371  *out
2372  << right << setw(5) << s
2373  << right << setw(20) << ij_act_change[s-1]
2374  << right << setw(10) << bnd_str(bnds[s-1]) << endl;
2375  }
2376 
2377  // Initialize the active set.
2378  try {
2379  act_set_.initialize( qp, num_act_change, ij_act_change, bnds
2380  , test_what == RUN_TESTS, salvage_init_schur_comp(), out, output_level );
2381  // If this throws a WrongInteriaUpdateExecption it will be
2382  // thrown clean out of here!
2383  }
2384  catch( const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
2385  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2386  *out
2387  << "\n*** Error in initializing schur complement\n"
2388  << excpt.what() << std::endl
2389  << "\nSetting num_act_change = 0 and proceeding with a cold start...\n";
2390  }
2391  act_set_.initialize( qp, num_act_change = 0, ij_act_change, bnds
2392  , test_what == RUN_TESTS, salvage_init_schur_comp(), out, output_level );
2393  }
2394 
2395  // Compute vo = inv(Ko) * fo
2396  Workspace<value_type> vo_ws(wss,qp.n_R()+qp.m());
2397  DVectorSlice vo(&vo_ws[0],vo_ws.size());
2398  V_InvMtV( &vo, qp.Ko(), BLAS_Cpp::no_trans, qp.fo() );
2399 
2400  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2401  *out
2402  << "\nSolution to the initial KKT system, vo = inv(Ko)*fo:\n\n||vo||inf = " << norm_inf(vo) << std::endl;
2403  }
2404  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2405  *out
2406  << "\nvo =\n" << vo;
2407  }
2408 
2409  // ////////////////////////////////////////////////
2410  // Remove constraints until we are dual feasible.
2411  //
2412  // Here we are essentially performing a primal algorithm where we are only
2413  // removing constraints. If the dual variables are not dual feasible then
2414  // we will remove the one with the largest scaled dual feasibility
2415  // violation then compute the dual variables again. Eventually we
2416  // will come to a point where we have a dual feasible point. If
2417  // we have to, we will remove all of the inequality constraints and then
2418  // this will by default be a dual feasible point (i.e. we picked all the
2419  // wrong inequality constraints).
2420  //
2421  // The difficulty here is in dealing with near degenerate constraints.
2422  // If a constraint is near degenerate then we would like to not drop
2423  // it since we may have to add it again later.
2424  //
2425  *iter = 0;
2426  *num_adds = 0;
2427  *num_drops = 0;
2428  // Print header for removing constraints
2429  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2430  *out
2431  << "\n***"
2432  << "\n*** Removing constriants until we are dual feasible"
2433  << "\n***\n"
2434  << "\n*** Start by removing constraints within the Schur complement first\n";
2435  }
2436  // Print summary header for max_viol and jd.
2437  if( (int)OUTPUT_ITER_SUMMARY <= (int)output_level
2438  && (int)output_level < (int)OUTPUT_ITER_QUANTITIES
2439  && num_act_change > 0 )
2440  {
2441  *out
2442  << "\nIf max_viol > 0 and jd != 0 then constraint jd will be dropped from the active set\n\n"
2443  << right << setw(dbl_w) << "max_viol"
2444  << right << setw(5) << "sd"
2445  << right << setw(5) << "jd" << endl
2446  << right << setw(dbl_w) << "--------------"
2447  << right << setw(5) << "----"
2448  << right << setw(5) << "----" << endl;
2449  }
2450  for( int k = num_act_change; k > 0; --k, ++(*iter) ) {
2451  // Check runtime
2452  if( timeout_return(&timer,out,output_level) )
2453  return MAX_RUNTIME_EXEEDED_FAIL;
2454  // Compute z_hat (z_hat = inv(S_hat)*(d_hat - U_hat'*vo))
2455  DVectorSlice z_hat = act_set_.z_hat();
2456  calc_z( act_set_.S_hat(), act_set_.d_hat(), act_set_.U_hat(), &vo
2457  , &z_hat );
2458  // Determine if we are dual feasible.
2459  value_type max_viol = 0.0; // max scaled violation of dual feasability.
2460  size_type jd = 0; // indice of constraint with max scaled violation.
2461  DVectorSlice::iterator
2462  z_itr = z_hat.begin();
2463  const size_type q_hat = act_set_.q_hat(); // Size of schur complement system.
2464  // Print header for s, z_hat(s), bnd(s), viol, max_viol and jd
2465  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2466  *out
2467  << "\nLooking for a constraint with the maximum dual infeasiblity to drop...\n\n"
2468  << right << setw(5) << "s"
2469  << right << setw(dbl_w) << "z_hat"
2470  << right << setw(20) << "bnd"
2471  << right << setw(dbl_w) << "viol"
2472  << right << setw(dbl_w) << "max_viol"
2473  << right << setw(5) << "jd" << endl
2474  << right << setw(5) << "----"
2475  << right << setw(dbl_w) << "--------------"
2476  << right << setw(20) << "--------------"
2477  << right << setw(dbl_w) << "--------------"
2478  << right << setw(dbl_w) << "--------------"
2479  << right << setw(5) << "----" << endl;
2480  }
2481  for( int s = 1; s <= q_hat; ++s, ++z_itr ) {
2482  int j = act_set_.ij_map(s);
2483  if( j > 0 ) {
2484  // This is for an active constraint not initially in Ko so z_hat(s) = mu(j)
2485  EBounds bnd = act_set_.bnd(s);
2486  value_type viol; // Get the amount of violation and fix near degeneracies
2487  const int dual_feas_status
2488  = correct_dual_infeas(
2489  j,bnd,0.0,act_set_.constr_norm(s),dual_infeas_tol(),DEGENERATE_MULT
2490  ,out,output_level,false
2491  ,"z_hat(s)", &(*z_itr), &viol );
2492  if( dual_feas_status < 0 && viol < max_viol ) {
2493  max_viol = viol;
2494  jd = j;
2495  }
2496  // Print row for s, z_hat(s), bnd(s), viol, max_viol and jd
2497  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2498  *out
2499  << right << setw(5) << s
2500  << right << setw(dbl_w) << *z_itr
2501  << right << setw(20) << bnd_str(bnd)
2502  << right << setw(dbl_w) << viol
2503  << right << setw(dbl_w) << max_viol
2504  << right << setw(5) << jd << endl;
2505  }
2506  }
2507  }
2508  // Print row of max_viol and jd
2509  if( (int)OUTPUT_ITER_SUMMARY <= (int)output_level
2510  && (int)output_level < (int)OUTPUT_ITER_QUANTITIES )
2511  {
2512  *out
2513  << right << setw(dbl_w) << max_viol
2514  << right << setw(5) << act_set_.s_map(jd)
2515  << right << setw(5) << jd << endl;
2516  }
2517  if( jd == 0 ) break; // We have a dual feasible point w.r.t. these constraints
2518  // Remove the active constraint with the largest scaled violation.
2519  act_set_.drop_constraint( jd, out, output_level, true, true );
2520  ++(*iter);
2521  ++(*num_drops);
2522  // Print U_hat, S_hat and d_hat.
2523  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2524  *out
2525  << "\nPrinting active set after dropping constraint jd = " << jd << " ...\n";
2526  dump_act_set_quantities( act_set_, *out );
2527  }
2528  }
2529 
2530  // Print how many constraints where removed from the schur complement
2531  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2532  *out
2533  << "\nThere where " << (*num_drops)
2534  << " constraints dropped from the schur complement from the initial guess of the active set.\n";
2535  }
2536 
2537  // Compute v
2538  Workspace<value_type> v_ws(wss,qp.n_R()+qp.m());
2539  DVectorSlice v(&v_ws[0],v_ws.size());
2540  if( act_set_.q_hat() > 0 ) {
2541  calc_v( qp.Ko(), &qp.fo(), act_set_.U_hat(), act_set_.z_hat(), &v );
2542  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2543  *out
2544  << "\nSolution to the system; v = inv(Ko)*(fo - U_hat*z_hat):\n"
2545  << "\n||v||inf = " << norm_inf(v) << std::endl;
2546  }
2547  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2548  *out
2549  << "\nv =\n" << v;
2550  }
2551  }
2552  else {
2553  v = vo;
2554  }
2555 
2556  // Set x
2557  set_x( act_set_, v, x );
2558  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2559  *out
2560  << "\nCurrent guess for unknowns x:\n\n||x||inf = " << norm_inf(*x) << std::endl;
2561  }
2562  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2563  *out
2564  << "\nx =\n" << *x;
2565  }
2566 
2567  //
2568  // Determine if any initially fixed variables need to be freed by checking mu_D_hat.
2569  //
2570  if( act_set_.q_D_hat() ) {
2571  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2572  *out << "\n*** Second, free initially fixed variables not in Ko\n";
2573  }
2574  const QPSchurPack::QP::i_x_X_map_t& i_x_X_map = act_set_.qp().i_x_X_map();
2575  const QPSchurPack::QP::x_init_t& x_init = act_set_.qp().x_init();
2576  // Print summary header for max_viol and jd.
2577  if( (int)OUTPUT_ITER_SUMMARY <= (int)output_level
2578  && (int)output_level < (int)OUTPUT_ITER_QUANTITIES )
2579  {
2580  *out
2581  << "\nIf max_viol > 0 and id != 0 then the variable x(id) will be freed from its initial bound\n\n"
2582  << right << setw(dbl_w) << "max_viol"
2583  << right << setw(5) << "kd"
2584  << right << setw(5) << "id" << endl
2585  << right << setw(dbl_w) << "--------------"
2586  << right << setw(5) << "----"
2587  << right << setw(5) << "----" << endl;
2588  }
2589  size_type q_D_hat = act_set_.q_D_hat(); // This will be deincremented
2590  while( q_D_hat > 0 ) {
2591  // Check runtime
2592  if( timeout_return(&timer,out,output_level) )
2593  return MAX_RUNTIME_EXEEDED_FAIL;
2594  // mu_D_hat = ???
2595  DVectorSlice mu_D_hat = act_set_.mu_D_hat();
2596  calc_mu_D( act_set_, *x, v, &mu_D_hat );
2597  // Determine if we are dual feasible.
2598  value_type max_viol = 0.0; // max scaled violation of dual feasability.
2599  int id = 0; // indice of variable with max scaled violation.
2600  size_type kd = 0;
2601  DVectorSlice::iterator
2602  mu_D_itr = mu_D_hat.begin();
2603  // Print header for k, mu_D_hat(k), bnd, viol, max_viol and id
2604  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2605  *out
2606  << "\nLooking for a variable bound with the max dual infeasibility to drop...\n\n"
2607  << right << setw(5) << "k"
2608  << right << setw(dbl_w) << "mu_D_hat"
2609  << right << setw(20) << "bnd"
2610  << right << setw(dbl_w) << "viol"
2611  << right << setw(dbl_w) << "max_viol"
2612  << right << setw(5) << "id" << endl
2613  << right << setw(5) << "----"
2614  << right << setw(dbl_w) << "--------------"
2615  << right << setw(20) << "--------------"
2616  << right << setw(dbl_w) << "--------------"
2617  << right << setw(dbl_w) << "--------------"
2618  << right << setw(5) << "----" << endl;
2619  }
2620  for( int k = 1; k <= q_D_hat; ++k, ++mu_D_itr ) {
2621  int
2622  i = i_x_X_map(act_set_.l_fxfx(k));
2623  EBounds
2624  bnd = x_init(i);
2625  value_type viol; // Get the amount of violation and fix near degeneracies
2626  const int dual_feas_status
2627  = correct_dual_infeas(
2628  i,bnd,0.0,1.0,dual_infeas_tol(),DEGENERATE_MULT
2629  ,out,output_level,false
2630  ,"mu_D_hat(k)", &(*mu_D_itr), &viol );
2631  if( dual_feas_status < 0 && viol < max_viol ) {
2632  max_viol = viol;
2633  kd = k;
2634  id = i;
2635  }
2636  // Print row for k, mu_D_hat(k), bnd, viol, max_viol and jd
2637  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2638  *out
2639  << right << setw(5) << k
2640  << right << setw(dbl_w) << *mu_D_itr
2641  << right << setw(20) << bnd_str(bnd)
2642  << right << setw(dbl_w) << viol
2643  << right << setw(dbl_w) << max_viol
2644  << right << setw(5) << id << endl;
2645  }
2646  }
2647  // Print row of max_viol and id
2648  if( (int)OUTPUT_ITER_SUMMARY <= (int)output_level
2649  && (int)output_level < (int)OUTPUT_ITER_QUANTITIES )
2650  {
2651  *out
2652  << right << setw(dbl_w) << max_viol
2653  << right << setw(5) << kd
2654  << right << setw(5) << id << endl;
2655  }
2656  if( id == 0 ) break; // We have a dual feasible point w.r.t. these variable bounds
2657  // Remove the active constraint with the largest scaled violation.
2658  act_set_.drop_constraint( -id, out, output_level, true, true );
2659  ++(*iter);
2660  ++(*num_adds);
2661  --q_D_hat;
2662  // Print U_hat, S_hat and d_hat.
2663  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2664  *out
2665  << "\nPrinting active set after freeing initially fixed variable id = " << id << " ...\n";
2666  dump_act_set_quantities( act_set_, *out );
2667  }
2668  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2669  *out
2670  << "\nSolution to the new KKT system; z_hat = inv(S_hat)*(d_hat - U_hat'*vo), v = inv(Ko)*(fo - U_hat*z_hat):\n";
2671  }
2672  // Compute z_hat (z_hat = inv(S_hat)*(d_hat - U_hat'*vo))
2673  calc_z( act_set_.S_hat(), act_set_.d_hat(), act_set_.U_hat(), &vo
2674  , &act_set_.z_hat() );
2675  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2676  *out
2677  << "\n||z_hat||inf = " << norm_inf(act_set_.z_hat()) << std::endl;
2678  }
2679  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2680  *out
2681  << "\nz_hat =\n" << act_set_.z_hat();
2682  }
2683  // Compute v
2684  calc_v( qp.Ko(), &qp.fo(), act_set_.U_hat(), act_set_.z_hat(), &v );
2685  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2686  *out
2687  << "\n||v||inf = " << norm_inf(v) << std::endl;
2688  }
2689  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2690  *out
2691  << "\nv =\n" << v;
2692  }
2693  // Set x
2694  set_x( act_set_, v, x );
2695  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2696  *out
2697  << "\nCurrent guess for unknowns x:\n\n||x||inf = " << norm_inf(*x) << std::endl;
2698  }
2699  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2700  *out
2701  << "\nx =\n" << *x;
2702  }
2703  }
2704  }
2705 
2706  // Print how many initially fixed variables where freed
2707  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2708  *out
2709  << "\nThere where " << (*num_adds)
2710  << " initially fixed variables not in Ko that were freed and added to the schur complement.\n";
2711  }
2712 
2713  // Run the primal dual algorithm
2714  size_type iter_refine_num_resid = 0, iter_refine_num_solves = 0;
2715  solve_return = qp_algo(
2716  PICK_VIOLATED_CONSTRAINT
2717  ,out, output_level, test_what
2718  ,vo, &act_set_, &v
2719  ,x, iter, num_adds, num_drops
2720  ,&iter_refine_num_resid, &iter_refine_num_solves
2721  ,&timer
2722  );
2723 
2724  if( solve_return != OPTIMAL_SOLUTION )
2725  set_x( act_set_, v, x );
2726 
2727  // Correct the sign of near degenerate multipliers in case it has not been done yet!
2728  if( solve_return != SUBOPTIMAL_POINT && act_set_.q_hat() ) {
2729  const size_type q_hat = act_set_.q_hat();
2730  DVectorSlice z_hat = act_set_.z_hat();
2731  for( size_type s = 1; s <= q_hat; ++s ) {
2732  const int j = act_set_.ij_map(s);
2733  value_type viol = 0.0;
2734  const EBounds bnd = act_set_.bnd(s);
2735  if(bnd == FREE)
2736  continue;
2737  const int dual_feas_status
2738  = correct_dual_infeas(
2739  j,bnd,0.0,1.0,dual_infeas_tol(),DEGENERATE_MULT
2740  ,out,output_level,true,"z_hat(s)",&z_hat(s),&viol );
2741  if( dual_feas_status < 0 ) {
2742  solve_return = SUBOPTIMAL_POINT;
2743  break;
2744  }
2745  }
2746  }
2747  if( solve_return != SUBOPTIMAL_POINT && act_set_.q_D_hat() ) {
2748  const GenPermMatrixSlice& Q_XD_hat = act_set_.Q_XD_hat();
2749  DVectorSlice mu_D_hat = act_set_.mu_D_hat();
2750  const QPSchurPack::QP::x_init_t& x_init = qp.x_init();
2751  for( GenPermMatrixSlice::const_iterator itr = Q_XD_hat.begin(); itr != Q_XD_hat.end(); ++itr ) {
2752  const int i = itr->row_i();
2753  value_type viol = 0.0;
2754  const EBounds bnd = x_init(i);
2755  TEUCHOS_TEST_FOR_EXCEPT( !( bnd != FREE ) );
2756  const int dual_feas_status
2757  = correct_dual_infeas(
2758  i,bnd,0.0,1.0,dual_infeas_tol(),DEGENERATE_MULT
2759  ,out,output_level,true,"mu_D_hat(k)",&mu_D_hat(itr->col_j()),&viol );
2760  if( dual_feas_status < 0 ) {
2761  solve_return = SUBOPTIMAL_POINT;
2762  break;
2763  }
2764  }
2765  }
2766 
2767  set_multipliers( act_set_, v, mu, lambda, lambda_breve );
2768 
2769  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2770  switch(solve_return) {
2771  case OPTIMAL_SOLUTION:
2772  *out
2773  << "\n*** Solution found!\n";
2774  break;
2775  case MAX_ITER_EXCEEDED:
2776  *out
2777  << "\n*** Maximum iterations exceeded!\n";
2778  break;
2779  case MAX_RUNTIME_EXEEDED_DUAL_FEAS:
2780  *out
2781  << "\n*** Maximum runtime exceeded!\n";
2782  break;
2783  case MAX_ALLOWED_STORAGE_EXCEEDED:
2784  *out
2785  << "\n*** The maxinum size of the schur complement has been exceeded!\n";
2786  break;
2787  case INFEASIBLE_CONSTRAINTS:
2788  *out
2789  << "\n*** The constraints are infeasible!\n";
2790  break;
2791  case NONCONVEX_QP:
2792  *out
2793  << "\n*** The QP appears to be nonconvex but will return current point anyway!\n";
2794  break;
2795  case DUAL_INFEASIBILITY:
2796  *out
2797  << "\n*** The dual variables are infeasible (numerical roundoff?)!\n";
2798  break;
2799  case SUBOPTIMAL_POINT:
2800  *out
2801  << "\n*** The current point is suboptimal but we will return it anyway!\n";
2802  break;
2803  default:
2805  }
2806  *out << "\nNumber of QP iteratons = " << *iter;
2807  *out << "\nNumber of iterative refinement residual calculations = " << iter_refine_num_resid;
2808  *out << "\nNumber of iterative refinement solves = " << iter_refine_num_solves << endl;
2809  *out << "\n||x||inf = " << norm_inf(*x);
2810  *out << "\nmu.nz() = " << mu->nz();
2811  *out << "\nmax(|mu(i)|) = " << norm_inf((*mu)());
2812  *out << "\nmin(|mu(i)|) = " << min_abs((*mu)());
2813  if(lambda)
2814  *out
2815  << "\nmax(|lambda(i)|) = " << norm_inf(*lambda)
2816  << "\nmin(|lambda(i)|) = " << min_abs(*lambda);
2817  if(lambda_breve)
2818  *out
2819  << "\nlambda_breve.nz() = " << lambda_breve->nz()
2820  << "\nmax(|lambda_breve(i)|) = " << norm_inf((*lambda_breve)())
2821  << "\nmin(|lambda_breve(i)|) = " << min_abs((*lambda_breve)());
2822  *out << std::endl;
2823  }
2824  // Print Solution x, lambda and mu
2825  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2826  *out << "\nx =\n" << (*x)();
2827  *out << "\nmu =\n" << (*mu)();
2828  if(lambda)
2829  *out << "\nlambda =\n" << (*lambda)();
2830  if(lambda_breve)
2831  *out << "\nlambda_breve =\n" << (*lambda_breve)();
2832  }
2833  // Print 'goodby' header.
2834  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2835  *out
2836  << "\n*** Leaving QPSchur::solve_qp(...)\n";
2837  }
2838 
2839  return solve_return;
2840 }
2841 
2843 {
2844  return act_set_;
2845 }
2846 
2847 // protected member functions for QPSchur
2848 
2850  EPDSteps next_step
2851  , std::ostream *out, EOutputLevel output_level, ERunTests test_what
2852  , const DVectorSlice& vo, ActiveSet* act_set, DVectorSlice* v
2853  , DVectorSlice* x, size_type* iter, size_type* num_adds, size_type* num_drops
2854  , size_type* iter_refine_num_resid, size_type* iter_refine_num_solves
2855  , StopWatchPack::stopwatch* timer
2856  )
2857 {
2858  using std::setw;
2859  using std::endl;
2860  using std::right;
2861  using BLAS_Cpp::no_trans;
2862  using BLAS_Cpp::trans;
2863  using DenseLinAlgPack::dot;
2864  using DenseLinAlgPack::norm_inf;
2865  using DenseLinAlgPack::Vt_S;
2866  using DenseLinAlgPack::V_mV;
2868  using DenseLinAlgPack::V_VmV;
2869  using LinAlgOpPack::Vp_V;
2870  using LinAlgOpPack::V_StV;
2871  using LinAlgOpPack::V_StMtV;
2873  using AbstractLinAlgPack::norm_inf;
2874  using AbstractLinAlgPack::V_MtV;
2876  using LinAlgOpPack::V_MtV;
2877  using LinAlgOpPack::V_InvMtV;
2878  using LinAlgOpPack::Vp_StMtV;
2880  using Teuchos::Workspace;
2882 
2883  // Print header for "Starting Primal-Dual Iterations"
2884  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
2885  *out
2886  << "\n*** Starting Primal-Dual Iterations ***\n";
2887  }
2888 
2889  const int dbl_min_w = 20;
2890  const int dbl_w = (out ? my_max(dbl_min_w,int(out->precision()+8)): 20 );
2891 
2892  try {
2893 
2895  &qp = act_set->qp();
2896  const size_type
2897  n = qp.n(),
2898  n_R = qp.n_R(),
2899  m = qp.m(),
2900  m_breve = qp.constraints().m_breve();
2901 
2902  Workspace<value_type>
2903  v_plus_ws(wss,v->dim()),
2904  z_hat_plus_ws(wss,(n-m)+(n-n_R)),
2905  p_v_ws(wss,v->dim());
2906  DVectorSlice
2907  v_plus(&v_plus_ws[0],v_plus_ws.size()),
2908  z_hat_plus,
2909  p_v(&p_v_ws[0],p_v_ws.size());
2910 
2911  const value_type
2912  inf = std::numeric_limits<value_type>::max();
2913  size_type itr; // move to somewhere else?
2914 
2915  // Put these here because they need to be remembered between iterations if a linearly
2916  // dependent constriant is dropped.
2917  size_type ja = 0; // + indice of violated constraint to add to active set
2918  size_type last_ja = 0; // + last_ja to be added
2919  value_type con_ja_val; // value of violated constraint.
2920  value_type b_a; // value of the violated bound
2921  value_type norm_2_constr; // norm of violated constraint
2922  EBounds bnd_ja; // bound of constraint ja which is violated.
2923  bool can_ignore_ja; // true if we can ignore a constraint if it is LD.
2924  bool assume_lin_dep_ja;
2925  value_type gamma_plus; // used to store the new multipler value for the added
2926  // constraint.
2927  const int summary_lines_counter_max = 15;
2928  int summary_lines_counter = 0;
2929  long int jd = 0; // + indice of constraint to delete from active set.
2930  // - indice of intially fixed variable to be freed
2931  long int last_jd = 0; // Last jd change to the active set
2932  value_type t_P; // Primal step length (constraint ja made active)
2933  value_type t_D; // Dual step length ( longest step without violating dual
2934  // feasibility of currently active constraints ).
2935  value_type beta; // +1 if multiplier of constraint being added is positive
2936  // -1 if multiplier of constraint being added is negative.
2937  bool warned_degeneracy = false; // Will warn the user if degeneracy
2938  // is detected.
2939  value_type dual_infeas_scale = 1.0; // Scaling for determining if a
2940  // Lagrange multiplier is near degenerate.
2941  bool return_to_init_fixed = false; // True if the constraint being added
2942  // to the active set is a variable returning to its orginally
2943  // fixed variable bound.
2944  bool using_iter_refinement = false; // Will be set to true if instability detected
2945 
2946  for( itr = 0; itr <= max_iter_; ++itr, ++(*iter) ) {
2947  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2948  *out
2949  << "\n************************************************"
2950  << "\n*** qp_iter = " << itr
2951  << "\n*** q_hat = " << act_set->q_hat() << std::endl;
2952  }
2953  bool schur_comp_update_failed = false;
2954  switch( next_step ) { // no break; statements in this switch statement.
2955  case PICK_VIOLATED_CONSTRAINT: {
2956  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2957  *out
2958  << "\n*** PICK_VIOLATED_CONSTRAINT\n";
2959  }
2960  // Check runtime
2961  if( timeout_return(timer,out,output_level) )
2962  return MAX_RUNTIME_EXEEDED_DUAL_FEAS;
2963  // Save the indice of the last constriant to be added!
2964  last_ja = ja;
2965  // Set parts of x that are not currently fixed and may have changed.
2966  // Also, we want set specifially set those variables that where
2967  // initially free and then latter fixed to their bounds so that
2968  // they will not be seen as violated.
2969  set_x( *act_set, *v, x );
2970  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2971  *out
2972  << "\n||x||inf = " << norm_inf(*x) << std::endl;
2973  }
2974  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
2975  *out
2976  << "\nx =\n" << *x;
2977  }
2978  if( test_what == RUN_TESTS ) {
2979  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2980  *out
2981  << "\nChecking current iteration quantities ...\n";
2982  }
2983  const char
2984  sep_line[] = "\n--------------------------------------------------------------------------------\n";
2985 // AbstractLinAlgPack::TestingPack::CompareDenseVectors comp_v;
2986 
2987  //
2988  // Check the optimality conditions of the augmented KKT system!
2989  //
2990 
2991  // ToDo: Implement
2992 
2993  //
2994  // Check mu_D_hat
2995  //
2996  if( act_set->q_D_hat() ) {
2997  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
2998  *out
2999  << sep_line
3000  << "\nComputing: mu_D_hat_calc = -Q_XD_hat'*g - Q_XD_hat'*G*x\n"
3001  << " - Q_XD_hat'*A*v(n_R+1:n_R+m) - Q_XD_hat'*A_bar*P_plus_hat*z_hat ...\n";
3002  }
3003  Workspace<value_type> mu_D_hat_calc_ws( wss, act_set->q_D_hat() );
3004  DVectorSlice mu_D_hat_calc( &mu_D_hat_calc_ws[0], mu_D_hat_calc_ws.size() );
3005  calc_mu_D( *act_set, *x, *v, &mu_D_hat_calc );
3006  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3007  *out
3008  << "\n||mu_D_hat_calc||inf = " << norm_inf(mu_D_hat_calc) << std::endl;
3009  }
3010  if( (int)output_level >= (int)OUTPUT_ACT_SET ) {
3011  *out
3012  << "\nmu_D_hat_calc =\n" << mu_D_hat_calc;
3013  }
3014  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3015  *out
3016  << "\nChecking mu_D_hat_calc == mu_D_hat\n";
3017  }
3018  DVector mu_D_hat_diff(mu_D_hat_calc.dim());
3019  LinAlgOpPack::V_VmV( &mu_D_hat_diff(), mu_D_hat_calc(), act_set->mu_D_hat() );
3020  const value_type
3021  mu_D_hat_err = norm_inf(mu_D_hat_diff()) / (1.0 + norm_inf(mu_D_hat_calc()));
3022  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3023  *out
3024  << "\n||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3025  << mu_D_hat_err << std::endl;
3026  }
3028  mu_D_hat_err >= error_tol(), TestFailed
3029  ,"QPSchur::qp_algo(...) : Error, "
3030  "||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3031  << mu_D_hat_err << " >= error_tol = " << error_tol()
3032  );
3033  if( mu_D_hat_err >= warning_tol() && (int)output_level >= (int)OUTPUT_ACT_SET ) {
3034  *out
3035  << "\nWarning! ||mu_D_hat_calc-mu_D_hat||inf/(1.0+||mu_D_hat_calc||inf) = "
3036  << mu_D_hat_err << " >= warning_tol = " << warning_tol() << std::endl;
3037  }
3038  }
3039  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3040  *out
3041  << sep_line;
3042  }
3043  }
3044  act_set->qp().constraints().pick_violated( *x, &ja, &con_ja_val
3045  , &b_a, &norm_2_constr, &bnd_ja, &can_ignore_ja );
3046  assume_lin_dep_ja = false; // Assume this initially.
3047  if( ja > 0 && act_set->is_init_fixed(ja) && qp.x_init()(ja) == bnd_ja )
3048  return_to_init_fixed = true;
3049  else
3050  return_to_init_fixed = false;
3051  // Print ja, bnd_ja, can_ignore_ja
3052  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3053  *out
3054  << "\nja = " << ja << endl;
3055  if(ja) {
3056  *out
3057  << "\ncon_ja_val = " << con_ja_val
3058  << "\nb_a = " << b_a
3059  << "\nnorm_2_constr = " << norm_2_constr
3060  << "\nbnd_ja = " << bnd_str(bnd_ja)
3061  << "\ncan_ignore_ja = " << bool_str(can_ignore_ja)
3062  << "\nreturn_to_init_fixed = " << bool_str(return_to_init_fixed)
3063  << endl
3064  ;
3065  }
3066  }
3067  // Print header for itr, nact, change (ADD, DROP), indice (ja, jb)
3068  //, bound (LOWER, UPPER, EQUALITY), violation (primal or dual), rank (LD,LI)
3069  if( (int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
3070  if( summary_lines_counter <= 0 ) {
3071  summary_lines_counter = summary_lines_counter_max;
3072  *out
3073  << endl
3074  << right << setw(6) << "itr"
3075  << right << setw(6) << "qhat"
3076  << right << setw(6) << "q(+)"
3077  << right << setw(6) << "q_F"
3078  << right << setw(6) << "q_C"
3079  << right << setw(6) << "q_D"
3080  << right << setw(8) << "change"
3081  << right << setw(9) << "type"
3082  << right << setw(6) << "j"
3083  << right << setw(10) << "bnd"
3084  << right << setw(dbl_w) << "viol, p_z(jd)"
3085  << right << setw(6) << "rank" << endl
3086  << right << setw(6) << "----"
3087  << right << setw(6) << "----"
3088  << right << setw(6) << "----"
3089  << right << setw(6) << "----"
3090  << right << setw(6) << "----"
3091  << right << setw(6) << "----"
3092  << right << setw(8) << "------"
3093  << right << setw(9) << "-------"
3094  << right << setw(6) << "----"
3095  << right << setw(10) << "--------"
3096  << right << setw(dbl_w) << "--------------"
3097  << right << setw(6) << "----" << endl;
3098  }
3099  }
3100  // Print first part of row for itr, q_hat, q(+), q_D, q_C, q_F, change, type, index, bound, violation
3101  if( (int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
3102  *out
3103  << right << setw(6) << itr // itr
3104  << right << setw(6) << act_set->q_hat() // q_hat
3105  << right << setw(6) << act_set->q_plus_hat() // q(+)
3106  << right << setw(6) << act_set->q_F_hat() // q_F
3107  << right << setw(6) << act_set->q_C_hat() // q_C
3108  << right << setw(6) << act_set->q_D_hat() // q_D
3109  << right << setw(8) << ( ja ? "ADD" : "-" ) // change
3110  << right << setw(9); // type
3111  if( ja == 0 ) {
3112  *out << "-";
3113  }
3114  else if( act_set->is_init_fixed(ja) ) {
3115  if( bnd_ja == qp.x_init()(ja) )
3116  *out << "X_F_X";
3117  else
3118  *out << "X_F_C";
3119  }
3120  else if( ja <= n ) {
3121  *out << "R_X";
3122  }
3123  else {
3124  *out << "GEN";
3125  }
3126  *out
3127  << right << setw(6) << ja // index
3128  << right << setw(10) << ( ja ? bnd_str(bnd_ja) : "-" ) // bound
3129  << right << setw(dbl_w); // violation
3130  if(ja)
3131  *out << (con_ja_val - b_a);
3132  else
3133  *out << "-";
3134  if(!ja)
3135  *out << right << setw(6) << "-" << endl; // rank for last iteration
3136  }
3137  bool found_solution = false;
3138  const size_type sa = act_set->s_map(ja);
3139  value_type scaled_viol = 0.0;
3140  if( ja == 0 ) {
3141  found_solution = true;
3142  }
3143  else if( sa != 0 || ( act_set->is_init_fixed(ja) && act_set->s_map(-ja) == 0 ) ) {
3144  const bool is_most_violated
3145  = (qp.constraints().pick_violated_policy() == QPSchurPack::Constraints::MOST_VIOLATED);
3146  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3147  *out
3148  << "\n\nWarning, we have picked the constriant a(" << ja << ") with violation\n"
3149  << "(a(ja)'*x - b_a) = (" << con_ja_val << " - " << b_a << ") = " << (con_ja_val - b_a)
3150  << "\nto add to the active set but it is already part of the active set.\n";
3151  }
3152  const EBounds act_bnd = ( sa != 0 ? act_set->bnd(sa) : qp.x_init()(ja) );
3153  if( act_bnd != bnd_ja ) {
3154  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3155  *out
3156  << "However, this is not the same bound as the active bound!\n";
3157  }
3158  const value_type
3159  act_b_a = qp.constraints().get_bnd(ja,act_bnd);
3160  if( act_bnd == LOWER && act_b_a > b_a || act_bnd == UPPER && act_b_a < b_a ) {
3161  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3162  *out
3163  << "\nError, c_L_bar(" << ja <<") = " << (act_b_a > b_a ? act_b_a : b_a)
3164  << " > c_U_bar(" << ja << ") = " << (act_b_a < b_a ? act_b_a : b_a) << ".\n";
3165  }
3166  }
3167  else {
3168  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
3169  }
3170  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3171  *out
3172  << "The constraints are infeasible! Terminating the QP algorithm!\n";
3173  }
3174  return INFEASIBLE_CONSTRAINTS;
3175  }
3176  else {
3177  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3178  *out
3179  << "This is the active bound so this is an indication of instability\n"
3180  << "in the calculations.\n";
3181  }
3182  }
3183  summary_lines_counter = 0;
3184  if( !is_most_violated ) {
3185  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3186  *out
3187  << "\nThis is not the most violated constraint so set\n"
3188  << "pick_violated_policy = MOST_VIOLATED ...\n";
3189  }
3190  summary_lines_counter = 0;
3191  qp.constraints().pick_violated_policy(QP::Constraints::MOST_VIOLATED);
3192  next_step = PICK_VIOLATED_CONSTRAINT;
3193  continue; // Go back and pick the most violated constraint
3194  }
3195  else if( !using_iter_refinement ) {
3196  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3197  *out
3198  << "\nThis is the most violated constraint but we are currently not using\n"
3199  << "iterative refinement so let's switch it on ...\n";
3200  }
3201  summary_lines_counter = 0;
3202  EIterRefineReturn status = iter_refine(
3203  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, act_set->q_hat() ? &act_set->d_hat() : NULL
3204  ,v, act_set->q_hat() ? &act_set->z_hat() : NULL
3205  ,iter_refine_num_resid, iter_refine_num_solves
3206  );
3207  if( status == ITER_REFINE_IMPROVED || status == ITER_REFINE_CONVERGED ) {
3208  using_iter_refinement = true; // Use iterative refinement from now on!
3209  next_step = PICK_VIOLATED_CONSTRAINT;
3210  continue; // Now go back
3211  }
3212  else {
3213  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3214  *out
3215  << "\nError, iterative refinement was not allowed or failed to improve the residuals.\n"
3216  << "Terminating the QP algorithm ...\n";
3217  }
3218  return SUBOPTIMAL_POINT;
3219  }
3220  }
3221  else {
3222  scaled_viol = std::fabs(con_ja_val - b_a)/(1.0 + std::fabs(con_ja_val));
3223  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3224  *out
3225  << "\n\nThis is the most violated constraint, we are using iterative refinement and\n"
3226  << "|a("<<ja<<")'*x - b_a| / (1 + |a("<<ja<<")'*x|) = |"<<con_ja_val<<" - "<<b_a
3227  << "| / (1 + |"<<con_ja_val<<"|) = "<<scaled_viol<< (scaled_viol<loose_feas_tol()?" < ":" > ")
3228  << "loose_feas_tol = "<<loose_feas_tol();
3229  }
3230  if( scaled_viol < loose_feas_tol() ) {
3231  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3232  *out
3233  << "\nTerminating the algorithm with a near optimal solution!\n";
3234  }
3235  found_solution = true;
3236  }
3237  else {
3238  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3239  *out
3240  << "\nError! The QP algorithm is terminated!\n";
3241  }
3242  return SUBOPTIMAL_POINT;
3243  }
3244  }
3245  }
3246  else if( act_set->all_dof_used_up()
3247  && (scaled_viol = std::fabs(con_ja_val - b_a)/(1.0 + std::fabs(con_ja_val))) < feas_tol() )
3248  {
3249  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3250  *out
3251  << "\nWith all the dof used up, the violated inequality constriant "
3252  << "a("<< ja << ")'*x statisfies:\n"
3253  << "|a("<<ja<<")'*x - b_a| / (1 + |a("<<ja<<")'*x|) = |" << con_ja_val << " - " << b_a
3254  << "| / (1 + |"<<con_ja_val<<"|) = "<<scaled_viol<<" < feas_tol = "<<feas_tol();
3255  }
3256  if( act_set->qp().constraints().pick_violated_policy() == QP::Constraints::MOST_VIOLATED ) {
3257  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3258  *out
3259  << "\nThis is the most violated constraint so we will consider this\n"
3260  << "a near degenerate constraint so we are done!\n";
3261  }
3262  found_solution = true;
3263  }
3264  else {
3265  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3266  *out
3267  << "\nThis is not the most violated constraint so set pick_violated_policy = "
3268  << "MOST_VIOLATED and pick another violated constraint ....\n";
3269  }
3270  act_set->qp().constraints().pick_violated_policy(QP::Constraints::MOST_VIOLATED);
3271  next_step = PICK_VIOLATED_CONSTRAINT;
3272  continue; // Go back and pick the most violated constraint
3273  }
3274  }
3275  if(found_solution) {
3276  //
3277  // Solution found! All of the inequality constraints are satisfied!
3278  //
3279  // Use iterative refinement if we need to
3280  if( iter_refine_at_solution() && !using_iter_refinement ) {
3281  // The user has requested iterative refinement at the solution
3282  // and we have not been using iterative refinement up to this point.
3283  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3284  *out
3285  << "\nWe think we have found the solution and are not currently using iterative refinement\n"
3286  << "and iter_refine_at_solution==true so perform iterative refinement ...\n";
3287  }
3288  using_iter_refinement = true;
3289  EIterRefineReturn status = iter_refine(
3290  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, act_set->q_hat() ? &act_set->d_hat() : NULL
3291  ,v, act_set->q_hat() ? &act_set->z_hat() : NULL
3292  ,iter_refine_num_resid, iter_refine_num_solves
3293  );
3294  switch(status) {
3295  case ITER_REFINE_ONE_STEP:
3296  case ITER_REFINE_IMPROVED:
3297  case ITER_REFINE_CONVERGED:
3298  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3299  *out
3300  << "\nIterative refinement may have altered the unknowns so go back and look for another violated constraint ...\n";
3301  }
3302  summary_lines_counter = 0;
3303  next_step = PICK_VIOLATED_CONSTRAINT;
3304  continue;
3305  case ITER_REFINE_NOT_PERFORMED:
3306  case ITER_REFINE_NOT_NEEDED:
3307  case ITER_REFINE_NOT_IMPROVED:
3308  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3309  *out
3310  << "\nIterative refinement did not alter the unknowns so exit with this solution...\n";
3311  }
3312  set_x( *act_set, *v, x );
3313  break;
3314  default:
3315  TEUCHOS_TEST_FOR_EXCEPT(true); // Local programming error only!
3316  }
3317  }
3318  if( iter_refine_at_solution() || using_iter_refinement ) {
3319  // The user has requested iterative refinement at the solution
3320  // or we have been using iterative refinement so we should
3321  // compute the lagrange multipliers for the initially fixed
3322  // variables that are still fixed.
3323  if( act_set->q_D_hat() ) {
3324  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3325  *out
3326  << "\nRecomputing final mu_D_hat at the solution ...\n";
3327  }
3328  calc_mu_D( *act_set, *x, *v, &act_set->mu_D_hat() );
3329  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3330  *out
3331  << "\n||mu_D_hat||inf = " << norm_inf(act_set->mu_D_hat()) << std::endl;
3332  }
3333  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3334  *out
3335  << "\nmu_D_hat =\n" << act_set->mu_D_hat();
3336  }
3337  }
3338  }
3339  return OPTIMAL_SOLUTION; // current point is optimal.
3340  }
3341  }
3342  case UPDATE_ACTIVE_SET: {
3343  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3344  *out
3345  << "\n*** UPDATE_ACTIVE_SET\n";
3346  }
3347  ++(*num_adds);
3348  if( act_set->all_dof_used_up() || act_set->is_init_fixed(ja) ) {
3349  // All of the degrees of freedom are currently used up so we must
3350  // assume that we must remove one of these currently active
3351  // constraints and replace it with the violated constraint.
3352  // In this case we know that this will make the schur
3353  // complement singular so let's just skip the update
3354  // and set this now. We also may be here if we are fixing
3355  // an initially fixed variable to some bound.
3356  assume_lin_dep_ja = true;
3357  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3358  if(act_set->all_dof_used_up()) {
3359  *out
3360  << "\nAll of the degrees of freedom are used up so "
3361  << "the constraint ja must be linearly dependent.\n"
3362  << "The augmented KKT system is definitely singular!\n";
3363  }
3364  else {
3365  *out
3366  << "\nThis is an initially fixed variable that was freed and "
3367  << "now is being fixed again.\n"
3368  << "The augmented KKT system could be singular or nonsingular, "
3369  << "we don't know at this point.\n";
3370  }
3371  *out
3372  << "\nThe schur complement for the new KKT system will not be "
3373  << "updated yet in case it is singular!\n";
3374  }
3375  }
3376  else {
3377  assume_lin_dep_ja = false;
3378  try {
3379  if(act_set->add_constraint( ja, bnd_ja, false, out, output_level, true ))
3380  summary_lines_counter = 0;
3381  else {
3382  // Print end of row for rank if the right print level
3383  if( (int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
3384  *out << right << setw(6) << "LI" << endl;
3385  out->flush();
3386  --summary_lines_counter;
3387  }
3388  }
3389  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3390  *out << "\nNew KKT system is nonsingular! (linearly independent (LI) constraints)\n";
3391  }
3392  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3393  *out
3394  << "\nPrinting active set after adding constraint ja = " << ja
3395  << " ...\n";
3396  dump_act_set_quantities( *act_set, *out );
3397  }
3398  }
3399  catch( const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
3400  // Constraint really is linearly dependent.
3401  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3402  *out
3403  << "\n\nSchur complement update failed:\n"
3404  << "(" << excpt.what() << ")\n"
3405  << "\nConstraint ja = " << ja << " appears to be linearly dependent!\n\n";
3406  }
3407  summary_lines_counter = 0;
3408  if( !(act_set->q_D_hat() + act_set->q_plus_hat()) ) {
3409  // Print omsg.
3410  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3411  *out
3412  << "\nQPSchur::qp_algo(...) : "
3413  << "Error, constraint j = "<< ja << " is linearly dependent "
3414  << "and there are no other constraints to drop.\n"
3415  << "The QP must be infeasible\n";
3416  }
3417  return INFEASIBLE_CONSTRAINTS;
3418  }
3419  assume_lin_dep_ja = true;
3420  schur_comp_update_failed = true;
3421  }
3422  catch( const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
3423  // Reduced Hessian has the wrong inertia
3424  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3425  *out
3426  << "\n\nSchur complement appears to have the wrong inertia:\n"
3427  << "(" << excpt.what() << ")\n"
3428  << "\nThe QP appears to be nonconvex.\n"
3429  << "\nWe have no choice but to terminate the primal-dual QP algorithm!\n";
3430  }
3431  return NONCONVEX_QP;
3432  }
3433  }
3434  if( assume_lin_dep_ja && can_ignore_ja ) {
3435  act_set->qp().constraints().ignore( ja );
3436  next_step = PICK_VIOLATED_CONSTRAINT;
3437  continue;
3438  }
3439  // Infer the sign of the multiplier for the new constraint being added
3440  beta = ( con_ja_val > b_a ? +1.0 : -1.0 );
3441  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3442  *out
3443  << "\nbeta = " << beta << endl;
3444  }
3445  }
3446  case COMPUTE_SEARCH_DIRECTION: {
3447  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3448  *out
3449  << "\n*** COMPUTE_SEARCH_DIRECTION\n";
3450  }
3451  const EtaVector e_ja( ja, n + m_breve );
3452  if( assume_lin_dep_ja ) {
3453  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3454  *out
3455  << "\nThe KKT system for the trial active set is assumed or known to be singular!\n"
3456  << "\nComputing the steps from:"
3457  << "\np_z_hat = inv(S_hat)*(-v_a+U_hat'*inv(Ko)*u_a), p_v = inv(Ko)*(-u_a-U_hat*p_z_hat) ...\n";
3458  }
3459  //
3460  // The schur complement is not updated so we must compute
3461  // p_z_hat and p_v explicitly.
3462  //
3463  // If all the degrees of freedom are used up then we know that the step for
3464  // the primal variables will be zero. However, if m > 0 then v and p_v also
3465  // contain terms for the Lagrange multipliers for the equality constriants
3466  // but we don't need to compute these during the algorithm.
3467  // Therefore we can just set p_v = 0 and save a solve with Ko.
3468  // If the QP is feasible then a constraint will be dropped, the
3469  // KKT system will be updated and then v_plus will be computed
3470  // at the next iteration and be used to compute v so all is good.
3471  //
3472  const bool all_dof_used_up = act_set->all_dof_used_up();
3473  if( act_set->is_init_fixed(ja) ) {
3474  //
3475  // Fix a varaible that was fixed and then freed.
3476  //
3477  // [ Ko U_hat ] [ p_v ] = [ 0 ]
3478  // [ U_hat' V_hat ] [ p_z_hat ] [ -v_a ]
3479  //
3480  // v_a = e(sa) <: R^q_hat : where sa = s_map(-ja)
3481  //
3482  // / 0 : if x_init(ja) == bnd_ja
3483  // d_a = |
3484  // \ b_a - b_X(l_x_X_map(ja)) : otherwise
3485  //
3486  // p_z_hat = -inv(S_hat)*v_a
3487  //
3488  // p_v = inv(Ko)*(-U_hat*p_z_hat)
3489  //
3490  // gamma_plus = (d_a - v_a' * z_hat ) / ( v_a' * p_z_hat )
3491  //
3492  const size_type
3493  sa = act_set->s_map(-int(ja)),
3494  la = act_set->qp().l_x_X_map()(ja);
3495  TEUCHOS_TEST_FOR_EXCEPT( !( sa ) );
3496  TEUCHOS_TEST_FOR_EXCEPT( !( la ) );
3497  // v_a = e(sa) <: R^q_hat
3498  Workspace<value_type> v_a_ws(wss,act_set->q_hat());
3499  DVectorSlice v_a(&v_a_ws[0],v_a_ws.size());
3500  v_a = 0.0;
3501  v_a(sa) = 1.0;
3502  // d_a
3503  const value_type
3504  d_a = ( bnd_ja == qp.x_init()(ja)
3505  ? 0.0
3506  : b_a - qp.b_X()(la) );
3507  // p_z_hat = -inv(S_hat) * v_a
3508  V_InvMtV( &act_set->p_z_hat(), act_set->S_hat(), no_trans, v_a() );
3509  Vt_S( &act_set->p_z_hat(), -1.0 );
3510  // p_v = inv(Ko)*(-U_hat*p_z_hat)
3511  if(!all_dof_used_up) {
3512  calc_v( qp.Ko(), NULL, act_set->U_hat(), act_set->p_z_hat(), &p_v() );
3513  }
3514  else {
3515  p_v = 0.0;
3516  }
3517  // Iterative refinement?
3518  if( using_iter_refinement ) {
3519  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3520  *out
3521  << "\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3522  }
3523  summary_lines_counter = 0;
3524  // [ Ko U_hat ] [ p_v ] = [ 0 ]
3525  // [ U_hat' V_hat ] [ p_z_hat ] [ -v_a ]
3526  EIterRefineReturn status = iter_refine(
3527  *act_set, out, output_level, 0.0, NULL, +1.0, &v_a
3528  ,&p_v, &act_set->p_z_hat()
3529  ,iter_refine_num_resid, iter_refine_num_solves
3530  );
3531  }
3532  // gamma_plus = ( d_a - v_a'*z_hat ) / ( v_a'*p_z_hat )
3533  if(!all_dof_used_up)
3534  gamma_plus = ( ( d_a - dot(v_a(),act_set->z_hat()) )
3535  / ( dot(v_a(),act_set->p_z_hat()) ) );
3536  else
3537  gamma_plus = beta * inf;
3538  }
3539  else {
3540  //
3541  // Add a constraint that is not an initially fixed
3542  // variable bound.
3543  //
3544  // [ Ko U_hat ] [ p_v ] = [ -u_a ]
3545  // [ U_hat' V_hat ] [ p_z_hat ] [ -v_a ]
3546  //
3547  // p_z_hat = inv(S_hat) * ( - v_a + U_hat' * inv(Ko) * u_a )
3548  //
3549  // p_v = inv(Ko) * ( -u_a - U_hat * p_z_hat )
3550  //
3551  // gamma_plus = ( d_a - u_a'*v - v_a'*z_hat ) / ( u_a'*p_v + v_a'*p_z_hat )
3552  //
3553  // ToDo: (9/25/00): Make u_a and v_a both sparse and combine the following code.
3554  //
3555  if( ja <= n ) {
3556  // Fix an initially free variable
3557  //
3558  // u_a = [ Q_R' * e(ja) ] <: R^(n_R+m)
3559  // [ 0 ]
3560  //
3561  // v_a = 0 <: R^(q_hat)
3562  //
3563  // d_a = b_a <: R
3564  //
3565  const size_type
3566  la = act_set->qp().Q_R().lookup_col_j(ja);
3567  TEUCHOS_TEST_FOR_EXCEPT( !( la ) );
3568  const EtaVector u_a = EtaVector(la,n_R+m);
3569  const value_type d_a = b_a;
3570  DVector t1;
3571  // t1 = inv(Ko) * u_a
3572  V_InvMtV( &t1, qp.Ko(), no_trans, u_a() );
3573  if( act_set->q_hat() ) {
3574  // t2 = U_hat'*t1
3575  DVector t2;
3576  V_MtV( &t2, act_set->U_hat(), trans, t1() );
3577  // p_z_hat = inv(S_hat) * t2
3578  V_InvMtV( &act_set->p_z_hat(), act_set->S_hat(), no_trans, t2() );
3579  // t1 = - u_a
3580  V_StV( &t1, -1.0, u_a() );
3581  // t1 += - U_hat * p_z_hat
3582  Vp_StMtV( &t1(), -1.0, act_set->U_hat(), no_trans, act_set->p_z_hat() );
3583  // p_v = inv(Ko) * t1
3584  if(!all_dof_used_up)
3585  V_InvMtV( &p_v, qp.Ko(), no_trans, t1() );
3586  else
3587  p_v = 0.0;
3588  }
3589  else {
3590  // p_v = -t1
3591  V_mV( &p_v, t1() );
3592  }
3593  // Iterative refinement?
3594  if( using_iter_refinement ) {
3595  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3596  *out
3597  << "\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3598  }
3599  summary_lines_counter = 0;
3600  // [ Ko U_hat ] [ p_v ] = [ -u_a ]
3601  // [ U_hat' V_hat ] [ p_z_hat ] [ 0 ]
3602  Workspace<value_type> dense_u_a_ws(wss,u_a().dim());
3603  DVectorSlice dense_u_a(&dense_u_a_ws[0],dense_u_a_ws.size());
3604  dense_u_a = 0.0; // Make a dense copy of u_a!
3605  dense_u_a(u_a().begin()->index()+u_a().offset()) = 1.0;
3606  EIterRefineReturn status = iter_refine(
3607  *act_set, out, output_level, +1.0, &dense_u_a, 0.0, NULL
3608  ,&p_v, &act_set->p_z_hat()
3609  ,iter_refine_num_resid, iter_refine_num_solves
3610  );
3611  summary_lines_counter = 0;
3612  }
3613  // gamma_plus = ( d_a - u_a'*v) / ( u_a'*p_v )
3614  if(!all_dof_used_up)
3615  gamma_plus = ( d_a - dot(u_a(),*v) ) / dot(u_a(),p_v());
3616  else
3617  gamma_plus = beta * inf;
3618  }
3619  else {
3620  // Add a general inequality (or equality) constraint
3621  //
3622  // u_a = [ Q_R' * A_bar * e(ja) ] <: R^(n_R + m)
3623  // [ 0 ]
3624  //
3625  // v_a = P_XF_hat' * A_bar * e_ja <: R^(q_hat)
3626  //
3627  // d_a = b_a - b_X' * (Q_X' * A_bar * e_ja) <: R
3628  //
3629  Workspace<value_type> u_a_ws( wss, n_R + m );
3630  DVectorSlice u_a( &u_a_ws[0], u_a_ws.size() );
3631  Workspace<value_type> v_a_ws( wss, act_set->q_hat() );
3632  DVectorSlice v_a( &v_a_ws[0], v_a_ws.size() );
3633  // u_a(1:n_R) = Q_R' * A_bar * e(ja)
3634  Vp_StPtMtV( &u_a(1,n_R), 1.0, qp.Q_R(), trans
3635  , qp.constraints().A_bar(), no_trans, e_ja(), 0.0 );
3636  // u_a(n_R+1:n_R+m) = 0.0
3637  if(m)
3638  u_a(n_R+1,n_R+m) = 0.0;
3639  // t0 = Q_X' * A_bar * e_ja
3640  Workspace<value_type> t0_ws( wss, n-n_R );
3641  DVectorSlice t0( &t0_ws[0], t0_ws.size() );
3642  if( n > n_R )
3643  Vp_StPtMtV( &t0(), 1.0, qp.Q_X(), trans
3644  , qp.constraints().A_bar(), no_trans, e_ja(), 0.0 );
3645  // d_a = b_a - b_X'*t0
3646  const value_type
3647  d_a = b_a - ( n > n_R ? dot( qp.b_X(), t0() ) : 0.0 );
3648  // t1 = inv(Ko) * u_a
3649  Workspace<value_type> t1_ws( wss, n_R + m );
3650  DVectorSlice t1( &t1_ws[0], t1_ws.size() );
3651  V_InvMtV( &t1, qp.Ko(), no_trans, u_a );
3652  if( act_set->q_hat() ) {
3653  // t2 = U_hat'*t1
3654  Workspace<value_type> t2_ws( wss, act_set->q_hat() );
3655  DVectorSlice t2( &t2_ws[0], t2_ws.size() );
3656  V_MtV( &t2, act_set->U_hat(), trans, t1() );
3657  // v_a = P_XF_hat' * A_bar * e_ja
3658  Vp_StPtMtV( &v_a(), 1.0, act_set->P_XF_hat(), trans
3659  , qp.constraints().A_bar(), no_trans, e_ja(), 0.0 );
3660  // t2 += -v_a
3661  Vp_StV( &t2(), -1.0, v_a() );
3662  // p_z_hat = inv(S_hat) * t2
3663  V_InvMtV( &act_set->p_z_hat(), act_set->S_hat(), no_trans, t2() );
3664  if(!all_dof_used_up) {
3665  // t1 = - u_a
3666  V_StV( &t1, -1.0, u_a() );
3667  // t1 += - U_hat * p_z_hat
3668  Vp_StMtV( &t1(), -1.0, act_set->U_hat(), no_trans, act_set->p_z_hat() );
3669  // p_v = inv(Ko) * t1
3670  V_InvMtV( &p_v, qp.Ko(), no_trans, t1() );
3671  }
3672  else {
3673  p_v = 0.0;
3674  }
3675  // Iterative refinement?
3676  if( using_iter_refinement ) {
3677  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3678  *out
3679  << "\n\nPerforming iterative refinement on p_v, p_z_hat system ...\n";
3680  }
3681  summary_lines_counter = 0;
3682  // [ Ko U_hat ] [ p_v ] = [ -u_a ]
3683  // [ U_hat' V_hat ] [ p_z_hat ] [ -v_a ]
3684  EIterRefineReturn status = iter_refine(
3685  *act_set, out, output_level, +1.0, &u_a, +1.0, act_set->q_hat() ? &v_a : NULL
3686  ,&p_v, act_set->q_hat() ? &act_set->p_z_hat() : NULL
3687  ,iter_refine_num_resid, iter_refine_num_solves
3688  );
3689  summary_lines_counter = 0;
3690  }
3691  // gamma_plus = ( d_a - u_a'*v - v_a'*z_hat ) / ( u_a'*p_v + v_a'*p_z_hat )
3692  if(!all_dof_used_up)
3693  gamma_plus = ( ( d_a - dot(u_a,*v) - dot(v_a(),act_set->z_hat()) )
3694  / ( dot(u_a,p_v()) + dot(v_a(),act_set->p_z_hat()) ) );
3695  else
3696  gamma_plus = beta * inf;
3697  }
3698  else {
3699  // p_v = -t1
3700  if(!all_dof_used_up)
3701  V_mV( &p_v, t1() );
3702  else
3703  p_v = 0.0;
3704  // gamma_plus = ( d_a - u_a'*v) / ( u_a'*p_v )
3705  if(!all_dof_used_up)
3706  gamma_plus = ( d_a - dot(u_a,*v) ) / dot(u_a,p_v());
3707  else
3708  gamma_plus = beta * inf;
3709  }
3710  }
3711  }
3712  if( schur_comp_update_failed && gamma_plus * beta < 0 ) {
3713  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3714  *out
3715  << "\nThe schur complement update failed and gamma_plus = " << gamma_plus << " is the wrong sign"
3716  << "\nso we will assume the sign error for (...)/+-0 was due to arbitrary roundoff"
3717  << "\nand therefore we will set gamma_plus = -gamma_plus\n";
3718  }
3719  gamma_plus = -gamma_plus;
3720  }
3721  // Print the steps p_v and p_z_hat
3722  if(act_set->q_hat()) {
3723  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3724  *out
3725  << "\n||p_z_hat||inf = " << norm_inf(act_set->p_z_hat()) << endl;
3726  }
3727  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3728  *out << "\np_z_hat =\n" << act_set->p_z_hat();
3729  }
3730  }
3731  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3732  *out
3733  << "\n||p_v||inf = " << norm_inf(p_v()) << endl;
3734  }
3735  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3736  *out << "\np_v =\n" << p_v();
3737  }
3738  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3739  *out
3740  << "\ngamma_plus = " << gamma_plus << endl;
3741  }
3742  // Compute step for mu_D_hat
3743  if( act_set->q_D_hat() ) {
3744  // Compute for steps of all the constraints in the current active set
3745  calc_p_mu_D( *act_set, p_v(), act_set->p_z_hat(), &ja, &act_set->p_mu_D_hat() );
3746  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3747  *out
3748  << "\n||p_mu_D_hat||inf = " << norm_inf(act_set->p_mu_D_hat()) << std::endl;
3749  }
3750  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3751  *out
3752  << "\np_mu_D_hat =\n" << act_set->p_mu_D_hat();
3753  }
3754  }
3755  }
3756  else {
3757  // The new schur complement is already updated so compute
3758  // the solution outright.
3759  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3760  *out
3761  << "\nThe KKT system for the new active set is or known to be nonsingular!\n"
3762  << "\nComputing the steps from:"
3763  << "\nz_hat_plus = inv(S_hat)*(d_hat-U_hat'vo), v_plus = inv(Ko)*(fo-U_hat*z_hat_plus) ...\n";
3764  }
3765 
3766  // Compute z_hat_plus, v_plus
3767 
3768  // z_hat_plus = inv(S_hat) * ( d_hat - U_hat' * vo )
3769  const size_type q_hat = act_set->q_hat();
3770  z_hat_plus.bind(DVectorSlice(&z_hat_plus_ws[0],q_hat));
3771  calc_z( act_set->S_hat(), act_set->d_hat(), act_set->U_hat(), &vo
3772  , &z_hat_plus );
3773  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3774  *out
3775  << "\n||z_hat_plus||inf = " << norm_inf(z_hat_plus()) << std::endl;
3776  }
3777  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3778  *out
3779  << "\nz_hat_plus =\n" << z_hat_plus();
3780  }
3781  // v_plus = inv(Ko) * (fo - U_hat * z_hat_plus)
3782  calc_v( qp.Ko(), &qp.fo(), act_set->U_hat(), z_hat_plus
3783  , &v_plus );
3784  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3785  *out
3786  << "\n||v_plus||inf = " << norm_inf(v_plus()) << std::endl;
3787  }
3788  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3789  *out
3790  << "\nv_plus =\n" << v_plus();
3791  }
3792  if( using_iter_refinement ) {
3793  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
3794  *out
3795  << "\n\nPerforming iterative refinement on v_plus, z_hat_plus system ...\n";
3796  }
3797  summary_lines_counter = 0;
3798  // [ Ko U_hat ] [ v_plus ] = [ fo ]
3799  // [ U_hat' V_hat ] [ z_hat_plus ] [ d_hat ]
3800  EIterRefineReturn status = iter_refine(
3801  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, &act_set->d_hat()
3802  ,&v_plus, &z_hat_plus
3803  ,iter_refine_num_resid, iter_refine_num_solves
3804  );
3805  }
3806  // Compute p_z_hat (change in z_hat w.r.t newly added constriant multiplier)
3807  DVectorSlice p_z_hat = act_set->p_z_hat();
3808  // p_z_hat = z_hat_plus - z_hat
3809  V_VmV( &p_z_hat(), z_hat_plus(), act_set->z_hat() );
3810  // p_v = v_plus - v
3811  V_VmV( &p_v(), v_plus(), *v );
3812  // p_mu_D_hat
3813  if( act_set->q_D_hat() )
3814  calc_p_mu_D( *act_set, p_v(), p_z_hat(), NULL, &act_set->p_mu_D_hat() );
3815  // gamma_plus
3816  const size_type sa = act_set->s_map(ja);
3817  if(sa) {
3818  // This is not an initially fixed variable that returned to its
3819  // initial. The multiplier for this constriant may not be the
3820  // last element if an ADD/DROP was performed on the last iteration
3821  // in order to get here where the DROP was an initially fixed variable
3822  // that was freed and therefore the KKT system was augmented so this
3823  // multiplier is not the last element of z_hat(...).
3824  gamma_plus = z_hat_plus(sa);
3825  }
3826  else {
3827  // This must be an initially fixed variable that returned to its
3828  // initial bound. This will be the last element even if an ADD/DROP
3829  // was just performed since a drop would only remove elements from
3830  // p_mu_D_hat, not add them.
3831  gamma_plus = act_set->p_mu_D_hat()(act_set->q_D_hat());
3832  }
3833  // p_z_hat = p_z_hat / gamma_plus
3834  Vt_S( &p_z_hat(), 1.0 / gamma_plus );
3835  // p_v = p_v / gamma_plus
3836  Vt_S( &p_v(), 1.0 / gamma_plus );
3837  // Print gama_plus, p_z_hat, p_v and p_mu_D_hat
3838  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3839  *out
3840  << "\ngamma_plus = " << gamma_plus << std::endl;
3841  }
3842  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3843  *out
3844  << "\n||p_z_hat||inf = " << norm_inf(p_z_hat()) << std::endl;
3845  }
3846  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3847  *out
3848  << "\np_z_hat =\n" << p_z_hat();
3849  }
3850  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3851  *out
3852  << "\n||p_v||inf = " << norm_inf(p_v()) << std::endl;
3853  }
3854  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3855  *out
3856  << "\np_v =\n" << p_v();
3857  }
3858  if( act_set->q_D_hat() ) {
3859  // p_mu_D_hat = p_mu_D_hat / gamma_plus
3860  Vt_S( &act_set->p_mu_D_hat(), 1.0 / gamma_plus );
3861  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3862  *out
3863  << "\n||p_mu_D_hat||inf =\n" << norm_inf(act_set->p_mu_D_hat()) << std::endl;
3864  }
3865  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES ) {
3866  *out
3867  << "\np_mu_D_hat =\n" << act_set->p_mu_D_hat();
3868  }
3869  }
3870  }
3871  }
3872  case COMPUTE_STEP_LENGTHS: {
3873 
3874  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3875  *out
3876  << "\n*** COMPUTE_STEP_LENGTHS\n";
3877  }
3878  // Compute the dual infeasibility scaling
3879  const size_type q_hat = act_set->q_hat();
3880  dual_infeas_scale = 1.0;
3881 // if( q_hat )
3882 // dual_infeas_scale = my_max( dual_infeas_scale, norm_inf( act_set->z_hat() ) );
3883 // if( m )
3884 // dual_infeas_scale = my_max( dual_infeas_scale, norm_inf( (*v)(n_R+1,n_R+m) ) );
3885 // if( act_set->q_D_hat() )
3886 // dual_infeas_scale = my_max( dual_infeas_scale, norm_inf( act_set->mu_D_hat() ) );
3887 
3888  // Primal step length, t_P = beta * gamma_plus, z_plus = [ z_hat_plus; gama_plus ].
3889  // Or constraint ja is linearly dependent in which case p_x is zero so
3890  // t_P is infinite.
3891  t_P = beta * gamma_plus; // Could be < 0
3892  if( t_P < 0.0 ) {
3893  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3894  *out
3895  << "\nWarning, A near degenerate inequality constraint ja = " << ja
3896  << " is being added that has the wrong sign with:\n"
3897  << " t_P = " << t_P << std::endl
3898  << " dual_infeas_scale = " << dual_infeas_scale << std::endl
3899  << " norm_2_constr = " << norm_2_constr << std::endl
3900  << " |t_P/(norm_2_constr*dual_infeas_scale)| = "
3901  << std::fabs(t_P/(norm_2_constr*dual_infeas_scale))
3902  << " <= dual_infeas_tol = " << dual_infeas_tol() << std::endl;
3903  }
3904  summary_lines_counter = 0;
3905  if( !using_iter_refinement ) {
3906  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3907  *out << "We are not using iterative refinement yet so turn it on"
3908  << "\nthen recompute the steps ...\n";
3909  }
3910  using_iter_refinement = true;
3911  next_step = COMPUTE_SEARCH_DIRECTION;
3912  continue;
3913  }
3914  else {
3915  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
3916  *out
3917  << "We are already using iterative refinement so the QP algorithm is terminated!\n";
3918  }
3919  return DUAL_INFEASIBILITY;
3920  }
3921  }
3922  t_P = beta * gamma_plus; // Now guaranteed to be > 0
3923  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
3924  *out
3925  << "\nt_P = " << t_P << endl;
3926  }
3927 
3929  // Dual step length. Largest step t that does not cause violation in
3930  // dual feasibility (i.e. lagrange multipliers for inequalities are
3931  // dual feasible, or primal optimal ).
3932  // lambda_hat_new = lambda_hat + beta * t_D * p_lambda_hat must be dual feasible.
3933  t_D = inf;
3934  jd = 0;
3935  value_type max_feas_viol = 0.0; // Remember the amount of violation.
3936  int j_degen = 0; // remember which (if any) constraint was near
3937  // degenerate and had an incorrect sign.
3938  EBounds bnd_jd; // The bound of the constraint to be dropped.
3939 
3940  // Search through Lagrange multipliers in z_hat
3941  if( act_set->q_hat() ) {
3942  DVectorSlice z_hat = act_set->z_hat();
3943  DVectorSlice p_z_hat = act_set->p_z_hat();
3944  DVectorSlice::iterator
3945  z_itr = z_hat.begin(),
3946  p_z_itr = p_z_hat.begin();
3947  const size_type
3948  qq = assume_lin_dep_ja || (!assume_lin_dep_ja && return_to_init_fixed)
3949  ? q_hat : q_hat - 1;
3950  // Print header for s, j, z_hat(s), p_z_hat(s), bnds(s), t, t_D, jd
3951  if( qq > 0 && (int)output_level >= (int)OUTPUT_ACT_SET ) {
3952  *out
3953  << "\nComputing the maximum step for multiplers for dual feasibility\n\n"
3954  << right << setw(5) << "s"
3955  << right << setw(5) << "j"
3956  << right << setw(dbl_w) << "z_hat"
3957  << right << setw(dbl_w) << "p_z_hat"
3958  << right << setw(20) << "bnd"
3959  << right << setw(dbl_w) << "t"
3960  << right << setw(dbl_w) << "t_D"
3961  << right << setw(5) << "jd" << endl
3962  << right << setw(5) << "----"
3963  << right << setw(5) << "----"
3964  << right << setw(dbl_w) << "--------------"
3965  << right << setw(dbl_w) << "--------------"
3966  << right << setw(20) << "--------------"
3967  << right << setw(dbl_w) << "--------------"
3968  << right << setw(dbl_w) << "--------------"
3969  << right << setw(5) << "----" << endl;
3970  }
3971  for( int s = 1; s <= qq; ++s, ++z_itr, ++p_z_itr) {
3972  int j = act_set->ij_map(s);
3973  if( j > 0 ) {
3974  namespace ns = QPSchurPack;
3975  EBounds bnd = act_set->bnd(s);
3976  // Print first part of row for s, j, z_hat(s), p_z_hat(s), bnds(s) ....
3977  if( (int)output_level >= (int)OUTPUT_ACT_SET ) {
3978  *out
3979  << right << setw(5) << s
3980  << right << setw(5) << j
3981  << right << setw(dbl_w) << *z_itr
3982  << right << setw(dbl_w) << *p_z_itr
3983  << right << setw(20) << bnd_str(bnd);
3984  }
3985  value_type t = inf;
3986  // Lookout for degeneracy.
3987  bool j_is_degen = false;
3988  value_type viol;
3989  const int dual_feas_status
3990  = correct_dual_infeas(
3991  j,bnd,t_P,1.0,dual_infeas_tol(),DEGENERATE_MULT
3992  ,out,output_level,true,"z_hat(s)",&(*z_itr),&viol
3993  ,"p_z_hat(s)",&(*p_z_itr),"z_hat_plus(s)"
3994  , (assume_lin_dep_ja ? NULL: &z_hat_plus(s) ) );
3995  if( dual_feas_status < 0 ) {
3996  if( !using_iter_refinement ) {
3997  if( (int)output_level >= (int)OUTPUT_BASIC_INFO )
3998  *out << "We are not using iterative refinement yet so turn it on"
3999  << "\nthen recompute the steps ...\n";
4000  using_iter_refinement = true;
4001  next_step = COMPUTE_SEARCH_DIRECTION;
4002  continue;
4003  }
4004  else {
4005  if( (int)output_level >= (int)OUTPUT_BASIC_INFO )
4006  *out << "We are already using iterative refinement so the QP algorithm is terminated!\n";
4007  return DUAL_INFEASIBILITY;
4008  }
4009  }
4010  else if( dual_feas_status == 0 ) {
4011  j_is_degen = true;
4012  }
4013  // If we get here either the dual variable was feasible or it
4014  // was near degenerate and was corrected!
4015  const value_type feas_viol = beta*(*p_z_itr);
4016  if( bnd == EQUALITY )
4017  ; // We don't care
4018  else if( bnd == LOWER && feas_viol <= 0.0 )
4019  ; // dual feasible for all t > 0
4020  else if( bnd == UPPER && feas_viol >= 0.0 )
4021  ; // dual feasible for all t > 0
4022  else {
4023  // finite t.
4024  t = -beta*(*z_itr)/(*p_z_itr);
4025  if( t < t_D ) { // remember minimum step length
4026  t_D = t;
4027  jd = j;
4028  if(j_is_degen) j_degen = j;
4029  max_feas_viol = feas_viol;
4030  bnd_jd = bnd;
4031  }
4032  }
4033  // Print rest of row for ... t, t_D, jd
4034  if( (int)output_level >= (int)OUTPUT_ACT_SET ) {
4035  *out
4036  << right << setw(dbl_w) << t
4037  << right << setw(dbl_w) << t_D
4038  << right << setw(5) << jd << endl;
4039  }
4040  }
4041  }
4042  }
4043 
4044  // Search through Lagrange multipliers in mu_D_hat
4045  if( act_set->q_D_hat() ) {
4046  const QPSchurPack::QP::x_init_t &x_init = qp.x_init();
4047  const QPSchurPack::QP::i_x_X_map_t &i_x_X_map = qp.i_x_X_map();
4048  const size_type q_D_hat = act_set->q_D_hat();
4049  DVectorSlice mu_D_hat = act_set->mu_D_hat();
4050  DVectorSlice p_mu_D_hat = act_set->p_mu_D_hat();
4051  const size_type
4052  qD = assume_lin_dep_ja && return_to_init_fixed ? q_D_hat-1 : q_D_hat;
4053  // Print header for k, i, mu_D_hat(k), p_mu_D_hat(k), x_init(k), t, t_D, jd
4054  if( qD > 0 && (int)output_level >= (int)OUTPUT_ACT_SET ) {
4055  *out
4056  << "\nComputing the maximum step for multiplers for dual feasibility\n\n"
4057  << right << setw(5) << "k"
4058  << right << setw(5) << "i"
4059  << right << setw(dbl_w) << "mu_D_hat"
4060  << right << setw(dbl_w) << "p_mu_D_hat"
4061  << right << setw(20) << "x_init"
4062  << right << setw(dbl_w) << "t"
4063  << right << setw(dbl_w) << "t_D"
4064  << right << setw(5) << "jd" << endl
4065  << right << setw(5) << "----"
4066  << right << setw(5) << "----"
4067  << right << setw(dbl_w) << "--------------"
4068  << right << setw(dbl_w) << "--------------"
4069  << right << setw(20) << "--------------"
4070  << right << setw(dbl_w) << "--------------"
4071  << right << setw(dbl_w) << "--------------"
4072  << right << setw(5) << "----" << endl;
4073  }
4074  GenPermMatrixSlice::const_iterator
4075  Q_XD_itr = act_set->Q_XD_hat().begin(),
4076  Q_XD_end = Q_XD_itr + qD;
4077  for( ; Q_XD_itr != Q_XD_end; ++Q_XD_itr ) {
4078  const size_type k = Q_XD_itr->col_j();
4079  const size_type i = Q_XD_itr->row_i();
4080  DVectorSlice::iterator
4081  mu_D_itr = mu_D_hat.begin() + (k-1),
4082  p_mu_D_itr = p_mu_D_hat.begin() + (k-1);
4083  const size_type l = act_set->l_fxfx(k);
4084  EBounds bnd = qp.x_init()(i);
4085  // Print first part of row for s, j, z_hat(s), p_z_hat(s), bnds(s) ....
4086  if( (int)output_level >= (int)OUTPUT_ACT_SET ) {
4087  *out
4088  << right << setw(5) << k
4089  << right << setw(5) << i
4090  << right << setw(dbl_w) << *mu_D_itr
4091  << right << setw(dbl_w) << *p_mu_D_itr
4092  << right << setw(20) << bnd_str(bnd);
4093  }
4094  value_type t = inf;
4095  // Lookout for degeneracy.
4096  bool j_is_degen = false;
4097  value_type viol;
4098  const int dual_feas_status
4099  = correct_dual_infeas(
4100  i,bnd,t_P,1.0,dual_infeas_tol(),DEGENERATE_MULT
4101  ,out,output_level,true,"mu_D_hat(k)",&(*mu_D_itr),&viol
4102  ,"p_mu_D_hat(k)",&(*p_mu_D_itr) );
4103  if( dual_feas_status < 0 ) {
4104  if( !using_iter_refinement ) {
4105  if( (int)output_level >= (int)OUTPUT_BASIC_INFO )
4106  *out << "We are not using iterative refinement yet so turn it on"
4107  << "\nthen recompute the steps ...\n";
4108  using_iter_refinement = true;
4109  next_step = COMPUTE_SEARCH_DIRECTION;
4110  continue;
4111  }
4112  else {
4113  if( (int)output_level >= (int)OUTPUT_BASIC_INFO )
4114  *out << "We are already using iterative refinement so the QP algorithm is terminated!\n";
4115  return DUAL_INFEASIBILITY;
4116  }
4117  }
4118  else if( dual_feas_status == 0 ) {
4119  j_is_degen = true;
4120  }
4121  // If we get here either the dual variable was feasible or it
4122  // was near degenerate and was corrected!
4123  const value_type feas_viol = beta*(*p_mu_D_itr);
4124  if( bnd == EQUALITY )
4125  ; // We don't care
4126  else if( bnd == LOWER && feas_viol <= 0.0 )
4127  ; // dual feasible for all t > 0
4128  else if( bnd == UPPER && feas_viol >= 0.0 )
4129  ; // dual feasible for all t > 0
4130  else {
4131  // finite t.
4132  t = -beta*(*mu_D_itr)/(*p_mu_D_itr);
4133  if( t < t_D ) { // remember minimum step length
4134  t_D = t;
4135  jd = -i;
4136  if(j_is_degen) j_degen = jd;
4137  max_feas_viol = feas_viol;
4138  bnd_jd = bnd;
4139  }
4140  }
4141  // Print rest of row for ... t, t_D, jd
4142  if( (int)output_level >= (int)OUTPUT_ACT_SET ) {
4143  *out
4144  << right << setw(dbl_w) << t
4145  << right << setw(dbl_w) << t_D
4146  << right << setw(5) << jd << endl;
4147  }
4148  }
4149  }
4150  // Print t_D, jd
4151  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
4152  *out
4153  << "\nt_D = " << t_D << endl
4154  << "jd = " << jd << endl;
4155  }
4156  if( jd == j_degen && jd != 0 && t_D < t_P ) {
4157  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4158  *out
4159  << "\nWarning, the near degenerate constraint j = "
4160  << jd << " which had the incorrect sign\nand was adjusted "
4161  << "was selected to be dropped from the active set.\n";
4162  }
4163  }
4164  // Print end of row for rank if the right print level
4165  if( assume_lin_dep_ja && !schur_comp_update_failed && (int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
4166  if( t_P < huge_primal_step() )
4167  *out << right << setw(6) << "LI" << endl;
4168  else
4169  *out << right << setw(6) << "LD" << endl;
4170  out->flush();
4171  --summary_lines_counter;
4172  }
4173  // Print start of row for itr, q_hat, q(+), q_D, q_C, q_F, change, type, index, bound, violation
4174  if( t_D < t_P && (int)output_level == (int)OUTPUT_ITER_SUMMARY ) {
4175  *out
4176  << right << setw(6) << itr // itr
4177  << right << setw(6) << act_set->q_hat() // q_hat
4178  << right << setw(6) << act_set->q_plus_hat() // q(+)
4179  << right << setw(6) << act_set->q_F_hat() // q_F
4180  << right << setw(6) << act_set->q_C_hat() // q_C
4181  << right << setw(6) << act_set->q_D_hat() // q_D
4182  << right << setw(8) << "DROP" // change
4183  << right << setw(9); // type
4184  if( jd < 0 ) {
4185  *out << "X_F";
4186  }
4187  else if( jd <= n ) {
4188  if( bnd_jd == qp.x_init()(jd) )
4189  *out << "X_F_C_F";
4190  else
4191  *out << "R_X_R";
4192  }
4193  else {
4194  *out << "GEN";
4195  }
4196  *out
4197  << right << setw(6) << jd // index
4198  << right << setw(10) << bnd_str(bnd_jd) // bound
4199  << right << setw(dbl_w) << max_feas_viol // violation
4200  << right << setw(6) << "LI" << endl; // rank (this should be true!)
4201  }
4202  }
4203  case TAKE_STEP: {
4204  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
4205  *out
4206  << "\n*** TAKE_STEP\n";
4207  }
4208  if( t_P >= huge_primal_step() && t_D >= huge_dual_step() ) {
4209  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4210  *out
4211  << "Error, QP is infeasible, inconsistent constraint a("<<ja<<") detected\n";
4212  }
4213  if( using_iter_refinement ) {
4214  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4215  *out
4216  << "We are already using iterative refinement so the QP algorithm is terminated!\n";
4217  }
4218  return INFEASIBLE_CONSTRAINTS;
4219  }
4220  else {
4221  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4222  *out << "We are not using iterative refinement yet so turn it on";
4223  if(assume_lin_dep_ja)
4224  *out << "\nthen pick another violated constriant to add ... \n";
4225  else
4226  *out << "\nthen recompute the steps ...\n";
4227  }
4228  summary_lines_counter = 0;
4229  last_jd = 0; // erase this memory!
4230  last_ja = 0; // ..
4231  using_iter_refinement = true;
4232  if(assume_lin_dep_ja) {
4233  EIterRefineReturn status = iter_refine(
4234  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, act_set->q_hat() ? &act_set->d_hat() : NULL
4235  ,v, act_set->q_hat() ? &act_set->z_hat() : NULL
4236  ,iter_refine_num_resid, iter_refine_num_solves
4237  );
4238  next_step = PICK_VIOLATED_CONSTRAINT;
4239  }
4240  else {
4241  // Iterative refinement will be performed there
4242  next_step = COMPUTE_SEARCH_DIRECTION;
4243  }
4244  continue;
4245  }
4246  }
4247  else if( t_P > t_D ) {
4248  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
4249  if( t_P >= huge_primal_step() ) {
4250  *out
4251  << "\n*** (b) Dual Step (t_P = " << t_P << " >= huge_primal_step = "
4252  << huge_primal_step() << ")\n";
4253  }
4254  else {
4255  *out
4256  << "\n*** (b) Partial Primal-Dual Step\n";
4257  }
4258  }
4259  // Check for cycling
4260  if( ja == last_jd && jd == last_ja ) {
4261  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4262  *out
4263  << "\n\nQPSchur::qp_algo(...) : Error, the constraint "
4264  << "a(" << ja << ") with violation\n"
4265  << "(a(ja)'*x - b_a) = (" << con_ja_val
4266  << " - " << b_a << ") = " << (con_ja_val - b_a) << "\n"
4267  << "we are adding to the active set and the constraint constriant\n"
4268  << "a(" << jd << ") we are dropping were just dropped and added respectively!\n"
4269  << "The algorithm is cycling!\n";
4270  }
4271  if( using_iter_refinement ) {
4272  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4273  *out
4274  << "We are already using iterative refinement so the QP algorithm is terminated!\n";
4275  }
4276  return SUBOPTIMAL_POINT;
4277  }
4278  else {
4279  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4280  *out << "We are not using iterative refinement yet so turn it on";
4281  if(assume_lin_dep_ja)
4282  *out << "\nthen pick another violated constriant to add ... \n";
4283  else
4284  *out << "\nthen recompute the steps ...\n";
4285  }
4286  summary_lines_counter = 0;
4287  last_jd = 0; // erase this memory!
4288  last_ja = 0; // ..
4289  using_iter_refinement = true;
4290  if(assume_lin_dep_ja) {
4291  EIterRefineReturn status = iter_refine(
4292  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, act_set->q_hat() ? &act_set->d_hat() : NULL
4293  ,v, act_set->q_hat() ? &act_set->z_hat() : NULL
4294  ,iter_refine_num_resid, iter_refine_num_solves
4295  );
4296  next_step = PICK_VIOLATED_CONSTRAINT;
4297  }
4298  else {
4299  // Iterative refinement will be performed there
4300  next_step = COMPUTE_SEARCH_DIRECTION;
4301  }
4302  continue;
4303  }
4304  }
4305  // Update the augmented KKT system
4306  try {
4307  if( assume_lin_dep_ja ) {
4308  if(act_set->drop_add_constraints( jd, ja, bnd_ja, true, out, output_level ))
4309  summary_lines_counter = 0;
4310  }
4311  else {
4312  if(act_set->drop_constraint( jd, out, output_level, true, true ))
4313  summary_lines_counter = 0;
4314  }
4315  }
4316  catch( const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
4317  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4318  *out
4319  << "\n\nSchur complement appears to be singular and should not be:\n"
4320  << excpt.what()
4321  << "\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4322  }
4323  return NONCONVEX_QP;
4324  }
4325  catch( const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
4326  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4327  *out
4328  << "\n\nSchur complement appears to have the wrong inertia:\n"
4329  << excpt.what()
4330  << "\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4331  }
4332  return NONCONVEX_QP;
4333  }
4334  // z_hat = z_hat + beta * t_D * p_z_hat
4335  if(act_set->q_hat())
4336  Vp_StV( &act_set->z_hat(), beta * t_D, act_set->p_z_hat() );
4337  // v = v + beta * t_D * p_v
4338  Vp_StV( v, beta * t_D, p_v() );
4339  // mu_D_hat = mu_D_hat + beta * t_D * p_mu_D_hat
4340  if(act_set->q_D_hat())
4341  Vp_StV( &act_set->mu_D_hat(), beta * t_D, act_set->p_mu_D_hat() );
4342 
4343  ++(*num_drops);
4344 
4345  if( (int)output_level >= (int)OUTPUT_ITER_STEPS )
4346  {
4347  *out
4348  << "\nUpdated primal and dual variables:\n"
4349  << "\n||v||inf = " << norm_inf(*v) << endl;
4350  if(act_set->q_hat()) {
4351  *out
4352  << "||z_hat||inf = " << norm_inf(act_set->z_hat()) << endl;
4353  }
4354  if(act_set->q_D_hat()) {
4355  *out
4356  << "max(|mu_D_hat(i)|) = " << norm_inf(act_set->mu_D_hat()) << endl
4357  << "min(|mu_D_hat(i)|) = " << min_abs(act_set->mu_D_hat()) << endl;
4358  }
4359  }
4360  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES )
4361  {
4362  *out << "\nv = \n" << *v << endl;
4363  if(assume_lin_dep_ja) {
4364  *out
4365  << "\nPrinting active set after dropping constraint jd = " << jd
4366  << " and adding constraint ja = " << ja << " ...\n";
4367  }
4368  else {
4369  *out
4370  << "\nPrinting active set after dropping constraint jd = " << jd
4371  << " ...\n";
4372  }
4373  dump_act_set_quantities( *act_set, *out );
4374  }
4375  last_jd = jd;
4376  assume_lin_dep_ja = false; // If we get here then we know these are true!
4377  next_step = COMPUTE_SEARCH_DIRECTION;
4378  continue;
4379  }
4380  else { // t_P < t_D
4381  if( (int)output_level >= (int)OUTPUT_ITER_STEPS ) {
4382  *out
4383  << "\n*** (c) Full Primal-Dual Step\n";
4384  }
4385  ++(*num_adds);
4386  if( !assume_lin_dep_ja ) {
4387  act_set->z_hat() = z_hat_plus;
4388  *v = v_plus;
4389  }
4390  else {
4391  bool threw_exception = false;
4392  try {
4393  if(act_set->add_constraint( ja, bnd_ja, true, out, output_level, true, true ))
4394  summary_lines_counter = 0;
4395  }
4396  catch( const MatrixSymAddDelUpdateable::SingularUpdateException& excpt ) {
4397  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4398  *out
4399  << "\n\nSchur complement appears to be singular and should not be:\n"
4400  << excpt.what() << std::endl;
4401  }
4402  threw_exception = true;
4403  }
4404  catch( const MatrixSymAddDelUpdateable::WrongInertiaUpdateException& excpt ) {
4405  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4406  *out
4407  << "\n\nSchur complement appears to have the wrong inertia:\n"
4408  << excpt.what() << std::endl;
4409  }
4410  threw_exception = true;
4411  }
4412  if( threw_exception ) {
4413  if( !using_iter_refinement ) {
4414  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4415  *out << "We are not using iterative refinement yet so turn it on and\n"
4416  << "go back and pick a new violated constraint to add to the active set ...\n";
4417  }
4418  using_iter_refinement = true;
4419  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4420  *out
4421  << "\n\nPerforming iterative refinement on v, z_hat system ...\n";
4422  }
4423  summary_lines_counter = 0;
4424  // [ Ko U_hat ] [ v ] = [ fo ]
4425  // [ U_hat' V_hat ] [ z_hat ] [ d_hat ]
4426  EIterRefineReturn status = iter_refine(
4427  *act_set, out, output_level, -1.0, &qp.fo(), -1.0, &act_set->d_hat()
4428  ,v, &act_set->z_hat()
4429  ,iter_refine_num_resid, iter_refine_num_solves
4430  );
4431  next_step = PICK_VIOLATED_CONSTRAINT;
4432  continue;
4433  }
4434  else {
4435  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4436  *out << "Darn, we are already using iterative refinement!"
4437  << "\nThe QP appears to be nonconvex and we therefore terminate the primal-dual QP algorithm!\n";
4438  }
4439  return NONCONVEX_QP;
4440  }
4441  }
4442  // z_hat = z_hat + beta * t_P * p_z_hat
4443  if(act_set->q_hat())
4444  Vp_StV( &act_set->z_hat(), beta * t_P, act_set->p_z_hat() );
4445  // v = v + beta * t_P * p_v
4446  Vp_StV( v, beta * t_P, p_v() );
4447  }
4448  // mu_D_hat = mu_D_hat + beta * t_P * p_mu_D_hat
4449  if(act_set->q_D_hat())
4450  Vp_StV( &act_set->mu_D_hat(), beta * t_P, act_set->p_mu_D_hat() );
4451 
4452 
4453  if( (int)output_level >= (int)OUTPUT_ITER_STEPS )
4454  {
4455  *out
4456  << "\nUpdated primal and dual variables:\n"
4457  << "\n||v||inf = " << norm_inf(*v) << endl;
4458  if(act_set->q_hat()) {
4459  *out
4460  << "||z_hat||inf = " << norm_inf(act_set->z_hat()) << endl;
4461  }
4462  if(act_set->q_D_hat()) {
4463  *out
4464  << "max(|mu_D_hat(i)|) = " << norm_inf(act_set->mu_D_hat()) << endl
4465  << "min(|mu_D_hat(i)|) = " << min_abs(act_set->mu_D_hat()) << endl;
4466  }
4467  }
4468  if( (int)output_level >= (int)OUTPUT_ITER_QUANTITIES )
4469  {
4470  *out << "\nv = \n" << *v << endl;
4471  if( assume_lin_dep_ja ) {
4472  *out
4473  << "\nPrinting active set after adding constraint ja = " << ja
4474  << " ...\n";
4475  dump_act_set_quantities( *act_set, *out );
4476  }
4477  else {
4478  if(act_set->q_hat())
4479  *out << "\nz_hat =\n" << act_set->z_hat();
4480  if(act_set->q_D_hat())
4481  *out << "\nmu_D_hat =\n" << act_set->mu_D_hat();
4482  }
4483  }
4484  assume_lin_dep_ja = false; // If we get here then we know these are true!
4485  next_step = PICK_VIOLATED_CONSTRAINT;
4486  continue;
4487  }
4488  }
4489  default:
4490  TEUCHOS_TEST_FOR_EXCEPT(true); // only a local programming error
4491  }
4492  }
4493 
4494  } // end try
4495  catch( std::exception& excpt ) {
4496  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4497  *out
4498  << "\n\n*** Caught a standard exception :\n"
4499  << excpt.what()
4500  << "\n*** Rethrowing the exception ...\n";
4501  }
4502  throw;
4503  }
4504  catch(...) {
4505  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4506  *out
4507  << "\n\n*** Caught an unknown exception. Rethrowing the exception ...\n";
4508  }
4509  throw;
4510  }
4511  // If you get here then the maximum number of QP iterations has been exceeded
4512  return MAX_ITER_EXCEEDED;
4513 }
4514 
4515 void QPSchur::set_x( const ActiveSet& act_set, const DVectorSlice& v, DVectorSlice* x )
4516 {
4517  using BLAS_Cpp::no_trans;
4518  using LinAlgOpPack::V_MtV;
4519  using LinAlgOpPack::Vp_MtV;
4520 
4521  // x = Q_R * v(1:n_R) + Q_X * b_X + P_XF_hat * z_hat
4522  V_MtV( x, act_set.qp().Q_R(), no_trans, v(1,act_set.qp().n_R()) );
4523  if( act_set.qp().n() > act_set.qp().n_R() )
4524  Vp_MtV( x, act_set.qp().Q_X(), no_trans, act_set.qp().b_X() );
4525  if( act_set.q_F_hat() )
4526  Vp_MtV( x, act_set.P_XF_hat(), no_trans, act_set.z_hat() );
4527 }
4528 
4530  const ActiveSet& act_set, const DVectorSlice& v
4531  ,SpVector* mu, DVectorSlice* lambda, SpVector* lambda_breve
4532  )
4533 {
4534  using BLAS_Cpp::no_trans;
4535  using LinAlgOpPack::V_MtV;
4536  using LinAlgOpPack::Vp_MtV;
4537  using LinAlgOpPack::Vp_StMtV;
4538  using LinAlgOpPack::V_MtV;
4539  using LinAlgOpPack::Vp_MtV;
4540  using AbstractLinAlgPack::V_MtV;
4541  using AbstractLinAlgPack::Vp_MtV;
4542  namespace GPMSTP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
4543  using Teuchos::Workspace;
4545 
4546  const size_type
4547  n = act_set.qp().n(),
4548  n_R = act_set.qp().n_R(),
4549  m = act_set.qp().m(),
4550  m_breve = act_set.qp().constraints().m_breve(),
4551  q_hat = act_set.q_hat();
4553  &x_init = act_set.qp().x_init();
4554 
4555  //
4556  // mu = P_plus_hat(1:n,:) * z_hat + Q_XD_hat * mu_D + (steps for initially fixed
4557  // variables fixed to the other bounds)
4558  //
4559  // lambda_breve = P_plus_hat(n+1:n+m_breve,:) * z_hat
4560  //
4561  typedef SpVector::element_type ele_t;
4562  mu->resize( n, n-m ); // Resize for the maxinum number
4563  lambda_breve->resize( m_breve, n-m ); // of active constraints possible.
4564  // mu += Q_XD_hat * mu_D_hat
4565  if( act_set.q_D_hat() )
4566  Vp_MtV( mu, act_set.Q_XD_hat(), no_trans, act_set.mu_D_hat() );
4567  // Set all the multipliers in z_hat
4568  if(q_hat){
4569  const DVectorSlice
4570  z_hat = act_set.z_hat();
4571  for( size_type s = 1; s <= q_hat; ++s ) {
4572  const int ij = act_set.ij_map(s);
4573  if(ij > 0) {
4574  const size_type j = ij;
4575  if( j <= n )
4576  mu->add_element(ele_t(j,z_hat(s)));
4577  else
4578  lambda_breve->add_element(ele_t(j-n,z_hat(s)));
4579  }
4580  }
4581  }
4582  mu->sort();
4583  lambda_breve->sort();
4584  // lambda = v(n_R+1,n_R+m)
4585  if( m ) {
4586  *lambda = v(n_R+1,n_R+m);
4587  }
4588 }
4589 
4590 bool QPSchur::timeout_return( StopWatchPack::stopwatch* timer, std::ostream *out, EOutputLevel output_level ) const
4591 {
4592  const value_type minutes = timer->read() / 60;
4593  if( minutes >= max_real_runtime() ) {
4594  if( (int)output_level >= (int)OUTPUT_BASIC_INFO ) {
4595  *out
4596  << "\n*** Runtime = " << minutes << " min >= max_real_runtime = " << max_real_runtime() << " min!\n"
4597  << "Must terminite the algorithm!\n";
4598  }
4599  return true;
4600  }
4601  return false;
4602 }
4603 
4606  const ActiveSet &act_set
4607  ,std::ostream *out
4608  ,EOutputLevel output_level
4609  ,const value_type ao
4610  ,const DVectorSlice *bo
4611  ,const value_type aa
4612  ,const DVectorSlice *ba
4613  ,DVectorSlice *v
4614  ,DVectorSlice *z
4615  ,size_type *iter_refine_num_resid
4616  ,size_type *iter_refine_num_solves
4617  )
4618 {
4619  using std::endl;
4620  using std::setw;
4621  using std::left;
4622  using std::right;
4623  using BLAS_Cpp::no_trans;
4624  using BLAS_Cpp::trans;
4625  using DenseLinAlgPack::norm_inf;
4627  using LinAlgOpPack::Vp_V;
4628  using LinAlgOpPack::Vp_StMtV;
4629  using LinAlgOpPack::V_InvMtV;
4630  using Teuchos::Workspace;
4632 
4633  typedef DenseLinAlgPack::value_type extra_value_type;
4634 
4635  const value_type small_num = std::numeric_limits<value_type>::min();
4636 
4637  const int int_w = 8;
4638  const char int_ul[] = "------";
4639  const int dbl_min_w = 20;
4640  const int dbl_w = ( out ? my_max(dbl_min_w,int(out->precision()+8)): 20 );
4641  const char dbl_ul[] = "------------------";
4642 
4643  const QPSchurPack::QP
4644  &qp = act_set.qp();
4645  const MatrixSymOpNonsing
4646  &Ko = qp.Ko(),
4647  &S_hat = act_set.S_hat();
4648  const MatrixOp
4649  &U_hat = act_set.U_hat();
4651  n = qp.n(),
4652  n_R = qp.n_R(),
4653  m = qp.m(),
4654  q_hat = act_set.q_hat();
4655  const DVectorSlice
4656  fo = qp.fo(),
4657  d_hat = (q_hat ? act_set.d_hat() : DVectorSlice());
4658 
4659  Workspace<extra_value_type>
4660  ext_ro_ws(wss,n_R+m),
4661  ext_ra_ws(wss,q_hat);
4663  ext_ro(&ext_ro_ws[0],ext_ro_ws.size()),
4664  ext_ra(ext_ra_ws.size()?&ext_ra_ws[0]:NULL,ext_ra_ws.size());
4665  Workspace<value_type>
4666  ro_ws(wss,n_R+m),
4667  ra_ws(wss,q_hat),
4668  t1_ws(wss,n_R+m),
4669  del_v_ws(wss,n_R+m),
4670  del_z_ws(wss,q_hat),
4671  v_itr_ws(wss,n_R+m),
4672  z_itr_ws(wss,q_hat);
4673  DVectorSlice
4674  ro(&ro_ws[0],ro_ws.size()),
4675  ra(ra_ws.size()?&ra_ws[0]:NULL,ra_ws.size()),
4676  t1(&t1_ws[0],t1_ws.size()),
4677  del_v(&del_v_ws[0],del_v_ws.size()),
4678  del_z(del_z_ws.size()?&del_z_ws[0]:NULL,del_z_ws.size()),
4679  v_itr(&v_itr_ws[0],v_itr_ws.size()),
4680  z_itr(z_itr_ws.size()?&z_itr_ws[0]:NULL,z_itr_ws.size());
4681 
4682  // Accumulate into temporary variables
4683  v_itr = *v;
4684  if(q_hat)
4685  z_itr = *z;
4686 
4687  // Print summary header
4688  if( out && (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4689  *out
4690  << "\nBeginning iterative refinement ..."
4691  << "\niter_refine_opt_tol = " << iter_refine_opt_tol()
4692  << ", iter_refine_feas_tol = " << iter_refine_feas_tol()
4693  << "\niter_refine_min_iter = " << iter_refine_min_iter()
4694  << ", iter_refine_max_iter = " << iter_refine_max_iter() << "\n\n";
4695  //
4696  *out
4697  << right << setw(int_w) << "ir_itr"
4698  << right << setw(dbl_w) << "roR_scaling"
4699  << right << setw(dbl_w) << "||roR||s"
4700  << left << setw(1) << " ";
4701  if(m) {
4702  *out
4703  << right << setw(dbl_w) << "rom_scaling"
4704  << right << setw(dbl_w) << "||rom||s"
4705  << left << setw(1) << " ";
4706  }
4707  if(q_hat) {
4708  *out
4709  << right << setw(dbl_w) << "ra_scaling"
4710  << right << setw(dbl_w) << "||ra||s"
4711  << left << setw(1) << " ";
4712  }
4713  *out
4714  << right << setw(dbl_w) << "||del_v||/||v||inf"
4715  << right << setw(dbl_w) << "||del_z||/||z||inf"
4716  << endl;
4717  //
4718  *out
4719  << right << setw(int_w) << int_ul
4720  << right << setw(dbl_w) << dbl_ul
4721  << right << setw(dbl_w) << dbl_ul
4722  << left << setw(1) << " ";
4723  if(m) {
4724  *out
4725  << right << setw(dbl_w) << dbl_ul
4726  << right << setw(dbl_w) << dbl_ul
4727  << left << setw(1) << " ";
4728  }
4729  if(q_hat) {
4730  *out
4731  << right << setw(dbl_w) << dbl_ul
4732  << right << setw(dbl_w) << dbl_ul
4733  << left << setw(1) << " ";
4734  }
4735  *out
4736  << right << setw(dbl_w) << dbl_ul
4737  << right << setw(dbl_w) << dbl_ul
4738  << endl;
4739  }
4740  //
4741  // Perform iterative refinement iterations
4742  //
4743  EIterRefineReturn return_status = ITER_REFINE_NOT_PERFORMED;
4744  value_type
4745  roR_nrm_o, rom_nrm_o, ra_nrm_o,
4746  roR_nrm, rom_nrm, ra_nrm;
4747  for( size_type iter_refine_k = 0;
4748  iter_refine_k < iter_refine_max_iter() && return_status != ITER_REFINE_CONVERGED;
4749  ++iter_refine_k)
4750  {
4751  //
4752  // Compute the residual (in extended precision?)
4753  //
4754  // [ ro ] = [ Ko U_hat ] [ v ] + [ ao*bo ]
4755  // [ ra ] [ U_hat' V_hat ] [ z ] [ aa*ba ]
4756  //
4757  value_type
4758  roR_scaling = 0.0,
4759  rom_scaling = 0.0,
4760  ra_scaling = 0.0;
4761  ++(*iter_refine_num_resid);
4762  calc_resid(
4763  act_set
4764  ,v_itr, z_itr
4765  ,ao
4766  ,bo
4767  ,&ext_ro
4768  ,&roR_scaling
4769  ,m ? &rom_scaling : NULL
4770  ,aa
4771  ,ba
4772  ,q_hat ? &ext_ra : NULL
4773  ,q_hat ? &ra_scaling : NULL
4774  );
4775  std::copy(ext_ro.begin(),ext_ro.end(),ro.begin()); // Convert back to standard precision
4776  if(q_hat) std::copy(ext_ra.begin(),ext_ra.end(),ra.begin());
4777  //
4778  // Calcuate convergence criteria
4779  //
4780  roR_nrm = norm_inf(ro(1,n_R));
4781  rom_nrm = (m ? norm_inf(ro(n_R+1,n_R+m)) : 0.0);
4782  ra_nrm = (q_hat ? norm_inf(ra) : 0.0);
4783  if( iter_refine_k == 0 ) {
4784  roR_nrm_o = roR_nrm;
4785  rom_nrm_o = rom_nrm;
4786  ra_nrm_o = rom_nrm;
4787  }
4788  const bool
4789  is_roR_conv = roR_nrm / (1.0 + roR_scaling) < iter_refine_opt_tol(),
4790  is_rom_conv = (m ? rom_nrm / (1.0 + rom_scaling) < iter_refine_feas_tol() : true ),
4791  is_ra_conv = (q_hat ? ra_nrm / (1.0 + ra_scaling) < iter_refine_feas_tol() : true ),
4792  is_conv = is_roR_conv && is_rom_conv && is_ra_conv;
4793  //
4794  // Print beginning of summary line for residuals
4795  //
4796  if( out && (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4797  *out
4798  << right << setw(int_w) << iter_refine_k
4799  << right << setw(dbl_w) << roR_scaling
4800  << right << setw(dbl_w) << (roR_nrm / (1.0 + roR_scaling))
4801  << left << setw(1) << (is_roR_conv ? "*" : " ");
4802  if(m) {
4803  *out
4804  << right << setw(dbl_w) << rom_scaling
4805  << right << setw(dbl_w) << (rom_nrm /(1.0 + rom_scaling))
4806  << left << setw(1) << (is_rom_conv ? "*" : " ");
4807  }
4808  if(q_hat) {
4809  *out
4810  << right << setw(dbl_w) << ra_scaling
4811  << right << setw(dbl_w) << (ra_nrm /(1.0 + ra_scaling))
4812  << left << setw(1) << (is_ra_conv ? "*" : " ");
4813  }
4814  }
4815  //
4816  // Check for convergence
4817  //
4818  if( iter_refine_k + 1 < iter_refine_min_iter() ) {
4819  // Keep on going even if we have converged to desired tolerances!
4820  }
4821  else if( is_conv ) {
4822  if( out && (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4823  *out
4824  << right << setw(dbl_w) << "-"
4825  << right << setw(dbl_w) << "-"
4826  << endl;
4827  }
4828  if( iter_refine_k == 0 )
4829  return_status = ITER_REFINE_NOT_NEEDED;
4830  else
4831  return_status = ITER_REFINE_CONVERGED;
4832  break;
4833  }
4834  // Make sure we have made progress
4835  if( roR_nrm_o < roR_nrm && rom_nrm_o < rom_nrm && ra_nrm_o < rom_nrm ) {
4836  return_status = ITER_REFINE_NOT_IMPROVED;
4837  break; // No progress was make in converging the equations!
4838  }
4839  //
4840  // Solve for the steps
4841  //
4842  // [ Ko U_hat ] [ del_v ] = [ ro ]
4843  // [ U_hat' V_hat ] [ del_z ] [ ra ]
4844  //
4845  ++(*iter_refine_num_solves);
4846  if( q_hat ) {
4847  // del_z = inv(S_hat)*(ra - U_hat'*inv(Ko)*ro)
4848  V_InvMtV( &t1, Ko, no_trans, ro );
4849  calc_z( act_set.S_hat(), ra, act_set.U_hat(), &t1, &del_z );
4850  }
4851  calc_v( Ko, &ro, U_hat, del_z, &del_v );
4852  //
4853  // Compute steps:
4854  //
4855  // v += -del_v
4856  // z += -del_z
4857  //
4858  Vp_StV( &v_itr, -1.0, del_v );
4859  if( q_hat )
4860  Vp_StV( &z_itr, -1.0, del_z );
4861  //
4862  // Print rest of summary line for steps
4863  //
4864  if( out && (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4865  *out
4866  << right << setw(dbl_w) << norm_inf(del_v) / (norm_inf(v_itr) + small_num)
4867  << right << setw(dbl_w) << norm_inf(del_z) / (norm_inf(z_itr) + small_num)
4868  << endl;
4869  }
4870  }
4871  if( iter_refine_max_iter() == 0 ) {
4872  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4873  *out
4874  << "\nWarning, iter_refine_max_iter == 0. Iterative refinement was not performed."
4875  << "\nLeaving the original solution intact ...\n";
4876  }
4877  return_status = ITER_REFINE_NOT_PERFORMED;
4878  }
4879  else {
4880  if( iter_refine_max_iter() == 1 ) {
4881  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4882  *out
4883  << "\nWarning, iter_refine_max_iter == 1. Only one step of iterative refinement"
4884  << "was performed and the step is taken which out checking the residual ...\n";
4885  }
4886  *v = v_itr;
4887  if(q_hat)
4888  *z = z_itr;
4889  return_status = ITER_REFINE_ONE_STEP;
4890  }
4891  else if( roR_nrm_o < roR_nrm && rom_nrm_o < rom_nrm && ra_nrm_o < rom_nrm ) {
4892  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4893  *out
4894  << "\nNo progress was made in reducing the residuals."
4895  << "\nLeaving the original solution intact ...\n";
4896  }
4897  return_status = ITER_REFINE_NOT_IMPROVED;
4898  }
4899  else {
4900  // The residuals were at least not increased so let's take the new solution
4901  *v = v_itr;
4902  if(q_hat)
4903  *z = z_itr;
4904  if( return_status != ITER_REFINE_CONVERGED && return_status != ITER_REFINE_NOT_NEEDED ) {
4905  if( (int)output_level >= (int)OUTPUT_ITER_SUMMARY ) {
4906  *out
4907  << "\nThe residuals were not converged but they were not increased either."
4908  << "\nTake the new point anyway ...\n";
4909  }
4910  return_status = ITER_REFINE_IMPROVED;
4911  }
4912  }
4913  }
4914  return return_status;
4915 }
4916 
4917 // private static member functions for QPSchur
4918 
4920  const ActiveSet& act_set, std::ostream& out
4921  ,bool print_S_hat
4922  )
4923 {
4924  using std::endl;
4925  using std::setw;
4926  using std::left;
4927  using std::right;
4928 
4929  const QPSchurPack::QP
4930  &qp = act_set.qp();
4932  &constraints = qp.constraints();
4933 
4934  const int int_w = 10;
4935  const char int_ul[] = "--------";
4936  const int dbl_min_w = 20;
4937  const int dbl_w = my_max(dbl_min_w,int(out.precision()+8));
4938  const char dbl_ul[] = "------------------";
4939 
4940  out << "\n*** Dumping the current active set ***\n"
4941  << "\nDimensions of the current active set:\n"
4942  << "\nn = " << right << setw(int_w) << qp.n() << " (Number of unknowns)"
4943  << "\nn_R = " << right << setw(int_w) << qp.n_R() << " (Number of initially free variables in Ko)"
4944  << "\nm = " << right << setw(int_w) << qp.m() << " (Number of initially fixed variables not in Ko)"
4945  << "\nm_breve = " << right << setw(int_w) << constraints.m_breve() << " (Number of extra general equality/inequality constriants)"
4946  << "\nq_hat = " << right << setw(int_w) << act_set.q_hat() << " (Number of augmentations to the initial KKT system Ko)"
4947  << "\nq_plus_hat = " << right << setw(int_w) << act_set.q_plus_hat() << " (Number of added variable bounds and general constraints)"
4948  << "\nq_F_hat = " << right << setw(int_w) << act_set.q_F_hat() << " (Number of initially fixed variables not at their initial bound)"
4949  << "\nq_C_hat = " << right << setw(int_w) << act_set.q_C_hat() << " (Number of initially fixed variables at the other bound)"
4950  << "\nq_D_hat = " << right << setw(int_w) << act_set.q_D_hat() << " (Number of initially fixed variables still fixed at initial bound)"
4951  << endl;
4952 
4953  // Print table of quantities in augmented KKT system
4954  out << "\nQuantities for augmentations to the initial KKT system:\n";
4955  const size_type q_hat = act_set.q_hat();
4956  out << endl
4957  << right << setw(int_w) << "s"
4958  << right << setw(int_w) << "ij_map(s)"
4959  << right << setw(int_w) << "bnd(s)"
4960  << right << setw(dbl_w) << "constr_norm(s)"
4961  << right << setw(dbl_w) << "d_hat(s)"
4962  << right << setw(dbl_w) << "z_hat(s)"
4963  << right << setw(dbl_w) << "p_z_hat(s)"
4964  << endl;
4965  out << right << setw(int_w) << int_ul
4966  << right << setw(int_w) << int_ul
4967  << right << setw(int_w) << int_ul
4968  << right << setw(dbl_w) << dbl_ul
4969  << right << setw(dbl_w) << dbl_ul
4970  << right << setw(dbl_w) << dbl_ul
4971  << right << setw(dbl_w) << dbl_ul
4972  << endl;
4973  {for( size_type s = 1; s <= q_hat; ++s ) {
4974  out << right << setw(int_w) << s
4975  << right << setw(int_w) << act_set.ij_map(s)
4976  << right << setw(int_w) << bnd_str(act_set.bnd(s))
4977  << right << setw(dbl_w) << act_set.constr_norm(s)
4978  << right << setw(dbl_w) << act_set.d_hat()(s)
4979  << right << setw(dbl_w) << act_set.z_hat()(s)
4980  << right << setw(dbl_w) << act_set.p_z_hat()(s)
4981  << endl;
4982  }}
4983 
4984  // Print P_XF_hat, P_FC_hat, P_plus_hat, U_hat and S_hat
4985  out << "\nP_XF_hat =\n" << act_set.P_XF_hat();
4986  out << "\nP_FC_hat =\n" << act_set.P_FC_hat();
4987  out << "\nP_plus_hat =\n" << act_set.P_plus_hat();
4988  out << "\nU_hat =\n" << act_set.U_hat();
4989  if(print_S_hat)
4990  out << "\nS_hat =\n" << act_set.S_hat();
4991 
4992  // Print table of multipliers for q_D_hat
4993  out << "\nQuantities for initially fixed variables which are still fixed at their initial bound:\n";
4994  const size_type q_D_hat = act_set.q_D_hat();
4995  out << endl
4996  << right << setw(int_w) << "k"
4997  << right << setw(int_w) << "l_fxfx(k)"
4998  << right << setw(dbl_w) << "mu_D_hat(k)"
4999  << right << setw(dbl_w) << "p_mu_D_hat(s)"
5000  << endl;
5001  out << right << setw(int_w) << int_ul
5002  << right << setw(int_w) << int_ul
5003  << right << setw(dbl_w) << dbl_ul
5004  << right << setw(dbl_w) << dbl_ul
5005  << endl;
5006  {for( size_type k = 1; k <= q_D_hat; ++k ) {
5007  out << right << setw(int_w) << k
5008  << right << setw(int_w) << act_set.l_fxfx(k)
5009  << right << setw(dbl_w) << act_set.mu_D_hat()(k)
5010  << right << setw(dbl_w) << act_set.p_mu_D_hat()(k)
5011  << endl;
5012  }}
5013 
5014  // Print Q_XD_hat
5015  out << "\nQ_XD_hat =\n" << act_set.Q_XD_hat();
5016 
5017  out << "\n*** End dump of current active set ***\n";
5018 }
5019 
5020 // QPSchurPack::QP
5021 
5022 void QPSchurPack::QP::dump_qp( std::ostream& out )
5023 {
5024  using std::endl;
5025  using std::setw;
5026  using std::left;
5027  using std::right;
5028  using Teuchos::Workspace;
5030 
5031  const Constraints
5032  &constraints = this->constraints();
5033 
5034  const size_type
5035  n = this->n(),
5036  n_R = this->n_R(),
5037  m = this->m(),
5038  m_breve = constraints.m_breve();
5039 
5040  out << "\n*** Original QP ***\n"
5041  << "\nn = " << n
5042  << "\nm = " << m
5043  << "\nm_breve = " << m_breve
5044  << endl;
5045  out << "\ng =\n" << g();
5046  out << "\nG =\n" << G();
5047  if(m) {
5048  out << "\nA =\n" << A();
5049  // Le'ts recover c from fo(n_R+1:n_R+m) = c - A' * Q_X * b_x
5050  throw std::logic_error(
5051  error_msg(__FILE__,__LINE__,"QPSchurPack::QP::dump_qp(...) : Error, "
5052  "m != not supported yet!"));
5053  // ToDo: Implement this when needed!
5054  }
5055  out << "\nA_bar =\n" << constraints.A_bar();
5056  // Get c_L_bar and c_U_bar
5057  DVector c_L_bar(n+m_breve), c_U_bar(n+m_breve);
5058  {for( size_type j = 1; j <= n+m_breve; ++j ){
5059  c_L_bar(j) = constraints.get_bnd(j,LOWER);
5060  c_U_bar(j) = constraints.get_bnd(j,UPPER);
5061  }}
5062  out << "\nc_L_bar =\n" << c_L_bar();
5063  out << "\nc_U_bar =\n" << c_U_bar();
5064 
5065  out << "\n*** Initial KKT system (fixed and free variables) ***\n"
5066  << "\nn_R = " << n_R
5067  << endl;
5068  out << "\nb_X =\n" << b_X();
5069  out << "\nQ_R =\n" << Q_R();
5070  out << "\nQ_X =\n" << Q_X();
5071  out << "\nKo =\n" << Ko();
5072  out << "\nfo =\n" << fo();
5073 }
5074 
5075 } // end namespace ConstrainedOptPack
virtual size_type n() const =0
bool drop_constraint(int jd, std::ostream *out, EOutputLevel output_level, bool force_refactorization=true, bool allow_any_cond=false)
Drop a constraint from the active set then refactorize the schur complement (if forced).
static value_type DEGENERATE_MULT
Value for near degenerate lagrange multipliers.
virtual const DVectorSlice fo() const =0
size_type s_map(int ij) const
Map from a constraint or initially fixed variable to a row and column in the schur complement S_bar...
virtual const MatrixOp & A_bar() const =0
int ij_map(size_type s) const
Returns -i for row & column of S_bar for an initially fixed variable left out of Ko that became free ...
QPSchur(const schur_comp_ptr_t &schur_comp=Teuchos::null, size_type max_iter=100, value_type max_real_runtime=1e+20, value_type feas_tol=1e-8, value_type loose_feas_tol=1e-6, value_type dual_infeas_tol=1e-12, value_type huge_primal_step=1e+20, value_type huge_dual_step=1e+20, value_type warning_tol=1e-10, value_type error_tol=1e-5, size_type iter_refine_min_iter=1, size_type iter_refine_max_iter=3, value_type iter_refine_opt_tol=1e-12, value_type iter_refine_feas_tol=1e-12, bool iter_refine_at_solution=true, bool salvage_init_schur_comp=true, MSADU::PivotTolerances pivot_tols=MSADU::PivotTolerances(1e-8, 1e-11, 1e-11))
virtual const x_init_t & x_init() const =0
Return the status of a variable initially.
void initialize(const MatrixSymOp *G, const MatrixOp *A, const MatrixOp *A_bar, const GenPermMatrixSlice *Q_R, const GenPermMatrixSlice *P_XF_hat, const GenPermMatrixSlice *P_plus_hat)
Initialize.
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
EBounds bnd(size_type s) const
Return which bound is active for the active constraint.
void Vp_StPtMtV(VectorMutable *v_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs2, BLAS_Cpp::Transp M_rhs2_trans, const Vector &v_rhs3, value_type beta=1.0)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual value_type get_bnd(size_type j, EBounds bnd) const =0
Return the bound for a constraint.
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
size_type q_F_hat() const
Return the number of variables that where initially fixed but are currently free or fixed to another ...
virtual void set_multipliers(const ActiveSet &act_set, const DVectorSlice &v, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve)
Map from the active set to the sparse multipliers for the inequality constraints. ...
MSADU::PivotTolerances pivot_tols() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
size_type q_C_hat() const
Return the number of variables that where initially fixed but are currently fixed to another bound...
const GenPermMatrixSlice & P_FC_hat() const
bool is_init_fixed(size_type j) const
Determine if a constriant was an initially fixed variable.
void refactorize_schur_comp()
Reinitialize the schur complement factorization for the current active set.
static void dump_act_set_quantities(const ActiveSet &act_set, std::ostream &out, bool print_S_hat=true)
Dump all the active set quantities for debugging.
virtual void pick_violated_policy(EPickPolicy pick_policy)=0
Set the policy used to pick a violated constraint.
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
virtual const l_x_X_map_t & l_x_X_map() const =0
Map from full x(i) to initially fixed x_X(l).
Represents and manages the active set for the QPSchur algorithm.
void syr2k(const MatrixOp &M, BLAS_Cpp::Transp M_trans, value_type alpha, const GenPermMatrixSlice &P1, BLAS_Cpp::Transp P1_trans, const GenPermMatrixSlice &P2, BLAS_Cpp::Transp P2_trans, value_type beta, MatrixSymOp *symwo_lhs)
void Mp_StPtMtP(MatrixOp *mwo_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs, const GenPermMatrixSlice &P_rhs2, BLAS_Cpp::Transp P_rhs2_trans)
virtual const DVectorSlice g() const =0
virtual const GenPermMatrixSlice & Q_X() const =0
(Q_X().ordered_by() == BY_ROW)
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
virtual size_type cols() const
bool add_constraint(size_type ja, EBounds bnd_ja, bool update_steps, std::ostream *out, EOutputLevel output_level, bool force_refactorization=true, bool allow_any_cond=false)
Add a constraint to the active set then refactorize the schur complemnt (if forced).
virtual const MatrixSymOpNonsing & Ko() const =0
virtual size_type m_breve() const =0
const GenPermMatrixSlice & P_plus_hat() const
const_iterator begin() const
Represents the QP to be solved by QPSchur {abstract}.
virtual void set_x(const ActiveSet &act_set, const DVectorSlice &v, DVectorSlice *x)
Set the values in x for all the variables.
size_t size_type
value_type constr_norm(size_type s) const
Returns ||a(j)||2 where j = ij_map(s).
virtual const MatrixOp & A() const =0
If m == 0 then don't call this, it may throw an exception or worse.
bool timeout_return(StopWatchPack::stopwatch *timer, std::ostream *out, EOutputLevel output_level) const
Determine if time has run out and if we should return.
size_type q_D_hat() const
Return the number of variables that where initially fixed and are still currently fixed to their inti...
virtual const GenPermMatrixSlice & Q_R() const =0
(Q_R().ordered_by() == BY_ROW)
virtual Constraints & constraints()=0
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
size_type q_plus_hat() const
Return the number of constraints from A_bar added to the active set.
virtual void ignore(size_type j)=0
Inform to ignore the jth constraint the next time pick_violated(...) is called.
virtual const i_x_X_map_t & i_x_X_map() const =0
Map from initially fixed x_X(l) to full x(i).
bool drop_add_constraints(int jd, size_type ja, EBounds bnd_ja, bool update_steps, std::ostream *out, EOutputLevel output_level)
Drop a constraint from, then add a constraint to the active set and refactorize the schur complement...
size_type q_hat() const
Return the total size of the schur complement.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
bool all_dof_used_up() const
Returns true if all the degrees of freedom of the QP are used up.
const ActiveSet & act_set() const
Return a reference to the active set object.
virtual size_type m() const =0
virtual ESolveReturn solve_qp(QP &qp, size_type num_act_change, const int ij_act_change[], const EBounds bnds[], std::ostream *out, EOutputLevel output_level, ERunTests test_what, DVectorSlice *x, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve, size_type *iter, size_type *num_adds, size_type *num_drops)
Solve a QP.
size_type l_fxfx(size_type k) const
Returns the indice of x_X(l) of the initially fixed variables that are still fixed at their original ...
virtual void dump_qp(std::ostream &out)
Dump the definition of the QP to a stream.
const GenPermMatrixSlice & Q_XD_hat() const
virtual ESolveReturn qp_algo(EPDSteps first_step, std::ostream *out, EOutputLevel output_level, ERunTests test_what, const DVectorSlice &vo, ActiveSet *act_set, DVectorSlice *v, DVectorSlice *x, size_type *iter, size_type *num_adds, size_type *num_drops, size_type *iter_refine_num_resid, size_type *iter_refine_num_solves, StopWatchPack::stopwatch *timer)
Run the algorithm from a dual feasible point.
ERunTests
Enumeration for if to run internal tests or not.
Represents the extra constraints in the QP to be satisfied by the schur complement QP solver QPSchur ...
const GenPermMatrixSlice & P_XF_hat() const
virtual const MatrixSymOp & G() const =0
EIterRefineReturn iter_refine(const ActiveSet &act_set, std::ostream *out, EOutputLevel output_level, const value_type ao, const DVectorSlice *bo, const value_type aa, const DVectorSlice *ba, DVectorSlice *v, DVectorSlice *z, size_type *iter_refine_num_resid, size_type *iter_refine_num_solves)
Perform iterative refinement on the augmented KKT system for the current active set.
void M_StMtInvMtM(MatrixSymOp *sym_gms_lhs, value_type alpha, const MatrixOp &mwo, BLAS_Cpp::Transp mwo_trans, const MatrixSymNonsing &mswof, MatrixSymNonsing::EMatrixDummyArg mwo_rhs)
Transp
void initialize(QP &qp, size_type num_act_change, const int ij_act_change[], const EBounds bnds[], bool test, bool salvage_init_schur_comp, std::ostream *out, EOutputLevel output_level)
Initialize with an additional active set.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
virtual void pick_violated(const DVectorSlice &x, size_type *j_viol, value_type *constr_val, value_type *viol_bnd_val, value_type *norm_2_constr, EBounds *bnd, bool *can_ignore) const =0
Pick a violated constraint.
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
virtual const DVectorSlice b_X() const =0
virtual size_type n_R() const =0
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)