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