MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_ConstraintsRelaxedStd.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 //
00042 // ToDo: 12/29/00: Consider scaling when determining if a
00043 // constraint is violated.  We should consider the size of
00044 // ||d(e[j](x))/d(x)||inf but this is expensive to compute
00045 // given the current interfaces.  We really need to rectify
00046 // this!
00047 //
00048 
00049 #include <assert.h>
00050 
00051 #include <limits>
00052 
00053 #include "ConstrainedOptPack_ConstraintsRelaxedStd.hpp"
00054 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00056 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00057 #include "AbstractLinAlgPack_MatrixOp.hpp"
00058 #include "AbstractLinAlgPack_VectorMutable.hpp"
00059 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00060 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00061 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00062 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00063 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00064 #include "Teuchos_Assert.hpp"
00065 
00066 namespace {
00067 
00068 ConstrainedOptPack::EBounds
00069 convert_bnd_type( int bnd_type )
00070 {
00071   switch(bnd_type) {
00072     case -1:
00073       return ConstrainedOptPack::LOWER;
00074     case 0:
00075       return ConstrainedOptPack::EQUALITY;
00076     case +1:
00077       return ConstrainedOptPack::UPPER;
00078     default:
00079       TEUCHOS_TEST_FOR_EXCEPT(true);
00080   }
00081   return ConstrainedOptPack::LOWER; // Never be called
00082 }
00083 
00084 // Get an element from a sparse vector and return zero if it does not exist
00085 AbstractLinAlgPack::value_type get_sparse_element(
00086   const AbstractLinAlgPack::SpVectorSlice& v
00087   ,AbstractLinAlgPack::size_type i
00088   )
00089 {
00090   const AbstractLinAlgPack::SpVectorSlice::element_type
00091     *ele_ptr = v.lookup_element(i);
00092   return ele_ptr ? ele_ptr->value() : 0.0;
00093 }
00094 
00095 } // end namespace
00096 
00097 namespace ConstrainedOptPack {
00098 namespace QPSchurPack {
00099 
00100 // members for ConstraintsRelaxedStd
00101 
00102 ConstraintsRelaxedStd::ConstraintsRelaxedStd()
00103   :inequality_pick_policy_(ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY)
00104   ,etaL_(0.0)
00105   ,dL_(NULL)
00106   ,dU_(NULL)
00107   ,eL_(NULL)
00108   ,eU_(NULL)
00109   ,Ed_(NULL)
00110   ,last_added_j_(0)
00111   ,last_added_bound_(0.0)
00112   ,last_added_bound_type_(FREE)
00113   ,next_undecomp_f_k_(0)
00114 {}
00115 
00116 void ConstraintsRelaxedStd::initialize(
00117   const VectorSpace::space_ptr_t   &space_d_eta
00118   ,value_type                      etaL
00119   ,const Vector              *dL
00120   ,const Vector              *dU
00121   ,const MatrixOp              *E
00122   ,BLAS_Cpp::Transp                trans_E
00123   ,const Vector              *b
00124   ,const Vector              *eL
00125   ,const Vector              *eU
00126   ,const MatrixOp              *F
00127   ,BLAS_Cpp::Transp                trans_F
00128   ,const Vector              *f
00129   ,size_type                       m_undecomp
00130   ,const size_type                 j_f_undecomp[]
00131   ,VectorMutable             *Ed
00132   ,bool                            check_F
00133   ,value_type                      bounds_tol
00134   ,value_type                      inequality_tol
00135   ,value_type                      equality_tol
00136   )
00137 {
00138   size_type
00139     nd   = space_d_eta->dim() - 1,
00140     m_in = 0,
00141     m_eq = 0;
00142 
00143   TEUCHOS_TEST_FOR_EXCEPT( !(  m_undecomp == (F ? f->dim() : 0)  ) ); // ToDo: support decomposed equalities in future.
00144 
00145   // Validate that the correct sets of constraints are selected
00146   TEUCHOS_TEST_FOR_EXCEPTION(
00147     dL && !dU, std::invalid_argument
00148     ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00149     "If dL!=NULL then dU!=NULL must also be true." );
00150   TEUCHOS_TEST_FOR_EXCEPTION(
00151     E && ( !b || !eL || !eU ), std::invalid_argument
00152     ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00153     "If E!=NULL then b!=NULL, eL!=NULL and eU!=NULL must also be true." );
00154   TEUCHOS_TEST_FOR_EXCEPTION(
00155     F && !f, std::invalid_argument
00156     ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00157     "If F!=NULL then f!=NULL must also be true." );
00158 
00159   // Validate input argument sizes
00160   if(dL) {
00161     const size_type dL_dim = dL->dim(), dU_dim = dU->dim();
00162     TEUCHOS_TEST_FOR_EXCEPTION(
00163       dL_dim != nd, std::invalid_argument
00164       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00165       "dL.dim() != d->dim()." );
00166     TEUCHOS_TEST_FOR_EXCEPTION(
00167       dU_dim != nd, std::invalid_argument 
00168       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00169       "dU.dim() != d->dim()." );
00170   }
00171   if(E) {
00172     const size_type
00173       E_rows = E->rows(),
00174       E_cols = E->cols(),
00175       opE_cols = BLAS_Cpp::cols( E_rows, E_cols, trans_E ),
00176       b_dim = b->dim(),
00177       eL_dim = eL->dim(),
00178       eU_dim = eU->dim(),
00179       Ed_dim = Ed ? Ed->dim() : 0;
00180     m_in = BLAS_Cpp::rows( E_rows, E_cols, trans_E );
00181     TEUCHOS_TEST_FOR_EXCEPTION(
00182       opE_cols != nd, std::invalid_argument
00183       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00184       "op(E).cols() != nd." );
00185     TEUCHOS_TEST_FOR_EXCEPTION(
00186       b_dim != m_in, std::invalid_argument
00187       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00188       "b->dim() != op(E).rows()." );
00189     TEUCHOS_TEST_FOR_EXCEPTION(
00190       eL_dim != m_in, std::invalid_argument
00191       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00192       "eL->dim() != op(E).rows()." );
00193     TEUCHOS_TEST_FOR_EXCEPTION(
00194       eU_dim != m_in, std::invalid_argument
00195       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00196       "eU->dim() != op(E).rows()." );
00197     TEUCHOS_TEST_FOR_EXCEPTION(
00198       Ed && Ed_dim != m_in, std::invalid_argument
00199       ,"ConstraintsRelaxedStd::initialize(...) : Error, "
00200       "Ed->dim() != op(E).rows()." );
00201   }
00202   if(F) {
00203     const size_type
00204       F_rows = F->rows(),
00205       F_cols = F->cols(),
00206       opF_cols = BLAS_Cpp::cols( F_rows, F_cols, trans_F ),
00207       f_dim = f->dim();
00208     m_eq = BLAS_Cpp::rows( F_rows, F_cols, trans_F );
00209     TEUCHOS_TEST_FOR_EXCEPTION(
00210       opF_cols != nd, std::invalid_argument
00211       ,"QPSolverRelaxed::solve_qp(...) : Error, "
00212       "op(F).cols() != nd." );
00213     TEUCHOS_TEST_FOR_EXCEPTION(
00214       f_dim != m_eq, std::invalid_argument
00215       ,"QPSolverRelaxed::solve_qp(...) : Error, "
00216       "f->dim() != op(F).rows()." );
00217   }
00218   
00219   // Initialize other members
00220   A_bar_.initialize(
00221     space_d_eta,m_in,m_eq,E,trans_E,b,F,trans_F,f,m_undecomp,j_f_undecomp);
00222   etaL_       = etaL;
00223   dL_         = dL;
00224   dU_         = dU;
00225   eL_         = eL;
00226   eU_         = eU;
00227   Ed_         = Ed;
00228   check_F_      = check_F;
00229   bounds_tol_     = bounds_tol;
00230   inequality_tol_   = inequality_tol;
00231   equality_tol_   = equality_tol;
00232   last_added_j_   = 0;  // No cached value.
00233   next_undecomp_f_k_  = m_undecomp ? 1 : 0; // Check the first undecomposed equality
00234 
00235 }
00236 
00237 const ConstraintsRelaxedStd::MatrixConstraints&
00238 ConstraintsRelaxedStd::A_bar_relaxed() const
00239 {
00240   return A_bar_;
00241 }
00242 
00243 // Overridden from Constraints
00244 
00245 size_type ConstraintsRelaxedStd::n() const
00246 {
00247   return A_bar_.nd() + 1;
00248 }
00249 
00250 size_type ConstraintsRelaxedStd::m_breve() const
00251 {
00252   return A_bar_.m_in() + A_bar_.m_eq();
00253 }
00254 
00255 const MatrixOp& ConstraintsRelaxedStd::A_bar() const
00256 {
00257   return A_bar_;
00258 }
00259 
00260 void ConstraintsRelaxedStd::pick_violated_policy( EPickPolicy pick_policy )
00261 {
00262   switch(pick_policy) {
00263     case ANY_VIOLATED:
00264       inequality_pick_policy_ = ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY;
00265       break;
00266     case MOST_VIOLATED:
00267       inequality_pick_policy_ = ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY;
00268       break;
00269     default:
00270       TEUCHOS_TEST_FOR_EXCEPT(true);
00271   }
00272 }
00273 
00274 Constraints::EPickPolicy
00275 ConstraintsRelaxedStd::pick_violated_policy() const
00276 {
00277   switch(inequality_pick_policy_) {
00278     case ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY:
00279       return ANY_VIOLATED;
00280     case ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY:
00281       return ANY_VIOLATED;
00282     case ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
00283       return MOST_VIOLATED;
00284     default:
00285       TEUCHOS_TEST_FOR_EXCEPT(true);
00286   }
00287   return ANY_VIOLATED;  // will never be executed
00288 }
00289 
00290 void ConstraintsRelaxedStd::pick_violated(
00291   const DVectorSlice& x_in, size_type* j_viol, value_type* constr_val
00292   ,value_type* viol_bnd_val, value_type* norm_2_constr, EBounds* bnd, bool* can_ignore
00293   ) const
00294 {
00295   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00296   using AbstractLinAlgPack::max_inequ_viol;
00297   using LinAlgOpPack::V_VmV;
00298 
00299   TEUCHOS_TEST_FOR_EXCEPTION(
00300     x_in.dim() != A_bar_.nd()+1, std::length_error
00301     ,"ConstraintsRelaxedStd::pick_violated(...) : Error, "
00302     "x is the wrong length" );
00303 
00304   const size_type
00305     nd = A_bar_.nd();
00306 
00307   // Get a version of x = [ d; eta ] in the correct vector object
00308   VectorSpace::vec_mut_ptr_t
00309     x = A_bar_.space_cols().create_member();
00310   (VectorDenseMutableEncap(*x)()) = x_in;
00311   VectorSpace::vec_mut_ptr_t
00312     d = x->sub_view(1,nd);
00313   const value_type
00314     eta = x->get_ele(nd+1);
00315 
00316   bool Ed_computed = false;
00317 
00318   // //////////////////////////////////////////////
00319   // Check the equality constraints first
00320   if( check_F_ && A_bar_.F() ) {
00321     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Update below code!
00322 /*
00323     // The basic strategy here is to go through all of the equality
00324     // constraints first and add all of the ones that are violated by
00325     // more than the set tolerance.  Those that are not sufficiently
00326     // violated are passed over but they are remembered for later.
00327     // Only when all of the constraints have been gone through once
00328     // will those passed over constraints be considered and then only
00329     // if ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY is selected.
00330     const GenPermMatrixSlice& P_u = A_bar_.P_u();
00331     size_type e_k_mat_row, e_k_mat_col = 1; // Mapping matrix for e(j)
00332     GenPermMatrixSlice e_k_mat; // ToDo: Change to EtaVector
00333     if( next_undecomp_f_k_ <= P_u.nz() ) {
00334       // This is our first pass through the undecomposed equalities.
00335       GenPermMatrixSlice::const_iterator P_u_itr, P_u_end; // only used if P_u is not identity
00336       if( !P_u.is_identity() ) {
00337         P_u_itr = P_u.begin() + (next_undecomp_f_k_ - 1);
00338         P_u_end = P_u.end();
00339       }
00340       for( ; next_undecomp_f_k_ <=  P_u.nz(); ) {
00341         size_type k = 0;
00342         if( !P_u.is_identity() ) {
00343           k = P_u_itr->row_i();
00344           ++P_u_itr;
00345         }
00346         else {
00347           k = next_undecomp_f_k_;
00348         }
00349         ++next_undecomp_f_k_;
00350         // evaluate the constraint e(k)'*[op(F)*d + (1-eta)*f]
00351         value_type
00352           Fd_k = 0.0,
00353           f_k = (*A_bar_.f())(k);
00354         e_k_mat.initialize(
00355           A_bar_.m_eq(),1,1,0,0,GPMSIP::BY_ROW_AND_COL
00356           ,&(e_k_mat_row=k),&e_k_mat_col,false );
00357         DVectorSlice Fd_k_vec(&Fd_k,1);
00358         AbstractLinAlgPack::Vp_StPtMtV(   // ToDo: Use transVtMtV(...) instead!
00359           &Fd_k_vec, 1.0, e_k_mat, BLAS_Cpp::trans
00360           ,*A_bar_.F(), A_bar_.trans_F(), d, 0.0 );
00361         const value_type
00362           err = ::fabs(Fd_k + (1.0 - eta)*f_k) / (1.0 + ::fabs(f_k));
00363         if( err > equality_tol() ) {
00364           *j_viol         = nd + 1 + A_bar_.m_in() + k;
00365           *constr_val     = Fd_k - eta*f_k;
00366           *norm_2_constr  = 1.0; // ToDo: Compute it some how?
00367           *bnd            = EQUALITY;
00368           *viol_bnd_val   = -f_k;
00369           *can_ignore     = false; // Given this careful algorithm this should be false
00370           // cache this
00371           last_added_j_     = *j_viol;
00372           last_added_bound_type_  = *bnd;
00373           last_added_bound_   = *viol_bnd_val;
00374           return;
00375         }
00376         else {
00377           passed_over_equalities_.push_back(k); // remember for later
00378         }
00379       }
00380     }
00381     else if(
00382       inequality_pick_policy() == ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY
00383       && passed_over_equalities_.size() > 0
00384       )
00385     {
00386       // Now look through the list of passed over equalities and see which one
00387       // is violated.  If a passed over constraint is added to the active set
00388       // then it is removed from this list.
00389       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00390     }
00391 */
00392   }
00393 
00394   // /////////////////////////////////////////////
00395   // Find the most violated variable bound.
00396 
00397   size_type       max_bound_viol_j       = 0;
00398   value_type      max_bound_viol         = 0.0;
00399   value_type      max_bound_d_viol       = 0.0;
00400   value_type      max_bound_dLU_viol     = 0.0;
00401   int             max_bound_viol_type    = -2;
00402   if( dL_ && ( dL_->nz() || dU_->nz() ) ) {
00403     // dL <= d <= dU
00404     max_inequ_viol(
00405       *d, *dL_, *dU_
00406       ,&max_bound_viol_j, &max_bound_viol
00407       ,&max_bound_d_viol, &max_bound_viol_type, &max_bound_dLU_viol 
00408       );
00409     if(  max_bound_viol > bounds_tol_ ) {
00410       // Set the return values
00411       *j_viol         = max_bound_viol_j;
00412       *constr_val     = max_bound_d_viol;
00413       *norm_2_constr  = 1.0; // This is correct
00414       *bnd            = convert_bnd_type(max_bound_viol_type);
00415       *viol_bnd_val   = max_bound_dLU_viol;
00416       *can_ignore     = false;
00417     }
00418     else {
00419       max_bound_viol_j = 0; // No variable bounds sufficiently violated.
00420     }
00421   }
00422 
00423   if( ( inequality_pick_policy_ == ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY
00424       ||  inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY )
00425     && max_bound_viol_j
00426     )
00427   {
00428     // A variable bound has been violated so lets just return it!
00429     last_added_j_           = *j_viol;
00430     last_added_bound_type_  = *bnd;
00431     last_added_bound_       = *viol_bnd_val;
00432     return;
00433   }
00434 
00435   // /////////////////////////////////////////////
00436   // Check the general inequalities
00437 
00438   size_type       max_inequality_viol_j        = 0;
00439   value_type      max_inequality_viol          = 0.0;
00440   value_type      max_inequality_e_viol        = 0.0;
00441   value_type      max_inequality_eLU_viol      = 0.0;
00442   int             max_inequality_viol_type     = -2;
00443 
00444   if( inequality_pick_policy_ == ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY ) {
00445     // Find the first general inequality violated by more than
00446     // the defined tolerance.
00447     TEUCHOS_TEST_FOR_EXCEPTION(
00448       true, std::logic_error
00449       ,"ConstraintsRelaxedStd::pick_violated(...) : Error,\n"
00450       "The option ADD_BOUNDS_THEN_FIRST_VIOLATED_INEQUALITY has not been implemented yet\n" );
00451   }
00452   else {
00453     // Find the most violated inequality constraint
00454     if( A_bar_.m_in() &&  ( eL_->nz() || eU_->nz() ) ) {
00455       // e = op(E)*d - b*eta
00456       VectorSpace::vec_mut_ptr_t e = eL_->space().create_member();
00457       LinAlgOpPack::V_MtV( e.get(), *A_bar_.E(), A_bar_.trans_E(), *d );
00458       if(Ed_) {
00459         *Ed_        = *e;
00460         Ed_computed = true;
00461       }
00462       LinAlgOpPack::Vp_StV( e.get(), -eta, *A_bar_.b() );
00463       // eL <= e <= eU
00464       max_inequ_viol(
00465         *e, *eL_, *eU_
00466         ,&max_inequality_viol_j, &max_inequality_viol
00467         ,&max_inequality_e_viol, &max_inequality_viol_type, &max_inequality_eLU_viol 
00468         );
00469       if( max_inequality_viol > max_bound_viol
00470         && max_inequality_viol > inequality_tol_ )
00471       {
00472         *j_viol         = max_inequality_viol_j + nd + 1; // offset into A_bar
00473         *constr_val     = max_inequality_e_viol;
00474         *norm_2_constr  = 1.0;  // This is not correct!
00475         *bnd            = convert_bnd_type(max_inequality_viol_type);
00476         *viol_bnd_val   = max_inequality_eLU_viol;
00477         *can_ignore     = false;
00478       }
00479       else {
00480         max_inequality_viol_j = 0; // No general inequality constraints sufficiently violated.
00481       }
00482     }
00483   }
00484 
00485   if( max_bound_viol_j || max_inequality_viol_j ) {
00486     // One of the constraints has been violated so just return it.
00487     last_added_j_           = *j_viol;
00488     last_added_bound_type_  = *bnd;
00489     last_added_bound_       = *viol_bnd_val;
00490     return;
00491   }
00492 
00493   // If we get here then no constraint was found that violated any of the tolerances.
00494   if( Ed_ && !Ed_computed ) {
00495     // Ed = op(E)*d
00496     LinAlgOpPack::V_MtV( Ed_, *A_bar_.E(), A_bar_.trans_E(), *d );
00497   }
00498   *j_viol     = 0;
00499   *constr_val   = 0.0;
00500   *viol_bnd_val = 0.0;
00501   *norm_2_constr  = 0.0;
00502   *bnd      = FREE;  // Meaningless
00503   *can_ignore   = false; // Meaningless
00504 }
00505 
00506 void ConstraintsRelaxedStd::ignore( size_type j )
00507 {
00508   TEUCHOS_TEST_FOR_EXCEPTION(
00509     true, std::logic_error
00510     ,"ConstraintsRelaxedStd::ignore(...) : Error, "
00511     "This operation is not supported yet!" );
00512 }
00513 
00514 value_type ConstraintsRelaxedStd::get_bnd( size_type j, EBounds bnd ) const
00515 {
00516   const value_type inf = std::numeric_limits<value_type>::max();
00517 
00518   TEUCHOS_TEST_FOR_EXCEPTION(
00519     j > A_bar_.cols(), std::range_error
00520     ,"ConstraintsRelaxedStd::get_bnd(j,bnd) : Error, "
00521     "j = " << j << " is not in range [1," << A_bar_.cols() << "]" );
00522 
00523   // See if this is the last constraint we added to the active set.
00524   if( j == last_added_j_ && bnd == last_added_bound_type_ ) {
00525     return last_added_bound_;
00526   }
00527 
00528   // Lookup the bound! (sparse lookup)
00529   size_type j_local = j;
00530   const SpVectorSlice::element_type *ele_ptr = NULL;
00531   if( j_local <= A_bar_.nd() ) {
00532     if(dL_) {
00533       switch( bnd ) {
00534         case EQUALITY:
00535         case LOWER:
00536           return dL_->get_ele(j_local);
00537         case UPPER:
00538           return dU_->get_ele(j_local);
00539         default:
00540           TEUCHOS_TEST_FOR_EXCEPT(true);
00541       }
00542     }
00543     else {
00544       return ( bnd == LOWER ? -1.0 : +1.0 ) * inf;
00545     }
00546   }
00547   else if( (j_local -= A_bar_.nd()) <= 1 ) {
00548     switch( bnd ) {
00549       case EQUALITY:
00550       case LOWER:
00551         return etaL_;
00552       case UPPER:
00553         return +inf;
00554       default:
00555         TEUCHOS_TEST_FOR_EXCEPT(true);
00556     }
00557   }
00558   else if( (j_local -= 1) <= A_bar_.m_in() ) {
00559     switch( bnd ) {
00560       case EQUALITY:
00561       case LOWER:
00562         return eL_->get_ele(j_local);
00563       case UPPER:
00564         return eU_->get_ele(j_local);
00565       default:
00566         TEUCHOS_TEST_FOR_EXCEPT(true);
00567     }
00568   }
00569   else if( (j_local -= A_bar_.m_in()) <= A_bar_.m_eq() ) {
00570     switch( bnd ) {
00571       case EQUALITY:
00572       case LOWER:
00573       case UPPER:
00574         return -A_bar_.f()->get_ele(j_local);
00575       default:
00576         TEUCHOS_TEST_FOR_EXCEPT(true);
00577     }
00578   }
00579   return 0.0; // will never be executed!
00580 }
00581 
00582 void ConstraintsRelaxedStd::cache_last_added(
00583   size_type last_added_j, value_type last_added_bound
00584   ,EBounds last_added_bound_type
00585   ) const
00586 {
00587   last_added_j_           = last_added_j;
00588   last_added_bound_       = last_added_bound;
00589   last_added_bound_type_  = last_added_bound_type;
00590 }
00591 
00592 // members for ConstraintsRelaxedStd::MatrixConstraints
00593 
00594 ConstraintsRelaxedStd::MatrixConstraints::MatrixConstraints()
00595   :nd_(0)
00596   ,m_in_(0)
00597   ,m_eq_(0)
00598   ,E_(NULL)
00599   ,trans_E_(BLAS_Cpp::no_trans)
00600   ,b_(NULL)
00601   ,F_(NULL)
00602   ,trans_F_(BLAS_Cpp::no_trans)
00603   ,f_(NULL)
00604   ,space_cols_(Teuchos::null)
00605   ,space_rows_(NULL,0)
00606 {}
00607 
00608 void ConstraintsRelaxedStd::MatrixConstraints::initialize(
00609   const VectorSpace::space_ptr_t   &space_d_eta
00610   ,const size_type                 m_in
00611   ,const size_type                 m_eq
00612   ,const MatrixOp                  *E
00613   ,BLAS_Cpp::Transp                trans_E
00614   ,const Vector                    *b
00615   ,const MatrixOp                  *F
00616   ,BLAS_Cpp::Transp                trans_F
00617   ,const Vector                    *f
00618   ,size_type                       m_undecomp
00619   ,const size_type                 j_f_undecomp[]
00620   )
00621 {
00622   namespace mmp = MemMngPack;
00623   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00624 
00625   const size_type nd = space_d_eta->dim() - 1;
00626 
00627   // Setup P_u
00628   const bool test_setup = true; // Todo: Make this an argument!
00629   if( m_undecomp > 0 && f->dim() > m_undecomp ) {
00630     P_u_row_i_.resize(m_undecomp);
00631     P_u_col_j_.resize(m_undecomp);
00632     const size_type
00633       *j_f_u = j_f_undecomp;  // This is sorted by row!
00634     row_i_t::iterator
00635       row_i_itr = P_u_row_i_.begin();
00636     col_j_t::iterator
00637       col_j_itr = P_u_col_j_.begin();
00638     for( size_type i = 1; i <= m_undecomp; ++i, ++j_f_u, ++row_i_itr, ++col_j_itr ) {
00639       *row_i_itr = *j_f_u; // This is sorted in asscending order!
00640       *col_j_itr = i;
00641     }
00642     P_u_.initialize(nd,m_undecomp,m_undecomp,0,0,GPMSIP::BY_ROW_AND_COL
00643       ,&P_u_row_i_[0],&P_u_col_j_[0],test_setup);
00644   }
00645   else if( m_undecomp > 0) { // Must be == f->dim()
00646     // Set to identity
00647     P_u_.initialize(m_undecomp,m_undecomp,GenPermMatrixSlice::IDENTITY_MATRIX);
00648   }
00649 
00650   // space_cols_
00651   space_cols_ = space_d_eta;
00652 
00653   // space_rows_
00654   VectorSpace::space_ptr_t  row_spaces[3];
00655   int num_row_spaces = 1;
00656   row_spaces[0] = space_d_eta;
00657   if(m_in)
00658     row_spaces[num_row_spaces++] = Teuchos::rcp(
00659       trans_E == BLAS_Cpp::no_trans ? &E->space_cols() : &E->space_rows()
00660       ,false
00661       );
00662   if(m_eq) {
00663     VectorSpace::space_ptr_t
00664       vs = Teuchos::rcp(
00665         trans_F == BLAS_Cpp::no_trans ? &F->space_cols() : &F->space_rows()
00666         ,false
00667         );
00668     if(m_undecomp)
00669       vs = vs->space(P_u_,BLAS_Cpp::trans);
00670     row_spaces[num_row_spaces++] = vs;
00671   }
00672   space_rows_.initialize( row_spaces, num_row_spaces, space_d_eta->small_vec_spc_fcty() );
00673 
00674   // Set the rest of the members
00675   nd_       = space_d_eta->dim() - 1;
00676   m_in_     = m_in;
00677   m_eq_     = m_eq;
00678   E_        = E;
00679   trans_E_  = trans_E;
00680   b_        = b;
00681   F_        = F;
00682   trans_F_  = trans_F;
00683   f_        = f;
00684 
00685 }
00686 
00687 // Overridden from MatrixOp
00688 
00689 const VectorSpace& ConstraintsRelaxedStd::MatrixConstraints::space_cols() const
00690 {
00691   return *space_cols_;
00692 }
00693 
00694 const VectorSpace& ConstraintsRelaxedStd::MatrixConstraints::space_rows() const
00695 {
00696   return space_rows_;
00697 }
00698 
00699 MatrixOp& ConstraintsRelaxedStd::MatrixConstraints::operator=(
00700   const MatrixOp& M
00701   )
00702 {
00703   // ToDo: Finish me
00704   TEUCHOS_TEST_FOR_EXCEPT(true);
00705   return *this;
00706 }
00707 
00708 /* 10/25/00 I don't think we need this function yet!
00709 void ConstraintsRelaxedStd::MatrixConstraints::Mp_StPtMtP(
00710   DMatrixSlice* C, value_type a
00711   ,const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00712   ,BLAS_Cpp::Transp M_trans
00713   ,const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00714   )const
00715 {
00716   using BLAS_Cpp::trans_not;
00717   using AbstractLinAlgPack::dot;
00718   using AbstractLinAlgPack::Vp_StMtV;
00719   using AbstractLinAlgPack::Vp_StPtMtV;
00720   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00721 
00722   //  
00723   //  A_bar = [  I   0  op(E')   op(F')  ]
00724   //          [  0   1   -b'      -f'    ]
00725   //
00726 
00727   const size_type
00728     E_start = nd() + 1 + 1,
00729     F_start = E_start + m_in(),
00730     F_end = F_start + m_eq() - 1;
00731   const Range1D
00732     d_rng = Range1D(1,nd()),
00733     E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
00734     F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
00735 
00736   // For this to work (as shown below) we need to have P1 sorted by
00737   // row if op(P1) = P1' or sorted by column if op(P1) = P1.
00738   // Also, we must have P1 sorted by
00739   // row if op(P2) = P2 or sorted by column if op(P2) = P2'
00740   // If P1 or P2 are not sorted properly, we will just use the default
00741   // implementation of this operation.
00742   if(   ( P1.ordered_by() == GPMSIP::BY_ROW && P1_trans == BLAS_Cpp::no_trans )
00743       ||  ( P1.ordered_by() == GPMSIP::BY_COL && P1_trans == BLAS_Cpp::trans )
00744       ||  ( P2.ordered_by() == GPMSIP::BY_ROW && P2_trans == BLAS_Cpp::trans )
00745       ||  ( P2.ordered_by() == GPMSIP::BY_COL && P2_trans == BLAS_Cpp::no_trans ) )
00746   {
00747     // Call the default implementation
00748     MatrixOp::Vp_StPtMtV(C,a,P1,P1_trans,M_trans,P2,P2_trans);
00749     return;
00750   }
00751 
00752   if( M_trans == BLAS_Cpp::no_trans ) {
00753     //
00754     // C += a * op(P1) * A_bar * op(P2)
00755     //
00756     //    = a * [ op(P11)  op(P12) ] * [ I  0  op(E')  op(F') ] * [ op(P21) ]
00757     //                                 [ 0  1    -b'    -f'   ]   [ op(P22) ]
00758     //                                                            [ op(P23) ]
00759     //                                                            [ op(P24) ]
00760     //
00761     // C +=   a*op(P11)*op(P21) + a*op(P21)*op(P22)
00762     //      + a*op(P11)*op(E')*op(P23) - a*op(P12)*b'*op(P23)
00763     //      + a*op(P11)*op(F')*op(P24) - a*op(P12)*f'*op(P24)
00764     //      
00765 
00766     TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: Implement This!
00767 
00768   }
00769   else {
00770     TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: Finish This!
00771   }
00772 }
00773 */
00774 
00775 void ConstraintsRelaxedStd::MatrixConstraints::Vp_StMtV(
00776   VectorMutable* y, value_type a, BLAS_Cpp::Transp trans_rhs1
00777   ,const Vector& x, value_type b
00778   ) const
00779 {
00780   TEUCHOS_TEST_FOR_EXCEPT( !(  !F_ || P_u_.cols() == f_->dim()  ) ); // ToDo: Add P_u when needed!
00781 
00782   namespace mmp = MemMngPack;
00783   using BLAS_Cpp::trans_not;
00784   using AbstractLinAlgPack::dot;
00785   using LinAlgOpPack::Vt_S;
00786   using LinAlgOpPack::Vp_StV;
00787   using LinAlgOpPack::Vp_StMtV;
00788 
00789   // ToDo: Replace with proper check!
00790 //  LinAlgOpPack::Vp_MtV_assert_sizes(y->dim(),rows(),cols(),trans_rhs1,x.dim());
00791 
00792   //  
00793   //  A_bar = [  I   0  op(E')   op(F')  ]
00794   //          [  0   1   -b'      -f'    ]
00795   //
00796 
00797   const size_type
00798     E_start = nd() + 1 + 1,
00799     F_start = E_start + m_in(),
00800     F_end = F_start + m_eq() - 1;
00801   const Range1D
00802     d_rng = Range1D(1,nd()),
00803     E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
00804     F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
00805 
00806   // y = b * y
00807   Vt_S( y, b );
00808   
00809   if( trans_rhs1 == BLAS_Cpp::no_trans ) {
00810     //
00811     // y += a* A_bar * x
00812     // 
00813     //   += a * [ I  0  op(E')  op(F') ] * [ x1 ]
00814     //          [ 0  1   -b'     -f'   ]   [ x2 ]
00815     //                                     [ x3 ]
00816     //                                     [ x4 ]
00817     //
00818     // [ y1 ]  += [ a * x1 + a * op(E') * x3 + a * op(F') * x4 ]
00819     // [ y2 ]     [ a * x2 - a * b' * x3     - a * f' * x4     ]
00820     //
00821     VectorMutable::vec_mut_ptr_t
00822       y1 = y->sub_view(d_rng);
00823     value_type
00824       y2 = y->get_ele(nd()+1);
00825     Vector::vec_ptr_t
00826       x1 = x.sub_view(d_rng);
00827     const value_type
00828       x2 = x.get_ele(nd()+1);
00829     Vector::vec_ptr_t
00830       x3 = m_in() ? x.sub_view(E_rng) : Teuchos::null,
00831       x4 = m_eq() ? x.sub_view(F_rng) : Teuchos::null;
00832     
00833     // [ y1 ]  += [ a * x1 + a * op(E') * x3 + a * op(F') * x4 ]
00834     Vp_StV( y1.get(), a, *x1 );
00835     if( m_in() )
00836       Vp_StMtV( y1.get(), a, *E(), trans_not( trans_E() ), *x3 );
00837     if( m_eq() )
00838       Vp_StMtV( y1.get(), a, *F(), trans_not( trans_F() ), *x4 );
00839     // [ y2 ]  += [ a * x2 - a * b' * x3     - a * f' * x4     ]
00840     y2 += a * x2;
00841     if( m_in() )
00842       y2 += - a * dot( *this->b(), *x3 );
00843     if( m_eq() )
00844       y2 += - a * dot( *f(), *x4 );
00845     y->set_ele(nd()+1,y2);
00846   }
00847   else if ( trans_rhs1 == BLAS_Cpp::trans ) {
00848     //
00849     // y += a* A_bar' * x
00850     // 
00851     //   += a * [ I       0 ] * [ x1 ]
00852     //          [ 0       1 ]   [ x2 ]
00853     //          [ op(E)  -b ]
00854     //          [ op(F)  -f ]
00855     //
00856     // [ y1 ]    [ a * x1                        ]
00857     // [ y2 ]    [                + a * x2       ]
00858     // [ y3 ] += [ a * op(E) * x1 - a * b * x2   ]
00859     // [ y4 ]    [ a * op(F) * x1 - a * f * x2   ]
00860     //
00861     VectorMutable::vec_mut_ptr_t
00862       y1 = y->sub_view(d_rng);
00863     value_type
00864       y2 = y->get_ele(nd()+1);
00865     VectorMutable::vec_mut_ptr_t
00866       y3 = m_in() ? y->sub_view(E_rng) : Teuchos::null,
00867       y4 = m_eq() ? y->sub_view(F_rng) : Teuchos::null;
00868     Vector::vec_ptr_t
00869       x1 = x.sub_view(d_rng);
00870     const value_type
00871       x2 = x.get_ele(nd()+1);
00872     // y1 += a * x1
00873     Vp_StV( y1.get(), a, *x1 );
00874     // y2 += a * x2
00875     y2 += a * x2;
00876     y->set_ele(nd()+1,y2);
00877     // y3 += a * op(E) * x1 - (a*x2) * b
00878     if( m_in() ) {
00879       Vp_StMtV( y3.get(), a, *E(), trans_E(), *x1 );
00880       Vp_StV( y3.get(), - a * x2, *this->b() );
00881     }
00882     // y4 += a * op(F) * x1 - (a*x2) * f
00883     if( m_eq() ) {
00884       Vp_StMtV( y4.get(), a, *F(), trans_F(), *x1 );
00885       Vp_StV( y4.get(), - a * x2, *f() );
00886     }
00887   }
00888   else {
00889     TEUCHOS_TEST_FOR_EXCEPT(true);  // Invalid trans value
00890   }
00891 }
00892 
00893 void ConstraintsRelaxedStd::MatrixConstraints::Vp_StPtMtV(
00894   VectorMutable* y_out, value_type a
00895   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00896   ,BLAS_Cpp::Transp M_trans
00897   ,const SpVectorSlice& x, value_type beta
00898   ) const
00899 {
00900   MatrixOp::Vp_StPtMtV(y_out,a,P,P_trans,M_trans,x,beta); // ToDo: Update below code!
00901   
00902 /*
00903 
00904   TEUCHOS_TEST_FOR_EXCEPT( !(  !F_ || P_u_.cols() == f_->dim()  ) ); // ToDo: Add P_u when needed!
00905 
00906   using BLAS_Cpp::no_trans;
00907   using BLAS_Cpp::trans;
00908   using BLAS_Cpp::trans_not;
00909    using DenseLinAlgPack::dot;
00910   using AbstractLinAlgPack::Vp_StMtV;
00911   using AbstractLinAlgPack::Vp_StPtMtV;
00912   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00913 
00914   AbstractLinAlgPack::VectorDenseMutableEncap y_d(*y_out);
00915   DVectorSlice *y = &y_d();
00916 
00917   DenseLinAlgPack::Vp_MtV_assert_sizes(
00918     y->dim(),P.rows(),P.cols(),P_trans
00919     ,BLAS_Cpp::rows( rows(), cols(), M_trans) );
00920   DenseLinAlgPack::Vp_MtV_assert_sizes(
00921     BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
00922     ,rows(),cols(),M_trans,x.dim());
00923 
00924   //  
00925   //  A_bar = [  I   0  op(E')   op(F')  ]
00926   //          [  0   1   -b'      -f'    ]
00927   //
00928 
00929   const size_type
00930     E_start = nd() + 1 + 1,
00931     F_start = E_start + m_in(),
00932     F_end = F_start + m_eq() - 1;
00933   const Range1D
00934     d_rng = Range1D(1,nd()),
00935     E_rng = m_in() ? Range1D(E_start,F_start-1) : Range1D(),
00936     F_rng = m_eq() ? Range1D(F_start,F_end) : Range1D();
00937 
00938   // For this to work (as shown below) we need to have P sorted by
00939   // row if op(P) = P' or sorted by column if op(P) = P.  If
00940   // P is not sorted properly, we will just use the default
00941   // implementation of this operation.
00942   if(   ( P.ordered_by() == GPMSIP::BY_ROW && P_trans == BLAS_Cpp::no_trans )
00943       ||  ( P.ordered_by() == GPMSIP::BY_COL && P_trans == BLAS_Cpp::trans )
00944     ||  ( P.ordered_by() == GPMSIP::UNORDERED ) )
00945   {
00946     // Call the default implementation
00947     //MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,beta);
00948     TEUCHOS_TEST_FOR_EXCEPT(true);
00949     return;
00950   }
00951 
00952   if( M_trans == BLAS_Cpp::no_trans ) {
00953     //
00954     // y = beta*y + a * op(P) * A_bar * x
00955     //
00956     //   = beta * y
00957     //   
00958     //    + a * [op(P1)  op(P2) ] * [ I  0  op(E')  op(F') ] * [ x1 ]
00959     //                              [ 0  1   -b'     -f'   ]   [ x2 ]
00960     //                                                         [ x3 ]
00961     //                                                         [ x4 ]
00962     //
00963     // y = beta*y + a*op(P1)*x1 + a*op(P1)*op(E')*x3 + a*op(P1)*op(F')*x4
00964     //     + a*op(P2)*x2 - a*op(P2)*b'*x3 - a*op(P2)*f'*x4
00965     //
00966     // Where:
00967     //   op(P1) = op(P)(:,1:nd)
00968     //   op(P2) = op(P)(:,nd+1:nd+1)
00969     //
00970 
00971     const GenPermMatrixSlice
00972       P1 = ( P.is_identity() 
00973            ? GenPermMatrixSlice( nd(), nd(), GenPermMatrixSlice::IDENTITY_MATRIX )
00974            : P.create_submatrix(d_rng,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00975         ),
00976       P2 = ( P.is_identity()
00977            ? GenPermMatrixSlice(
00978              P_trans == no_trans ? nd() : 1
00979              , P_trans == no_trans ? 1 : nd()
00980              , GenPermMatrixSlice::ZERO_MATRIX )
00981            : P.create_submatrix(Range1D(nd()+1,nd()+1),P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
00982         );
00983 
00984     const SpVectorSlice
00985       x1 = x(d_rng);
00986     const value_type
00987       x2 = get_sparse_element(x,nd()+1);
00988     const SpVectorSlice
00989       x3 = m_in() ? x(E_rng) : SpVectorSlice(NULL,0,0,0),
00990       x4 = m_eq() ? x(F_rng) : SpVectorSlice(NULL,0,0,0);
00991 
00992     // y = beta*y + a*op(P1)*x1
00993     Vp_StMtV( y, a, P1, P_trans, x1, beta );
00994     // y += a*op(P1)*op(E')*x3
00995     if( m_in() && P1.nz() )
00996       LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, *E(), trans_not(trans_E()), x3 );
00997     // y += a*op(P1)*op(F')*x4
00998     if( m_eq() && P1.nz() )
00999       LinAlgOpPack::Vp_StPtMtV( y, a, P1, P_trans, *F(), trans_not(trans_F()), x4 );
01000     //
01001     // y += a*op(P2)*x2 - a*op(P2)*b'*x3 - a*op(P2)*f'*x4
01002     //   += a * op(P2) * ( x2 + b'*x3 - f'*x4 )
01003     //   
01004     //   ==>
01005     //   
01006     // y(i) +=  a * ( x2 - b'*x3 - f'*x4 )
01007     //   
01008     if( P2.nz() ){
01009       TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
01010       const size_type
01011         i = P_trans == BLAS_Cpp::no_trans
01012           ? P2.begin()->row_i() : P2.begin()->col_j();
01013       value_type
01014         &y_i = (*y)(i) += a * x2;
01015       if(m_in())
01016         y_i += -a * dot(*this->b(),x3);
01017       if(m_eq())
01018         y_i += -a * dot(*this->f(),x4);
01019     }
01020   }
01021   else if ( M_trans == BLAS_Cpp::trans ) {
01022     //
01023     // y = beta*y + a*op(P)*A_bar'*x
01024     // 
01025     //   = beta*y
01026     //   
01027     //    + a * [ P1  P2  P3  P4 ] * [ I       0 ] * [ x1 ]
01028     //                               [ 0       1 ]   [ x2 ]
01029     //                               [ op(E)  -b ]
01030     //                               [ op(F)  -f ]
01031     //
01032     // y = beta*y + a*P1*x1 + a*P2*x2 + a*P3*op(E)*x1 - a*P3*b*x2
01033     //     + a*P4*op(F)*x1 - a*P4*f*x2
01034     //
01035     // Where:
01036     //   P1 = op(P)(:,1:nd)
01037     //   P2 = op(P)(:,nd+1:nd+1)
01038     //   P3 = op(P)(:,nd+1+1:nd+1+m_in)
01039     //   P4 = op(P)(:,nd+1+m_in+1:nd+1+m_in+m_eq)
01040     //
01041 
01042     TEUCHOS_TEST_FOR_EXCEPT( !(  !P.is_identity()  ) ); // We should never have this!
01043 
01044     size_type off = 0;
01045     const GenPermMatrixSlice
01046       P1 = P.create_submatrix(Range1D(off+1,off+nd())
01047                   ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL);
01048     off += nd();
01049     const GenPermMatrixSlice
01050       P2 = P.create_submatrix(Range1D(off+1,off+1)
01051                   ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL);
01052     off += 1;
01053     const GenPermMatrixSlice
01054       P3 = m_in()
01055         ? P.create_submatrix(Range1D(off+1,off+m_in())
01056                    ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
01057         : GenPermMatrixSlice();
01058     off += m_in();
01059     const GenPermMatrixSlice
01060       P4 = m_eq()
01061         ? P.create_submatrix(Range1D(off+1,off+m_eq())
01062                    ,P_trans==trans?GPMSIP::BY_ROW:GPMSIP::BY_COL)
01063         : GenPermMatrixSlice();
01064 
01065     const SpVectorSlice
01066       x1 = x(d_rng);
01067     const value_type
01068       x2 = get_sparse_element(x,nd()+1);
01069 
01070     // y = beta*y + a*op(P1)*x1
01071     Vp_StMtV( y, a, P1, P_trans, x1, beta );
01072     // y += a*op(P2)*x2
01073     if( P2.nz() ){
01074       TEUCHOS_TEST_FOR_EXCEPT( !( P2.nz() == 1 ) );
01075       (*y)( P_trans == BLAS_Cpp::no_trans
01076           ? P2.begin()->row_i() : P2.begin()->col_j() )
01077         += a * x2;
01078     }
01079     if( m_in() && P3.nz() ) {
01080       // y += a*P3*op(E)*x1
01081       LinAlgOpPack::Vp_StPtMtV( y, a, P3, P_trans, *E(), trans_E(), x1 );
01082       // y += (-a*x2)*P3*b
01083       AbstractLinAlgPack::Vp_StMtV( y, - a * x2, P3, P_trans, *this->b() );
01084     }
01085     if( m_eq() && P4.nz() ) {
01086       // y += a*P4*op(F)*x1
01087       LinAlgOpPack::Vp_StPtMtV( y, a, P4, P_trans, *F(), trans_F(), x1 );
01088       // y += (-a*x2)*P4*f
01089       AbstractLinAlgPack::Vp_StMtV( y, - a * x2, P4, P_trans, *this->f() );
01090     }
01091     
01092   }
01093   else {
01094     TEUCHOS_TEST_FOR_EXCEPT(true);  // Invalid trans value
01095   }
01096 */
01097 }
01098 
01099 } // end namespace QPSchurPack 
01100 } // end namespace ConstrainedOptPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines