MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConstrainedOptPack_ConstraintsRelaxedStd.cpp
Go to the documentation of this file.
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 // ToDo: 12/29/00: Consider scaling when determining if a
43 // constraint is violated. We should consider the size of
44 // ||d(e[j](x))/d(x)||inf but this is expensive to compute
45 // given the current interfaces. We really need to rectify
46 // this!
47 //
48 
49 #include <assert.h>
50 
51 #include <limits>
52 
64 #include "Teuchos_Assert.hpp"
65 
66 namespace {
67 
69 convert_bnd_type( int bnd_type )
70 {
71  switch(bnd_type) {
72  case -1:
74  case 0:
76  case +1:
78  default:
80  }
81  return ConstrainedOptPack::LOWER; // Never be called
82 }
83 
84 // Get an element from a sparse vector and return zero if it does not exist
85 AbstractLinAlgPack::value_type get_sparse_element(
88  )
89 {
91  *ele_ptr = v.lookup_element(i);
92  return ele_ptr ? ele_ptr->value() : 0.0;
93 }
94 
95 } // end namespace
96 
97 namespace ConstrainedOptPack {
98 namespace QPSchurPack {
99 
100 // members for ConstraintsRelaxedStd
101 
103  :inequality_pick_policy_(ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY)
104  ,etaL_(0.0)
105  ,dL_(NULL)
106  ,dU_(NULL)
107  ,eL_(NULL)
108  ,eU_(NULL)
109  ,Ed_(NULL)
110  ,last_added_j_(0)
111  ,last_added_bound_(0.0)
112  ,last_added_bound_type_(FREE)
113  ,next_undecomp_f_k_(0)
114 {}
115 
117  const VectorSpace::space_ptr_t &space_d_eta
118  ,value_type etaL
119  ,const Vector *dL
120  ,const Vector *dU
121  ,const MatrixOp *E
122  ,BLAS_Cpp::Transp trans_E
123  ,const Vector *b
124  ,const Vector *eL
125  ,const Vector *eU
126  ,const MatrixOp *F
127  ,BLAS_Cpp::Transp trans_F
128  ,const Vector *f
129  ,size_type m_undecomp
130  ,const size_type j_f_undecomp[]
131  ,VectorMutable *Ed
132  ,bool check_F
133  ,value_type bounds_tol
134  ,value_type inequality_tol
135  ,value_type equality_tol
136  )
137 {
138  size_type
139  nd = space_d_eta->dim() - 1,
140  m_in = 0,
141  m_eq = 0;
142 
143  TEUCHOS_TEST_FOR_EXCEPT( !( m_undecomp == (F ? f->dim() : 0) ) ); // ToDo: support decomposed equalities in future.
144 
145  // Validate that the correct sets of constraints are selected
147  dL && !dU, std::invalid_argument
148  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
149  "If dL!=NULL then dU!=NULL must also be true." );
150  TEUCHOS_TEST_FOR_EXCEPTION(
151  E && ( !b || !eL || !eU ), std::invalid_argument
152  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
153  "If E!=NULL then b!=NULL, eL!=NULL and eU!=NULL must also be true." );
154  TEUCHOS_TEST_FOR_EXCEPTION(
155  F && !f, std::invalid_argument
156  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
157  "If F!=NULL then f!=NULL must also be true." );
158 
159  // Validate input argument sizes
160  if(dL) {
161  const size_type dL_dim = dL->dim(), dU_dim = dU->dim();
162  TEUCHOS_TEST_FOR_EXCEPTION(
163  dL_dim != nd, std::invalid_argument
164  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
165  "dL.dim() != d->dim()." );
166  TEUCHOS_TEST_FOR_EXCEPTION(
167  dU_dim != nd, std::invalid_argument
168  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
169  "dU.dim() != d->dim()." );
170  }
171  if(E) {
172  const size_type
173  E_rows = E->rows(),
174  E_cols = E->cols(),
175  opE_cols = BLAS_Cpp::cols( E_rows, E_cols, trans_E ),
176  b_dim = b->dim(),
177  eL_dim = eL->dim(),
178  eU_dim = eU->dim(),
179  Ed_dim = Ed ? Ed->dim() : 0;
180  m_in = BLAS_Cpp::rows( E_rows, E_cols, trans_E );
181  TEUCHOS_TEST_FOR_EXCEPTION(
182  opE_cols != nd, std::invalid_argument
183  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
184  "op(E).cols() != nd." );
185  TEUCHOS_TEST_FOR_EXCEPTION(
186  b_dim != m_in, std::invalid_argument
187  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
188  "b->dim() != op(E).rows()." );
189  TEUCHOS_TEST_FOR_EXCEPTION(
190  eL_dim != m_in, std::invalid_argument
191  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
192  "eL->dim() != op(E).rows()." );
193  TEUCHOS_TEST_FOR_EXCEPTION(
194  eU_dim != m_in, std::invalid_argument
195  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
196  "eU->dim() != op(E).rows()." );
197  TEUCHOS_TEST_FOR_EXCEPTION(
198  Ed && Ed_dim != m_in, std::invalid_argument
199  ,"ConstraintsRelaxedStd::initialize(...) : Error, "
200  "Ed->dim() != op(E).rows()." );
201  }
202  if(F) {
203  const size_type
204  F_rows = F->rows(),
205  F_cols = F->cols(),
206  opF_cols = BLAS_Cpp::cols( F_rows, F_cols, trans_F ),
207  f_dim = f->dim();
208  m_eq = BLAS_Cpp::rows( F_rows, F_cols, trans_F );
209  TEUCHOS_TEST_FOR_EXCEPTION(
210  opF_cols != nd, std::invalid_argument
211  ,"QPSolverRelaxed::solve_qp(...) : Error, "
212  "op(F).cols() != nd." );
213  TEUCHOS_TEST_FOR_EXCEPTION(
214  f_dim != m_eq, std::invalid_argument
215  ,"QPSolverRelaxed::solve_qp(...) : Error, "
216  "f->dim() != op(F).rows()." );
217  }
218 
219  // Initialize other members
221  space_d_eta,m_in,m_eq,E,trans_E,b,F,trans_F,f,m_undecomp,j_f_undecomp);
222  etaL_ = etaL;
223  dL_ = dL;
224  dU_ = dU;
225  eL_ = eL;
226  eU_ = eU;
227  Ed_ = Ed;
228  check_F_ = check_F;
229  bounds_tol_ = bounds_tol;
230  inequality_tol_ = inequality_tol;
231  equality_tol_ = equality_tol;
232  last_added_j_ = 0; // No cached value.
233  next_undecomp_f_k_ = m_undecomp ? 1 : 0; // Check the first undecomposed equality
234 
235 }
236 
239 {
240  return A_bar_;
241 }
242 
243 // Overridden from Constraints
244 
246 {
247  return A_bar_.nd() + 1;
248 }
249 
251 {
252  return A_bar_.m_in() + A_bar_.m_eq();
253 }
254 
255 const MatrixOp& ConstraintsRelaxedStd::A_bar() const
256 {
257  return A_bar_;
258 }
259 
261 {
262  switch(pick_policy) {
263  case ANY_VIOLATED:
264  inequality_pick_policy_ = ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY;
265  break;
266  case MOST_VIOLATED:
267  inequality_pick_policy_ = ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY;
268  break;
269  default:
271  }
272 }
273 
276 {
277  switch(inequality_pick_policy_) {
279  return ANY_VIOLATED;
281  return ANY_VIOLATED;
283  return MOST_VIOLATED;
284  default:
286  }
287  return ANY_VIOLATED; // will never be executed
288 }
289 
291  const DVectorSlice& x_in, size_type* j_viol, value_type* constr_val
292  ,value_type* viol_bnd_val, value_type* norm_2_constr, EBounds* bnd, bool* can_ignore
293  ) const
294 {
295  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
297  using LinAlgOpPack::V_VmV;
298 
300  x_in.dim() != A_bar_.nd()+1, std::length_error
301  ,"ConstraintsRelaxedStd::pick_violated(...) : Error, "
302  "x is the wrong length" );
303 
304  const size_type
305  nd = A_bar_.nd();
306 
307  // Get a version of x = [ d; eta ] in the correct vector object
308  VectorSpace::vec_mut_ptr_t
309  x = A_bar_.space_cols().create_member();
310  (VectorDenseMutableEncap(*x)()) = x_in;
311  VectorSpace::vec_mut_ptr_t
312  d = x->sub_view(1,nd);
313  const value_type
314  eta = x->get_ele(nd+1);
315 
316  bool Ed_computed = false;
317 
318  // //////////////////////////////////////////////
319  // Check the equality constraints first
320  if( check_F_ && A_bar_.F() ) {
321  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Update below code!
322 /*
323  // The basic strategy here is to go through all of the equality
324  // constraints first and add all of the ones that are violated by
325  // more than the set tolerance. Those that are not sufficiently
326  // violated are passed over but they are remembered for later.
327  // Only when all of the constraints have been gone through once
328  // will those passed over constraints be considered and then only
329  // if ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY is selected.
330  const GenPermMatrixSlice& P_u = A_bar_.P_u();
331  size_type e_k_mat_row, e_k_mat_col = 1; // Mapping matrix for e(j)
332  GenPermMatrixSlice e_k_mat; // ToDo: Change to EtaVector
333  if( next_undecomp_f_k_ <= P_u.nz() ) {
334  // This is our first pass through the undecomposed equalities.
335  GenPermMatrixSlice::const_iterator P_u_itr, P_u_end; // only used if P_u is not identity
336  if( !P_u.is_identity() ) {
337  P_u_itr = P_u.begin() + (next_undecomp_f_k_ - 1);
338  P_u_end = P_u.end();
339  }
340  for( ; next_undecomp_f_k_ <= P_u.nz(); ) {
341  size_type k = 0;
342  if( !P_u.is_identity() ) {
343  k = P_u_itr->row_i();
344  ++P_u_itr;
345  }
346  else {
347  k = next_undecomp_f_k_;
348  }
349  ++next_undecomp_f_k_;
350  // evaluate the constraint e(k)'*[op(F)*d + (1-eta)*f]
351  value_type
352  Fd_k = 0.0,
353  f_k = (*A_bar_.f())(k);
354  e_k_mat.initialize(
355  A_bar_.m_eq(),1,1,0,0,GPMSIP::BY_ROW_AND_COL
356  ,&(e_k_mat_row=k),&e_k_mat_col,false );
357  DVectorSlice Fd_k_vec(&Fd_k,1);
358  AbstractLinAlgPack::Vp_StPtMtV( // ToDo: Use transVtMtV(...) instead!
359  &Fd_k_vec, 1.0, e_k_mat, BLAS_Cpp::trans
360  ,*A_bar_.F(), A_bar_.trans_F(), d, 0.0 );
361  const value_type
362  err = ::fabs(Fd_k + (1.0 - eta)*f_k) / (1.0 + ::fabs(f_k));
363  if( err > equality_tol() ) {
364  *j_viol = nd + 1 + A_bar_.m_in() + k;
365  *constr_val = Fd_k - eta*f_k;
366  *norm_2_constr = 1.0; // ToDo: Compute it some how?
367  *bnd = EQUALITY;
368  *viol_bnd_val = -f_k;
369  *can_ignore = false; // Given this careful algorithm this should be false
370  // cache this
371  last_added_j_ = *j_viol;
372  last_added_bound_type_ = *bnd;
373  last_added_bound_ = *viol_bnd_val;
374  return;
375  }
376  else {
377  passed_over_equalities_.push_back(k); // remember for later
378  }
379  }
380  }
381  else if(
382  inequality_pick_policy() == ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY
383  && passed_over_equalities_.size() > 0
384  )
385  {
386  // Now look through the list of passed over equalities and see which one
387  // is violated. If a passed over constraint is added to the active set
388  // then it is removed from this list.
389  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
390  }
391 */
392  }
393 
394  // /////////////////////////////////////////////
395  // Find the most violated variable bound.
396 
397  size_type max_bound_viol_j = 0;
398  value_type max_bound_viol = 0.0;
399  value_type max_bound_d_viol = 0.0;
400  value_type max_bound_dLU_viol = 0.0;
401  int max_bound_viol_type = -2;
402  if( dL_ && ( dL_->nz() || dU_->nz() ) ) {
403  // dL <= d <= dU
405  *d, *dL_, *dU_
406  ,&max_bound_viol_j, &max_bound_viol
407  ,&max_bound_d_viol, &max_bound_viol_type, &max_bound_dLU_viol
408  );
409  if( max_bound_viol > bounds_tol_ ) {
410  // Set the return values
411  *j_viol = max_bound_viol_j;
412  *constr_val = max_bound_d_viol;
413  *norm_2_constr = 1.0; // This is correct
414  *bnd = convert_bnd_type(max_bound_viol_type);
415  *viol_bnd_val = max_bound_dLU_viol;
416  *can_ignore = false;
417  }
418  else {
419  max_bound_viol_j = 0; // No variable bounds sufficiently violated.
420  }
421  }
422 
423  if( ( inequality_pick_policy_ == ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY
424  || inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY )
425  && max_bound_viol_j
426  )
427  {
428  // A variable bound has been violated so lets just return it!
429  last_added_j_ = *j_viol;
430  last_added_bound_type_ = *bnd;
431  last_added_bound_ = *viol_bnd_val;
432  return;
433  }
434 
435  // /////////////////////////////////////////////
436  // Check the general inequalities
437 
438  size_type max_inequality_viol_j = 0;
439  value_type max_inequality_viol = 0.0;
440  value_type max_inequality_e_viol = 0.0;
441  value_type max_inequality_eLU_viol = 0.0;
442  int max_inequality_viol_type = -2;
443 
444  if( inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY ) {
445  // Find the first general inequality violated by more than
446  // the defined tolerance.
448  true, std::logic_error
449  ,"ConstraintsRelaxedStd::pick_violated(...) : Error,\n"
450  "The option ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY has not been implemented yet\n" );
451  }
452  else {
453  // Find the most violated inequality constraint
454  if( A_bar_.m_in() && ( eL_->nz() || eU_->nz() ) ) {
455  // e = op(E)*d - b*eta
456  VectorSpace::vec_mut_ptr_t e = eL_->space().create_member();
457  LinAlgOpPack::V_MtV( e.get(), *A_bar_.E(), A_bar_.trans_E(), *d );
458  if(Ed_) {
459  *Ed_ = *e;
460  Ed_computed = true;
461  }
462  LinAlgOpPack::Vp_StV( e.get(), -eta, *A_bar_.b() );
463  // eL <= e <= eU
465  *e, *eL_, *eU_
466  ,&max_inequality_viol_j, &max_inequality_viol
467  ,&max_inequality_e_viol, &max_inequality_viol_type, &max_inequality_eLU_viol
468  );
469  if( max_inequality_viol > max_bound_viol
470  && max_inequality_viol > inequality_tol_ )
471  {
472  *j_viol = max_inequality_viol_j + nd + 1; // offset into A_bar
473  *constr_val = max_inequality_e_viol;
474  *norm_2_constr = 1.0; // This is not correct!
475  *bnd = convert_bnd_type(max_inequality_viol_type);
476  *viol_bnd_val = max_inequality_eLU_viol;
477  *can_ignore = false;
478  }
479  else {
480  max_inequality_viol_j = 0; // No general inequality constraints sufficiently violated.
481  }
482  }
483  }
484 
485  if( max_bound_viol_j || max_inequality_viol_j ) {
486  // One of the constraints has been violated so just return it.
487  last_added_j_ = *j_viol;
488  last_added_bound_type_ = *bnd;
489  last_added_bound_ = *viol_bnd_val;
490  return;
491  }
492 
493  // If we get here then no constraint was found that violated any of the tolerances.
494  if( Ed_ && !Ed_computed ) {
495  // Ed = op(E)*d
497  }
498  *j_viol = 0;
499  *constr_val = 0.0;
500  *viol_bnd_val = 0.0;
501  *norm_2_constr = 0.0;
502  *bnd = FREE; // Meaningless
503  *can_ignore = false; // Meaningless
504 }
505 
507 {
509  true, std::logic_error
510  ,"ConstraintsRelaxedStd::ignore(...) : Error, "
511  "This operation is not supported yet!" );
512 }
513 
515 {
517 
519  j > A_bar_.cols(), std::range_error
520  ,"ConstraintsRelaxedStd::get_bnd(j,bnd) : Error, "
521  "j = " << j << " is not in range [1," << A_bar_.cols() << "]" );
522 
523  // See if this is the last constraint we added to the active set.
524  if( j == last_added_j_ && bnd == last_added_bound_type_ ) {
525  return last_added_bound_;
526  }
527 
528  // Lookup the bound! (sparse lookup)
529  size_type j_local = j;
530  const SpVectorSlice::element_type *ele_ptr = NULL;
531  if( j_local <= A_bar_.nd() ) {
532  if(dL_) {
533  switch( bnd ) {
534  case EQUALITY:
535  case LOWER:
536  return dL_->get_ele(j_local);
537  case UPPER:
538  return dU_->get_ele(j_local);
539  default:
541  }
542  }
543  else {
544  return ( bnd == LOWER ? -1.0 : +1.0 ) * inf;
545  }
546  }
547  else if( (j_local -= A_bar_.nd()) <= 1 ) {
548  switch( bnd ) {
549  case EQUALITY:
550  case LOWER:
551  return etaL_;
552  case UPPER:
553  return +inf;
554  default:
556  }
557  }
558  else if( (j_local -= 1) <= A_bar_.m_in() ) {
559  switch( bnd ) {
560  case EQUALITY:
561  case LOWER:
562  return eL_->get_ele(j_local);
563  case UPPER:
564  return eU_->get_ele(j_local);
565  default:
567  }
568  }
569  else if( (j_local -= A_bar_.m_in()) <= A_bar_.m_eq() ) {
570  switch( bnd ) {
571  case EQUALITY:
572  case LOWER:
573  case UPPER:
574  return -A_bar_.f()->get_ele(j_local);
575  default:
577  }
578  }
579  return 0.0; // will never be executed!
580 }
581 
583  size_type last_added_j, value_type last_added_bound
584  ,EBounds last_added_bound_type
585  ) const
586 {
587  last_added_j_ = last_added_j;
588  last_added_bound_ = last_added_bound;
589  last_added_bound_type_ = last_added_bound_type;
590 }
591 
592 // members for ConstraintsRelaxedStd::MatrixConstraints
593 
595  :nd_(0)
596  ,m_in_(0)
597  ,m_eq_(0)
598  ,E_(NULL)
599  ,trans_E_(BLAS_Cpp::no_trans)
600  ,b_(NULL)
601  ,F_(NULL)
602  ,trans_F_(BLAS_Cpp::no_trans)
603  ,f_(NULL)
604  ,space_cols_(Teuchos::null)
605  ,space_rows_(NULL,0)
606 {}
607 
609  const VectorSpace::space_ptr_t &space_d_eta
610  ,const size_type m_in
611  ,const size_type m_eq
612  ,const MatrixOp *E
613  ,BLAS_Cpp::Transp trans_E
614  ,const Vector *b
615  ,const MatrixOp *F
616  ,BLAS_Cpp::Transp trans_F
617  ,const Vector *f
618  ,size_type m_undecomp
619  ,const size_type j_f_undecomp[]
620  )
621 {
622  namespace mmp = MemMngPack;
623  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
624 
625  const size_type nd = space_d_eta->dim() - 1;
626 
627  // Setup P_u
628  const bool test_setup = true; // Todo: Make this an argument!
629  if( m_undecomp > 0 && f->dim() > m_undecomp ) {
630  P_u_row_i_.resize(m_undecomp);
631  P_u_col_j_.resize(m_undecomp);
632  const size_type
633  *j_f_u = j_f_undecomp; // This is sorted by row!
634  row_i_t::iterator
635  row_i_itr = P_u_row_i_.begin();
636  col_j_t::iterator
637  col_j_itr = P_u_col_j_.begin();
638  for( size_type i = 1; i <= m_undecomp; ++i, ++j_f_u, ++row_i_itr, ++col_j_itr ) {
639  *row_i_itr = *j_f_u; // This is sorted in asscending order!
640  *col_j_itr = i;
641  }
642  P_u_.initialize(nd,m_undecomp,m_undecomp,0,0,GPMSIP::BY_ROW_AND_COL
643  ,&P_u_row_i_[0],&P_u_col_j_[0],test_setup);
644  }
645  else if( m_undecomp > 0) { // Must be == f->dim()
646  // Set to identity
647  P_u_.initialize(m_undecomp,m_undecomp,GenPermMatrixSlice::IDENTITY_MATRIX);
648  }
649 
650  // space_cols_
651  space_cols_ = space_d_eta;
652 
653  // space_rows_
654  VectorSpace::space_ptr_t row_spaces[3];
655  int num_row_spaces = 1;
656  row_spaces[0] = space_d_eta;
657  if(m_in)
658  row_spaces[num_row_spaces++] = Teuchos::rcp(
659  trans_E == BLAS_Cpp::no_trans ? &E->space_cols() : &E->space_rows()
660  ,false
661  );
662  if(m_eq) {
663  VectorSpace::space_ptr_t
664  vs = Teuchos::rcp(
665  trans_F == BLAS_Cpp::no_trans ? &F->space_cols() : &F->space_rows()
666  ,false
667  );
668  if(m_undecomp)
669  vs = vs->space(P_u_,BLAS_Cpp::trans);
670  row_spaces[num_row_spaces++] = vs;
671  }
672  space_rows_.initialize( row_spaces, num_row_spaces, space_d_eta->small_vec_spc_fcty() );
673 
674  // Set the rest of the members
675  nd_ = space_d_eta->dim() - 1;
676  m_in_ = m_in;
677  m_eq_ = m_eq;
678  E_ = E;
679  trans_E_ = trans_E;
680  b_ = b;
681  F_ = F;
682  trans_F_ = trans_F;
683  f_ = f;
684 
685 }
686 
687 // Overridden from MatrixOp
688 
690 {
691  return *space_cols_;
692 }
693 
695 {
696  return space_rows_;
697 }
698 
700  const MatrixOp& M
701  )
702 {
703  // ToDo: Finish me
705  return *this;
706 }
707 
708 /* 10/25/00 I don't think we need this function yet!
709 void ConstraintsRelaxedStd::MatrixConstraints::Mp_StPtMtP(
710  DMatrixSlice* C, value_type a
711  ,const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
712  ,BLAS_Cpp::Transp M_trans
713  ,const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
714  )const
715 {
716  using BLAS_Cpp::trans_not;
717  using AbstractLinAlgPack::dot;
718  using AbstractLinAlgPack::Vp_StMtV;
719  using AbstractLinAlgPack::Vp_StPtMtV;
720  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
721 
722  //
723  // A_bar = [ I 0 op(E') op(F') ]
724  // [ 0 1 -b' -f' ]
725  //
726 
727  const size_type
728  E_start = nd() + 1 + 1,
729  F_start = E_start + m_in(),
730  F_end = F_start + m_eq() - 1;
731  const Range1D
732  d_rng = Range1D(1,nd()),
733  E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
734  F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
735 
736  // For this to work (as shown below) we need to have P1 sorted by
737  // row if op(P1) = P1' or sorted by column if op(P1) = P1.
738  // Also, we must have P1 sorted by
739  // row if op(P2) = P2 or sorted by column if op(P2) = P2'
740  // If P1 or P2 are not sorted properly, we will just use the default
741  // implementation of this operation.
742  if( ( P1.ordered_by() == GPMSIP::BY_ROW && P1_trans == BLAS_Cpp::no_trans )
743  || ( P1.ordered_by() == GPMSIP::BY_COL && P1_trans == BLAS_Cpp::trans )
744  || ( P2.ordered_by() == GPMSIP::BY_ROW && P2_trans == BLAS_Cpp::trans )
745  || ( P2.ordered_by() == GPMSIP::BY_COL && P2_trans == BLAS_Cpp::no_trans ) )
746  {
747  // Call the default implementation
748  MatrixOp::Vp_StPtMtV(C,a,P1,P1_trans,M_trans,P2,P2_trans);
749  return;
750  }
751 
752  if( M_trans == BLAS_Cpp::no_trans ) {
753  //
754  // C += a * op(P1) * A_bar * op(P2)
755  //
756  // = a * [ op(P11) op(P12) ] * [ I 0 op(E') op(F') ] * [ op(P21) ]
757  // [ 0 1 -b' -f' ] [ op(P22) ]
758  // [ op(P23) ]
759  // [ op(P24) ]
760  //
761  // C += a*op(P11)*op(P21) + a*op(P21)*op(P22)
762  // + a*op(P11)*op(E')*op(P23) - a*op(P12)*b'*op(P23)
763  // + a*op(P11)*op(F')*op(P24) - a*op(P12)*f'*op(P24)
764  //
765 
766  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement This!
767 
768  }
769  else {
770  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish This!
771  }
772 }
773 */
774 
776  VectorMutable* y, value_type a, BLAS_Cpp::Transp trans_rhs1
777  ,const Vector& x, value_type b
778  ) const
779 {
780  TEUCHOS_TEST_FOR_EXCEPT( !( !F_ || P_u_.cols() == f_->dim() ) ); // ToDo: Add P_u when needed!
781 
782  namespace mmp = MemMngPack;
783  using BLAS_Cpp::trans_not;
785  using LinAlgOpPack::Vt_S;
786  using LinAlgOpPack::Vp_StV;
788 
789  // ToDo: Replace with proper check!
790 // LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),rows(),cols(),trans_rhs1,x.dim());
791 
792  //
793  // A_bar = [ I 0 op(E') op(F') ]
794  // [ 0 1 -b' -f' ]
795  //
796 
797  const size_type
798  E_start = nd() + 1 + 1,
799  F_start = E_start + m_in(),
800  F_end = F_start + m_eq() - 1;
801  const Range1D
802  d_rng = Range1D(1,nd()),
803  E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
804  F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
805 
806  // y = b * y
807  Vt_S( y, b );
808 
809  if( trans_rhs1 == BLAS_Cpp::no_trans ) {
810  //
811  // y += a* A_bar * x
812  //
813  // += a * [ I 0 op(E') op(F') ] * [ x1 ]
814  // [ 0 1 -b' -f' ] [ x2 ]
815  // [ x3 ]
816  // [ x4 ]
817  //
818  // [ y1 ] += [ a * x1 + a * op(E') * x3 + a * op(F') * x4 ]
819  // [ y2 ] [ a * x2 - a * b' * x3 - a * f' * x4 ]
820  //
821  VectorMutable::vec_mut_ptr_t
822  y1 = y->sub_view(d_rng);
823  value_type
824  y2 = y->get_ele(nd()+1);
825  Vector::vec_ptr_t
826  x1 = x.sub_view(d_rng);
827  const value_type
828  x2 = x.get_ele(nd()+1);
829  Vector::vec_ptr_t
830  x3 = m_in() ? x.sub_view(E_rng) : Teuchos::null,
831  x4 = m_eq() ? x.sub_view(F_rng) : Teuchos::null;
832 
833  // [ y1 ] += [ a * x1 + a * op(E') * x3 + a * op(F') * x4 ]
834  Vp_StV( y1.get(), a, *x1 );
835  if( m_in() )
836  Vp_StMtV( y1.get(), a, *E(), trans_not( trans_E() ), *x3 );
837  if( m_eq() )
838  Vp_StMtV( y1.get(), a, *F(), trans_not( trans_F() ), *x4 );
839  // [ y2 ] += [ a * x2 - a * b' * x3 - a * f' * x4 ]
840  y2 += a * x2;
841  if( m_in() )
842  y2 += - a * dot( *this->b(), *x3 );
843  if( m_eq() )
844  y2 += - a * dot( *f(), *x4 );
845  y->set_ele(nd()+1,y2);
846  }
847  else if ( trans_rhs1 == BLAS_Cpp::trans ) {
848  //
849  // y += a* A_bar' * x
850  //
851  // += a * [ I 0 ] * [ x1 ]
852  // [ 0 1 ] [ x2 ]
853  // [ op(E) -b ]
854  // [ op(F) -f ]
855  //
856  // [ y1 ] [ a * x1 ]
857  // [ y2 ] [ + a * x2 ]
858  // [ y3 ] += [ a * op(E) * x1 - a * b * x2 ]
859  // [ y4 ] [ a * op(F) * x1 - a * f * x2 ]
860  //
861  VectorMutable::vec_mut_ptr_t
862  y1 = y->sub_view(d_rng);
863  value_type
864  y2 = y->get_ele(nd()+1);
865  VectorMutable::vec_mut_ptr_t
866  y3 = m_in() ? y->sub_view(E_rng) : Teuchos::null,
867  y4 = m_eq() ? y->sub_view(F_rng) : Teuchos::null;
868  Vector::vec_ptr_t
869  x1 = x.sub_view(d_rng);
870  const value_type
871  x2 = x.get_ele(nd()+1);
872  // y1 += a * x1
873  Vp_StV( y1.get(), a, *x1 );
874  // y2 += a * x2
875  y2 += a * x2;
876  y->set_ele(nd()+1,y2);
877  // y3 += a * op(E) * x1 - (a*x2) * b
878  if( m_in() ) {
879  Vp_StMtV( y3.get(), a, *E(), trans_E(), *x1 );
880  Vp_StV( y3.get(), - a * x2, *this->b() );
881  }
882  // y4 += a * op(F) * x1 - (a*x2) * f
883  if( m_eq() ) {
884  Vp_StMtV( y4.get(), a, *F(), trans_F(), *x1 );
885  Vp_StV( y4.get(), - a * x2, *f() );
886  }
887  }
888  else {
889  TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid trans value
890  }
891 }
892 
894  VectorMutable* y_out, value_type a
895  ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
896  ,BLAS_Cpp::Transp M_trans
897  ,const SpVectorSlice& x, value_type beta
898  ) const
899 {
900  MatrixOp::Vp_StPtMtV(y_out,a,P,P_trans,M_trans,x,beta); // ToDo: Update below code!
901 
902 /*
903 
904  TEUCHOS_TEST_FOR_EXCEPT( !( !F_ || P_u_.cols() == f_->dim() ) ); // ToDo: Add P_u when needed!
905 
906  using BLAS_Cpp::no_trans;
907  using BLAS_Cpp::trans;
908  using BLAS_Cpp::trans_not;
909  using DenseLinAlgPack::dot;
910  using AbstractLinAlgPack::Vp_StMtV;
911  using AbstractLinAlgPack::Vp_StPtMtV;
912  namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
913 
914  AbstractLinAlgPack::VectorDenseMutableEncap y_d(*y_out);
915  DVectorSlice *y = &y_d();
916 
917  DenseLinAlgPack::Vp_MtV_assert_sizes(
918  y->dim(),P.rows(),P.cols(),P_trans
919  ,BLAS_Cpp::rows( rows(), cols(), M_trans) );
920  DenseLinAlgPack::Vp_MtV_assert_sizes(
921  BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
922  ,rows(),cols(),M_trans,x.dim());
923 
924  //
925  // A_bar = [ I 0 op(E') op(F') ]
926  // [ 0 1 -b' -f' ]
927  //
928 
929  const size_type
930  E_start = nd() + 1 + 1,
931  F_start = E_start + m_in(),
932  F_end = F_start + m_eq() - 1;
933  const Range1D
934  d_rng = Range1D(1,nd()),
935  E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
936  F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
937 
938  // For this to work (as shown below) we need to have P sorted by
939  // row if op(P) = P' or sorted by column if op(P) = P. If
940  // P is not sorted properly, we will just use the default
941  // implementation of this operation.
942  if( ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::no_trans )
943  || ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::trans )
944  || ( P.ordered_by() == GPMSIP::UNORDERED ) )
945  {
946  // Call the default implementation
947  //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,beta);
948  TEUCHOS_TEST_FOR_EXCEPT(true);
949  return;
950  }
951 
952  if( M_trans == BLAS_Cpp::no_trans ) {
953  //
954  // y = beta*y + a * op(P) * A_bar * x
955  //
956  // = beta * y
957  //
958  // + a * [op(P1) op(P2) ] * [ I 0 op(E') op(F') ] * [ x1 ]
959  // [ 0 1 -b' -f' ] [ x2 ]
960  // [ x3 ]
961  // [ x4 ]
962  //
963  // y = beta*y + a*op(P1)*x1 + a*op(P1)*op(E')*x3 + a*op(P1)*op(F')*x4
964  // + a*op(P2)*x2 - a*op(P2)*b'*x3 - a*op(P2)*f'*x4
965  //
966  // Where:
967  // op(P1) = op(P)(:,1:nd)
968  // op(P2) = op(P)(:,nd+1:nd+1)
969  //
970 
971  const GenPermMatrixSlice
972  P1 = ( P.is_identity()
973  ? GenPermMatrixSlice( nd(), nd(), GenPermMatrixSlice::IDENTITY_MATRIX )
974  : P.create_submatrix(d_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
975  ),
976  P2 = ( P.is_identity()
977  ? GenPermMatrixSlice(
978  P_trans == no_trans ? nd() : 1
979  , P_trans == no_trans ? 1 : nd()
980  , GenPermMatrixSlice::ZERO_MATRIX )
981  : P.create_submatrix(Range1D(nd()+1,nd()+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
982  );
983 
984  const SpVectorSlice
985  x1 = x(d_rng);
986  const value_type
987  x2 = get_sparse_element(x,nd()+1);
988  const SpVectorSlice
989  x3 = m_in() ? x(E_rng) : SpVectorSlice(NULL,0,0,0),
990  x4 = m_eq() ? x(F_rng) : SpVectorSlice(NULL,0,0,0);
991 
992  // y = beta*y + a*op(P1)*x1
993  Vp_StMtV( y, a, P1, P_trans, x1, beta );
994  // y += a*op(P1)*op(E')*x3
995  if( m_in() && P1.nz() )
996  LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, *E(), trans_not(trans_E()), x3 );
997  // y += a*op(P1)*op(F')*x4
998  if( m_eq() && P1.nz() )
999  LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, *F(), trans_not(trans_F()), x4 );
1000  //
1001  // y += a*op(P2)*x2 - a*op(P2)*b'*x3 - a*op(P2)*f'*x4
1002  // += a * op(P2) * ( x2 + b'*x3 - f'*x4 )
1003  //
1004  // ==>
1005  //
1006  // y(i) += a * ( x2 - b'*x3 - f'*x4 )
1007  //
1008  if( P2.nz() ){
1009  TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
1010  const size_type
1011  i = P_trans == BLAS_Cpp::no_trans
1012  ? P2.begin()->row_i() : P2.begin()->col_j();
1013  value_type
1014  &y_i = (*y)(i) += a * x2;
1015  if(m_in())
1016  y_i += -a * dot(*this->b(),x3);
1017  if(m_eq())
1018  y_i += -a * dot(*this->f(),x4);
1019  }
1020  }
1021  else if ( M_trans == BLAS_Cpp::trans ) {
1022  //
1023  // y = beta*y + a*op(P)*A_bar'*x
1024  //
1025  // = beta*y
1026  //
1027  // + a * [ P1 P2 P3 P4 ] * [ I 0 ] * [ x1 ]
1028  // [ 0 1 ] [ x2 ]
1029  // [ op(E) -b ]
1030  // [ op(F) -f ]
1031  //
1032  // y = beta*y + a*P1*x1 + a*P2*x2 + a*P3*op(E)*x1 - a*P3*b*x2
1033  // + a*P4*op(F)*x1 - a*P4*f*x2
1034  //
1035  // Where:
1036  // P1 = op(P)(:,1:nd)
1037  // P2 = op(P)(:,nd+1:nd+1)
1038  // P3 = op(P)(:,nd+1+1:nd+1+m_in)
1039  // P4 = op(P)(:,nd+1+m_in+1:nd+1+m_in+m_eq)
1040  //
1041 
1042  TEUCHOS_TEST_FOR_EXCEPT( !( !P.is_identity() ) ); // We should never have this!
1043 
1044  size_type off = 0;
1045  const GenPermMatrixSlice
1046  P1 = P.create_submatrix(Range1D(off+1,off+nd())
1047  ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL);
1048  off += nd();
1049  const GenPermMatrixSlice
1050  P2 = P.create_submatrix(Range1D(off+1,off+1)
1051  ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL);
1052  off += 1;
1053  const GenPermMatrixSlice
1054  P3 = m_in()
1055  ? P.create_submatrix(Range1D(off+1,off+m_in())
1056  ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
1057  : GenPermMatrixSlice();
1058  off += m_in();
1059  const GenPermMatrixSlice
1060  P4 = m_eq()
1061  ? P.create_submatrix(Range1D(off+1,off+m_eq())
1062  ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
1063  : GenPermMatrixSlice();
1064 
1065  const SpVectorSlice
1066  x1 = x(d_rng);
1067  const value_type
1068  x2 = get_sparse_element(x,nd()+1);
1069 
1070  // y = beta*y + a*op(P1)*x1
1071  Vp_StMtV( y, a, P1, P_trans, x1, beta );
1072  // y += a*op(P2)*x2
1073  if( P2.nz() ){
1074  TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
1075  (*y)( P_trans == BLAS_Cpp::no_trans
1076  ? P2.begin()->row_i() : P2.begin()->col_j() )
1077  += a * x2;
1078  }
1079  if( m_in() && P3.nz() ) {
1080  // y += a*P3*op(E)*x1
1081  LinAlgOpPack::Vp_StPtMtV( y, a, P3, P_trans, *E(), trans_E(), x1 );
1082  // y += (-a*x2)*P3*b
1083  AbstractLinAlgPack::Vp_StMtV( y, - a * x2, P3, P_trans, *this->b() );
1084  }
1085  if( m_eq() && P4.nz() ) {
1086  // y += a*P4*op(F)*x1
1087  LinAlgOpPack::Vp_StPtMtV( y, a, P4, P_trans, *F(), trans_F(), x1 );
1088  // y += (-a*x2)*P4*f
1089  AbstractLinAlgPack::Vp_StMtV( y, - a * x2, P4, P_trans, *this->f() );
1090  }
1091 
1092  }
1093  else {
1094  TEUCHOS_TEST_FOR_EXCEPT(true); // Invalid trans value
1095  }
1096 */
1097 }
1098 
1099 } // end namespace QPSchurPack
1100 } // end namespace ConstrainedOptPack
AbstractLinAlgPack::size_type size_type
void Vp_StMtV(DVectorSlice *vs_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2, value_type beta=1.0)
vs_lhs = alpha * op(mwo_rhs1) * vs_rhs2 + beta * vs_lhs.
void f()
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
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)
v_lhs = alpha * op(P_rhs1) * op(M_rhs2) * v_rhs3 + beta * v_rhs
void cache_last_added(size_type last_added_j, value_type last_added_bound, EBounds last_added_bound_type) const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return rows of a possible transposed matrix.
void Vp_StPtMtV(VectorMutable *vs_lhs, value_type alpha, const GenPermMatrixSlice &P_rhs1, BLAS_Cpp::Transp P_rhs1_trans, BLAS_Cpp::Transp M_rhs2_trans, const SpVectorSlice &sv_rhs3, value_type beta) const
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
Here the next violated constraint to add to the active set is selected.
Transposed.
Not transposed.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void initialize(const VectorSpace::space_ptr_t &space_d_eta, const size_type m_in, const size_type m_eq, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, size_type m_undecomp, const size_type j_f_undecomp[])
Initialize.
const MatrixOp & A_bar() const
Represents the constraints matrix.
const LAPACK_C_Decl::f_int & M
bool max_inequ_viol(const AbstractLinAlgPack::Vector &v, const AbstractLinAlgPack::Vector &vL, const AbstractLinAlgPack::Vector &vU, AbstractLinAlgPack::size_type *max_viol_i, AbstractLinAlgPack::value_type *max_viol, AbstractLinAlgPack::value_type *v_i, int *bnd_type, AbstractLinAlgPack::value_type *vLU_i)
Compute the maximum violation from a set of inequality constraints vL <= v <= vU. ...
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)
void Vp_StMtV(VectorMutable *vs_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1, const Vector &vs_rhs2, value_type beta) const
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
SparseVectorSlice< SparseElement< index_type, value_type > > SpVectorSlice
void initialize(const VectorSpace::space_ptr_t &space_d_eta, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const Vector *eL, const Vector *eU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, size_type m_undecomp, const size_type j_f_undecomp[], VectorMutable *Ed, bool check_F=true, value_type bounds_tol=1e-10, value_type inequality_tol=1e-10, value_type equality_tol=1e-10)
Initialize constriants.
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
const char inf
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
Transp
TRANS.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
Return columns of a possible transposed matrix.
void V_VmV(VectorMutable *v_lhs, const V1 &V1_rhs1, const V2 &V2_rhs2)
v_lhs = V_rhs1 - V_rhs2.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)