MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_LineSearchWatchDog_Step.cpp
Go to the documentation of this file.
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <ostream>
45 #include <typeinfo>
46 
54 #include "ConstrainedOptPack/src/VectorWithNorms.h"
55 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
60 
61 namespace {
62  const int NORMAL_LINE_SEARCH = -1;
63 }
64 
65 namespace LinAlgOpPack {
67 }
68 
70  const direct_line_search_ptr_t& direct_line_search
71  , const merit_func_ptr_t& merit_func
72  , value_type eta
73  , value_type opt_kkt_err_threshold
74  , value_type feas_kkt_err_threshold
75  )
76  :
77  direct_line_search_(direct_line_search)
78  , merit_func_(merit_func)
79  , eta_(eta)
80  , opt_kkt_err_threshold_(opt_kkt_err_threshold)
81  , feas_kkt_err_threshold_(feas_kkt_err_threshold)
82  , watch_k_(NORMAL_LINE_SEARCH)
83 {}
84 
86  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
87 {
92 
93  using LinAlgOpPack::Vp_V;
94  using LinAlgOpPack::V_MtV;
95 
97 
98  NLPAlgo &algo = rsqp_algo(_algo);
99  NLPAlgoState &s = algo.rsqp_state();
100  NLP &nlp = algo.nlp();
101 
102  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
103  std::ostream& out = algo.track().journal_out();
104  out << std::boolalpha;
105 
106  // print step header.
107  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
109  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
110  }
111 
112  // /////////////////////////////////////////
113  // Set references to iteration quantities
114  //
115  // Set k+1 first then go back to get k to ensure
116  // we have backward storage.
117 
118  DVector
119  &x_kp1 = s.x().set_k(+1).v();
120  value_type
121  &f_kp1 = s.f().set_k(+1);
122  DVector
123  &c_kp1 = s.c().set_k(+1).v();
124 
125  const value_type
126  &f_k = s.f().get_k(0);
127  const DVector
128  &c_k = s.c().get_k(0).v();
129  const DVector
130  &x_k = s.x().get_k(0).v();
131  const DVector
132  &d_k = s.d().get_k(0).v();
133  value_type
134  &alpha_k = s.alpha().get_k(0);
135 
136  // /////////////////////////////////////
137  // Compute Dphi_k, phi_kp1 and phi_k
138 
139  // Dphi_k
140  const value_type
141  Dphi_k = merit_func().deriv();
142  if( Dphi_k >= 0 ) {
143  throw LineSearchFailure( "LineSearch2ndOrderCorrect_Step::do_step(...) : "
144  "Error, d_k is not a descent direction for the merit function " );
145  }
146 
147  // ph_kp1
148  value_type
149  &phi_kp1 = s.phi().set_k(+1) = merit_func().value( f_kp1, c_kp1 );
150 
151  // Must compute phi(x) at the base point x_k since the penalty parameter may have changed.
152  const value_type
153  &phi_k = s.phi().set_k(0) = merit_func().value( f_k, c_k );
154 
155  // //////////////////////////////////////
156  // Setup the calculation merit function
157 
158  // Here f_kp1, and c_kp1 are updated at the same time the
159  // line search is being performed.
160  nlp.set_f( &f_kp1 );
161  nlp.set_c( &c_kp1 );
162  MeritFuncCalcNLP
163  phi_calc( &merit_func(), &nlp );
164 
165  // ////////////////////////////////
166  // Use Watchdog near the solution
167 
168  if( watch_k_ == NORMAL_LINE_SEARCH ) {
169  const value_type
170  opt_kkt_err_k = s.opt_kkt_err().get_k(0),
171  feas_kkt_err_k = s.feas_kkt_err().get_k(0);
172  if( opt_kkt_err_k <= opt_kkt_err_threshold() && feas_kkt_err_k <= feas_kkt_err_threshold() ) {
173  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
174  out << "\nopt_kkt_err_k = " << opt_kkt_err_k << " <= opt_kkt_err_threshold = "
175  << opt_kkt_err_threshold() << std::endl
176  << "\nfeas_kkt_err_k = " << feas_kkt_err_k << " <= feas_kkt_err_threshold = "
177  << feas_kkt_err_threshold() << std::endl
178  << "\nSwitching to watchdog linesearch ...\n";
179  }
180  watch_k_ = 0;
181  }
182  }
183 
184  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
185  out << "\nTrial point:\n"
186  << "phi_k = " << phi_k << std::endl
187  << "Dphi_k = " << Dphi_k << std::endl
188  << "phi_kp1 = " << phi_kp1 << std::endl;
189  }
190 
191  bool ls_success = true,
192  step_return = true;
193 
194  switch( watch_k_ ) {
195  case 0:
196  {
197  // Take a full step
198  const value_type phi_cord = phi_k + eta() * Dphi_k;
199  const bool accept_step = phi_kp1 <= phi_cord;
200 
201  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
202  out << "\n*** Zeroth watchdog iteration:\n"
203  << "\nphi_kp1 = " << phi_kp1 << ( accept_step ? " <= " : " > " )
204  << "phi_k + eta * Dphi_k = " << phi_cord << std::endl;
205  }
206 
207  if( phi_kp1 > phi_cord ) {
208  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
209  out << "\nAccept this increase for now but watch out next iteration!\n";
210  }
211  // Save this initial point
212  xo_ = x_k;
213  fo_ = f_k;
214  nrm_co_ = norm_inf( c_k );
215  do_ = d_k;
216  phio_ = phi_k;
217  Dphio_ = Dphi_k;
218  phiop1_ = phi_kp1;
219  // Slip the update of the penalty parameter
220  const value_type mu_k = s.mu().get_k(0);
221  s.mu().set_k(+1) = mu_k;
222  // Move on to the next step in the watchdog procedure
223  watch_k_ = 1;
224  }
225  else {
226  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
227  out << "\nAll is good!\n";
228  }
229  // watch_k_ stays 0
230  }
231  step_return = true;
232  break;
233  }
234  case 1:
235  {
236  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
237  out << "\n*** First watchdog iteration:\n"
238  << "\nDo a line search to determine x_kp1 = x_k + alpha_k * d_k ...\n";
239  }
240  // Now do a line search but and we require some type of reduction
241  const DVectorSlice xd[2] = { x_k(), d_k() };
242  MeritFuncCalc1DQuadratic phi_calc_1d( phi_calc, 1, xd, &x_kp1() );
243  ls_success = direct_line_search().do_line_search( phi_calc_1d, phi_k
244  , &alpha_k, &phi_kp1
245  , (int)olevel >= (int)PRINT_ALGORITHM_STEPS ?
246  &out : static_cast<std::ostream*>(0) );
247 
248  // If the linesearch failed then the rest of the tests will catch this.
249 
250  value_type phi_cord = 0;
251  bool test1, test2;
252 
253  if( ( test1 = ( phi_k <= phio_ ) )
254  || ( test2 = phi_kp1 <= ( phi_cord = phio_ + eta() * Dphio_ ) ) )
255  {
256  // We will accept this step and and move on.
257  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
258  out
259  << "\nphi_k = " << phi_k << ( test1 ? " <= " : " > " )
260  << "phi_km1 = " << phio_ << std::endl
261  << "phi_kp1 = " << phi_kp1 << ( test2 ? " <= " : " > " )
262  << "phi_km1 + eta * Dphi_km1 = " << phi_cord << std::endl
263  << "This is a sufficent reduction so reset watchdog.\n";
264  }
265  watch_k_ = 0;
266  step_return = true;
267  }
268  else if ( ! ( test1 = ( phi_kp1 <= phio_ ) ) ) {
269  // Even this reduction is no good!
270  // Go back to original point and do a linesearch from there.
271  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
272  out
273  << "\nphi_kp1 = " << phi_kp1 << " > phi_km1 = " << phio_ << std::endl
274  << "This is not good reduction in phi so do linesearch from x_km1\n"
275  << "\n* Go back to x_km1: x_kp1 = x_k - alpha_k * d_k ...\n";
276  }
277 
278  // Go back from x_k to x_km1 for iteration k:
279  //
280  // x_kp1 = x_km1
281  // x_kp1 = x_k - alpha_km1 * d_km1
282  //
283  // A negative sign for alpha is an indication that we are backtracking.
284  //
285  s.alpha().set_k(0) = -1.0;
286  s.d().set_k(0).v() = do_;
287  s.f().set_k(+1) = fo_;
288 
289  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
290  out << "Output iteration k ...\n"
291  << "k = k+1\n";
292  }
293 
294  // Output these iteration quantities
295  algo.track().output_iteration( algo ); // k
296  // Transition to iteration k+1
297  s.next_iteration();
298 
299  // Take the step from x_k = x_km2 to x_kp1 for iteration k (k+1):
300  //
301  // x_kp1 = x_km2 + alpha_n * d_km2
302  // x_kp1 = x_k + alpha_n * d_km1
303  // x_kp1 = x_n
304  //
305  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
306  out << "\n* Take the step from x_k = x_km2 to x_kp1 for iteration k (k+1)\n"
307  << "Find: x_kp1 = x_k + alpha_k * d_k = x_km2 + alpha_k * d_km2\n ...\n";
308  }
309 
310  // alpha_k = 1.0
311  value_type &alpha_k = s.alpha().set_k(0) = 1.0;
312 
313  // /////////////////////////////////////
314  // Compute Dphi_k and phi_k
315 
316  // x_k
317  const DVector &x_k = xo_;
318 
319  // d_k
320  const DVector &d_k = s.d().set_k(0).v() = do_;
321 
322  // Dphi_k
323  const value_type &Dphi_k = Dphio_;
324 
325  // phi_k
326  const value_type &phi_k = s.phi().set_k(0) = phio_;
327 
328  // Here f_kp1, and c_kp1 are updated at the same time the
329  // line search is being performed.
330  algo.nlp().set_f( &s.f().set_k(+1) );
331  algo.nlp().set_c( &s.c().set_k(+1).v() );
332  phi_calc.set_nlp( algo.get_nlp() );
333 
334  // ////////////////////////////////////////
335  // Compute x_xp1 and ph_kp1 for full step
336 
337  // x_kp1 = x_k + alpha_k * d_k
338  DVector &x_kp1 = s.x().set_k(+1).v();
339  V_VpV( &x_kp1, x_k, d_k );
340 
341  // phi_kp1
342  value_type &phi_kp1 = s.phi().set_k(+1) = phiop1_;
343 
344  const DVectorSlice xd[2] = { x_k(), d_k() };
345  MeritFuncCalc1DQuadratic phi_calc_1d( phi_calc, 1, xd, &x_kp1() );
346  ls_success = direct_line_search().do_line_search(
347  phi_calc_1d, phi_k
348  , &alpha_k, &phi_kp1
349  , (int)olevel >= (int)PRINT_ALGORITHM_STEPS ?
350  &out : static_cast<std::ostream*>(0) );
351 
352  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
353  out << "\nOutput iteration k (k+1) ...\n"
354  << "k = k+1 (k+2)\n"
355  << "Reinitialize watchdog algorithm\n";
356  }
357 
358  // Output these iteration quantities
359  algo.track().output_iteration( algo ); // (k+1)
360  // Transition to iteration k+1 (k+2)
361  s.next_iteration();
362 
363  watch_k_ = 0; // Reinitialize the watchdog
364 
365  // Any update for k (k+2) should use the last updated value
366  // which was for k-2 (k) since there is not much info for k-1 (k+1).
367  // Be careful here and make sure this is square with other steps.
368 
369  algo.do_step_next( EvalNewPoint_name );
370  step_return = false; // Redirect control
371  }
372  else {
373  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
374  out
375  << "phi_kp1 = " << phi_kp1 << " <= phi_km1 = " << phio_ << std::endl
376  << "\nAccept this step but do a linesearch next iteration!\n";
377  }
378  // Slip the update of the penalty parameter
379  const value_type mu_k = s.mu().get_k(0);
380  s.mu().set_k(+1) = mu_k;
381  // Do the last stage of the watchdog procedure next iteration.
382  watch_k_ = 2;
383  step_return = true;
384  }
385  break;
386  }
387  case NORMAL_LINE_SEARCH:
388  case 2:
389  {
390  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
391  if( watch_k_ == 2 ) {
392  out << "\n*** Second watchdog iteration:\n"
393  << "Do a line search to determine x_kp1 = x_k + alpha_k * d_k ...\n";
394  }
395  else {
396  out << "\n*** Normal linesearch:\n"
397  << "Do a line search to determine x_kp1 = x_k + alpha_k * d_k ...\n";
398  }
399  }
400 
401  const DVectorSlice xd[2] = { x_k(), d_k() };
402  MeritFuncCalc1DQuadratic phi_calc_1d( phi_calc, 1, xd, &x_kp1() );
403  ls_success = direct_line_search().do_line_search( phi_calc_1d, phi_k
404  , &alpha_k, &phi_kp1
405  , (int)olevel >= (int)PRINT_ALGORITHM_STEPS ?
406  &out : static_cast<std::ostream*>(0) );
407 
408  if( watch_k_ == 2 )
409  watch_k_ = 0;
410 
411  step_return = true;
412  break;
413  }
414  default:
415  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local programming error
416  }
417 
418  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
419  out << "\nalpha = " << s.alpha().get_k(0) << "\n";
420  out << "\nphi_kp1 = " << s.phi().get_k(+1) << "\n";
421  }
422 
423  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
424  out << "\nd_k = \n" << s.d().get_k(0)();
425  out << "\nx_kp1 = \n" << s.x().get_k(+1)();
426  }
427 
428  if( !ls_success )
429  throw LineSearchFailure("LineSearchWatchDog_Step::do_step(): Line search failure");
430 
431  return step_return;
432 
433 }
434 
435 void MoochoPack::LineSearchWatchDog_Step::print_step( const Algorithm& algo
436  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
437  , std::ostream& out, const std::string& L ) const
438 {
439  out << L << "*** Use the Watchdog linesearch when near solution.\n"
440  << L << "default: opt_kkt_err_threshold = 0.0\n"
441  << L << " feas_kkt_err_threshold = 0.0\n"
442  << L << " eta = 1.0e-4\n"
443  << L << " watch_k = NORMAL_LINE_SEARCH\n"
444  << L << "begin definition of NLP merit function phi.value(f(x),c(x)):\n";
445 
446  merit_func().print_merit_func( out, L + " " );
447 
448  out << L << "end definition\n"
449  << L << "Dphi_k = phi.deriv()\n"
450  << L << "if Dphi_k >= 0 then\n"
451  << L << " throw line_search_failure\n"
452  << L << "end\n"
453  << L << "phi_kp1 = phi_k.value(f_kp1,c_kp1)\n"
454  << L << "phi_k = phi.value(f_k,c_k)\n"
455  << L << "if watch_k == NORMAL_LINE_SEARCH then\n"
456  << L << " if opt_kkt_err <= opt_kkt_err_threshold\n"
457  << L << " and feas_kkt_err <= feas_kkt_err_threshold then\n"
458  << L << " *** Start using watchdog from now on!\n"
459  << L << " watch_k = 0\n"
460  << L << " end\n"
461  << L << "end\n"
462  << L << "if watch_k == 0 then\n"
463  << L << " *** Zeroth watchdog iteration\n"
464  << L << " if phi_kp1 >= phi_k + eta * Dphi_k then\n"
465  << L << " *** Accept this increase for now but watch out next iteration!\n"
466  << L << " *** Save the first point\n"
467  << L << " xo = x_k\n"
468  << L << " fo = f_k\n"
469  << L << " nrm_co = norm_inf_c_k\n"
470  << L << " do = d_k\n"
471  << L << " phio = phi_k\n"
472  << L << " Dphio = Dphi_k\n"
473  << L << " phiop1 = phi_kp1\n"
474  << L << " *** Skip the update of the penalty parameter next iteration.\n"
475  << L << " mu_kp1 = mu_k\n"
476  << L << " *** Continue with next step in watchdog\n"
477  << L << " watch_k = 1\n"
478  << L << " else\n"
479  << L << " *** This is a good step so take it!\n"
480  << L << " end\n"
481  << L << "else if watch_k == 1 then\n"
482  << L << " *** First watchdog iteration\n"
483  << L << " Do line search for: x_kp1 = x_k + alpha_k + d_k\n"
484  << L << " -> alpha_k, x_kp1, f_kp1, c_kp1, phi_kp1\n"
485  << L << " if ( phi_k <= phio ) or ( phi_kp1 <= phio + eta * Dphio ) then\n"
486  << L << " *** We will accept this step and reinitialize the watchdog\n"
487  << L << " watch_k = 0\n"
488  << L << " else if ( phi_kp1 > phio ) then\n"
489  << L << " *** This reduction is no good!\n"
490  << L << " *** Go back from x_k to x_km1 for this iteration (k)\n"
491  << L << " alpha_k = -1.0\n"
492  << L << " d_k = do\n"
493  << L << " f_kp1 = fo\n"
494  << L << " Output this iteration (k)\n"
495  << L << " k = k+1\n"
496  << L << " *** Go from x_k = x_km2 to x_kp1 for this iteration (k+1)\n"
497  << L << " alpha_k = 1\n"
498  << L << " x_k = xo\n"
499  << L << " d_k = do\n"
500  << L << " Dphi_k = Dphio\n"
501  << L << " phi_k = phio\n"
502  << L << " x_kp1 = x_k + d_k\n"
503  << L << " phi_kp1 = phiop1\n"
504  << L << " Do line search for: x_kp1 = x_k + alpha_k + d_k\n"
505  << L << " -> alpha_k, x_kp1, f_kp1, c_kp1, phi_kp1\n"
506  << L << " Output this iteration (k+1)\n"
507  << L << " k = k+1\n"
508  << L << " *** Any updates for k (k+2) should use the last updated value\n"
509  << L << " *** which was for k-2 (k) since there is not much info for k-1 (k+1).\n"
510  << L << " *** Be careful here and make sure this works with other steps.\n"
511  << L << " goto EvalNewPoint\n"
512  << L << " else\n"
513  << L << " *** Accept this reduction but do a linesearch next iteration!\n"
514  << L << " *** Skip the update of the penalty parameter next iteration.\n"
515  << L << " mu_kp1 = mu_k\n"
516  << L << " *** Continue with next step in watchdog\n"
517  << L << " watch_k = 2\n"
518  << L << " end\n"
519  << L << "else if ( watch_k == 2 ) then\n"
520  << L << " *** Second watchdog iteration\n"
521  << L << " Do line search for: x_kp1 = x_k + alpha_k + d_k\n"
522  << L << " -> alpha_k, x_kp1, f_kp1, c_kp1, phi_kp1\n"
523  << L << " *** Reset the watchdog algorithm\n"
524  << L << " watch_k = 1\n"
525  << L << "else if ( watch_k == NORMAL_LINE_SEARCH ) then\n"
526  << L << " Do line search for: x_kp1 = x_k + alpha_k + d_k\n"
527  << L << " -> alpha_k, x_kp1, f_kp1, c_kp1, phi_kp1\n"
528  << L << " begin direct line search : \""
529  << typeName(direct_line_search()) << "\"\n";
530 
531  direct_line_search().print_algorithm( out, L + " " );
532 
533  out
534  << L << " end direct line search\n"
535  << L << "end\n"
536  << L << "if maximum number of linesearch iterations are exceeded then\n"
537  << L << " throw line_search_failure\n"
538  << L << "end\n";
539 }
540 
541 #endif // 0
LineSearchWatchDog_Step(const direct_line_search_ptr_t &direct_line_search=0, const merit_func_ptr_t &merit_func=0, value_type eta=1e-4, value_type opt_kkt_err_threshold=1e-1, value_type feas_kkt_err_threshold=1e-3)
std::string typeName(const T &t)
void Vp_StV(DVectorSlice *vs_lhs, value_type alpha, const DVectorSlice &vs_rhs)
vs_lhs += alpha * vs_rhs (BLAS xAXPY)
void V_VpV(DVector *v_lhs, const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
v_lhs = alpha (elementwise)
void print_vector_change_stats(const DVectorSlice &x, const char x_name[], const DVectorSlice &d, const char d_name[], std::ostream &out)
Compute statistics for change in a vector and output to a stream.
EJournalOutputLevel
enum for journal output.
void V_VpV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
std::ostream * out
void print_algorithm_step(const Algorithm &algo, Algorithm::poss_type step_poss, EDoStepType type, Algorithm::poss_type assoc_step_poss, std::ostream &out)
Prints to 'out' the algorithm step.
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)
v_lhs = alpha * op(M_rhs1) * v_rhs2 + beta * v_lhs (BLAS xGEMV)
value_type norm_inf(const SparseVectorSlice< T_Ele > &sv_rhs)
result = ||sv_rhs||inf (BLAS IxAMAX)
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
value_type norm_inf(const DVectorSlice &vs_rhs)
result = ||vs_rhs||infinity (BLAS IxAMAX)
NLPAlgo & rsqp_algo(Algorithm &algo)
Convert from a Algorithm to a NLPAlgo.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void Vp_V(VectorMutable *v_lhs, const V &V_rhs)
v_lhs += V_rhs.