ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSolverRelaxedQPKWIK.cpp
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 
00043 #include "Moocho_ConfigDefs.hpp"
00044 
00045 
00046 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
00047 
00048 
00049 #include <assert.h>
00050 
00051 #include <vector>
00052 
00053 #include "ConstrainedOptPack_QPSolverRelaxedQPKWIK.hpp"
00054 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00055 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00056 #include "AbstractLinAlgPack_EtaVector.hpp"
00057 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00058 #include "AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
00059 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00060 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00061 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00062 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00063 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00064 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00065 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00066 #include "Teuchos_dyn_cast.hpp"
00067 #include "Teuchos_Assert.hpp"
00068 
00069 namespace QPKWIKNEW_CppDecl {
00070 
00071 // Declarations that will link to the fortran object file.
00072 // These may change for different platforms
00073 
00074 using FortranTypes::f_int;      // INTEGER
00075 using FortranTypes::f_real;     // REAL
00076 using FortranTypes::f_dbl_prec;   // DOUBLE PRECISION
00077 using FortranTypes::f_logical;    // LOGICAL
00078 
00079 // ////////////////////////////////////////////////////
00080 // Declarations to link with Fortran QPKWIK procedures
00081 
00082 namespace Fortran {
00083 extern "C" {
00084 
00085 FORTRAN_FUNC_DECL_UL(void,QPKWIKNEW,qpkwiknew) (
00086   const f_int& N, const f_int& M1, const f_int& M2, const f_int& M3
00087   ,const f_dbl_prec GRAD[], f_dbl_prec UINV[], const f_int& LDUINV
00088   ,const f_int IBND[], const f_dbl_prec BL[], const f_dbl_prec BU[]
00089   ,const f_dbl_prec A[], const f_int& LDA, const f_dbl_prec YPY[]
00090   ,const f_int& IYPY, const f_int& WARM, f_dbl_prec NUMPARAM[], const f_int& MAX_ITER
00091   ,f_dbl_prec X[], f_int* NACTSTORE, f_int IACTSTORE[], f_int* INF
00092   ,f_int* NACT, f_int IACT[], f_dbl_prec UR[], f_dbl_prec* EXTRA
00093   ,f_int* ITER, f_int* NUM_ADDS, f_int* NUM_DROPS
00094   ,f_int ISTATE[], const f_int& LRW, f_dbl_prec RW[]
00095   );
00096 
00097 FORTRAN_FUNC_DECL_UL_(f_int,QPKWIKNEW_LISTATE,qpkwiknew_listate) (
00098   const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3);
00099 
00100 FORTRAN_FUNC_DECL_UL_(f_int,QPKWIKNEW_LRW,qpkwiknew_lrw) (
00101   const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3);
00102 
00103 } // end extern "C"
00104 } // end namespace Fortran
00105 
00106 // //////////////////////////////////
00107 // QPKWIK interface functions
00108 
00109 // Solve a QP using QPKWIK.
00110 //
00111 // See the Fortran file for documentation.  C++ programs should use this interface.
00112 inline
00113 void qpkwiknew ( 
00114   const f_int& n, const f_int& m1, const f_int& m2, const f_int& m3
00115   ,const f_dbl_prec grad[], f_dbl_prec uinv[], const f_int& lduinv
00116   ,const f_int ibnd[], const f_dbl_prec bl[], const f_dbl_prec bu[]
00117   ,const f_dbl_prec a[], const f_int& lda, const f_dbl_prec ypy[]
00118   ,const f_int& iypy, const f_int& warm, f_dbl_prec numparam[], const f_int& max_iter
00119   ,f_dbl_prec x[], f_int* nactstore, f_int iactstore[], f_int* inf
00120   ,f_int* nact, f_int iact[], f_dbl_prec ur[], f_dbl_prec* extra
00121   ,f_int* iter, f_int* num_adds, f_int* num_drops
00122   ,f_int istate[], const f_int& lrw, f_dbl_prec rw[]
00123   )
00124 {
00125   Fortran::FORTRAN_FUNC_CALL_UL(QPKWIKNEW,qpkwiknew) (
00126     n, m1, m2, m3, grad, uinv, lduinv
00127     , ibnd, bl, bu, a, lda, ypy, iypy, warm, numparam, max_iter, x, nactstore
00128     , iactstore, inf, nact, iact, ur, extra, iter, num_adds, num_drops, istate
00129     , lrw, rw
00130     );
00131 }
00132 
00133 // Get the length of the integer state variables
00134 inline
00135 f_int qpkwiknew_listate(const f_int& n, const f_int& m1, const f_int& m2
00136             , const f_int& m3)
00137 {
00138   return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LISTATE,qpkwiknew_listate) (n, m1, m2, m3);
00139 }
00140 
00141 // Get the length of the real (double precision) workspace
00142 inline
00143 f_int qpkwiknew_lrw(const f_int& n, const f_int& m1, const f_int& m2
00144           , const f_int& m3)
00145 {
00146   return Fortran::FORTRAN_FUNC_CALL_UL_(QPKWIKNEW_LRW,qpkwiknew_lrw) (n, m1, m2, m3);
00147 }
00148 
00149 } // end namespace QPKWIKNEW_CppDecl
00150 
00151 // /////////////////////////////////////
00152 // Local helpers
00153 
00154 namespace {
00155 
00156 template< class T >
00157 inline
00158 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00159 
00160 using FortranTypes::f_int;
00161 typedef DenseLinAlgPack::value_type value_type;
00162 
00163 enum EConstraintType { NU_L, NU_U, GAMA_L, GAMA_U, LAMBDA, RELAXATION };
00164 char constraint_type_name[6][15] = { "NU_L", "NU_U", "GAMA_L", "GAMA_U", "LAMBDA", "RELAXATION" };
00165 
00166 EConstraintType constraint_type( const f_int m1, const f_int m2, const f_int m3, const f_int j )
00167 {
00168   if     (1 <= j       && j <= m1      ) return NU_L;
00169   else if(m1+1 <= j    && j <= m1+m2     ) return GAMA_L;
00170   else if(m1+m2+1 <= j   && j <= 2*m1+m2   ) return NU_U;
00171   else if(2*m1+m2+1 <= j   && j <= 2*m1+2*m2   ) return GAMA_U;
00172   else if(2*m1+2*m2+1 <= j && j <= 2*m1+2*m2+m3) return LAMBDA;
00173   else if( j == 2*m1+2*m2+m3 + 1         ) return RELAXATION;
00174   TEUCHOS_TEST_FOR_EXCEPT(true);
00175   return NU_L;  // should never be exectuted
00176 }
00177 
00178 f_int constraint_index( const f_int m1, const f_int m2, const f_int m3, const f_int ibnd[]
00179             , const EConstraintType type, const f_int j )
00180 {
00181   switch(type) {
00182     case NU_L   : return ibnd[j-1];
00183     case GAMA_L   : return j-m1;
00184     case NU_U   : return ibnd[j-m1-m2-1];
00185     case GAMA_U   : return j-2*m1-m2;
00186     case LAMBDA   : return j-2*m1-2*m2;
00187     case RELAXATION : return 0;
00188   }
00189   TEUCHOS_TEST_FOR_EXCEPT(true);
00190   return 0; // should never be exectuted
00191 }
00192 
00193 } // end namespace
00194 
00195 // ///////////////////////////////////////
00196 // Members for QPSolverRelaxedQPKWIK
00197 
00198 namespace ConstrainedOptPack {
00199 
00200 QPSolverRelaxedQPKWIK::QPSolverRelaxedQPKWIK(
00201   value_type        max_qp_iter_frac
00202   ,value_type       infinite_bound
00203   )
00204   :max_qp_iter_frac_(max_qp_iter_frac)
00205   ,infinite_bound_(infinite_bound)
00206   ,N_(0)
00207   ,M1_(0)
00208   ,M2_(0)
00209   ,M3_(0)
00210 {
00211   NUMPARAM_[0] = 1e-10; // SMALL
00212   NUMPARAM_[1] = 1e-20; // VSMALL
00213   NUMPARAM_[2] = 1e+20; // VLARGE
00214 }
00215 
00216 QPSolverRelaxedQPKWIK::~QPSolverRelaxedQPKWIK()
00217 {
00218   this->release_memory();
00219 }
00220 
00221 // Overridden from QPSolverRelaxed
00222 
00223 QPSolverStats
00224 QPSolverRelaxedQPKWIK::get_qp_stats() const
00225 {
00226   return qp_stats_;
00227 }
00228 
00229 void QPSolverRelaxedQPKWIK::release_memory()
00230 {
00231   // Todo: resize to zero all the workspace!
00232 }
00233 
00234 QPSolverStats::ESolutionType
00235 QPSolverRelaxedQPKWIK::imp_solve_qp(
00236   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00237   ,const Vector& g, const MatrixSymOp& G
00238   ,value_type etaL
00239   ,const Vector* dL, const Vector* dU
00240   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00241   ,const Vector* eL, const Vector* eU
00242   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00243   ,value_type* obj_d
00244   ,value_type* eta, VectorMutable* d
00245   ,VectorMutable* nu
00246   ,VectorMutable* mu, VectorMutable* Ed
00247   ,VectorMutable* lambda, VectorMutable* Fd
00248   )
00249 {
00250   using Teuchos::dyn_cast;
00251   using DenseLinAlgPack::nonconst_tri_ele;
00252   using LinAlgOpPack::dot;
00253   using LinAlgOpPack::V_StV;
00254   using LinAlgOpPack::assign;
00255   using LinAlgOpPack::V_StV;
00256   using LinAlgOpPack::V_MtV;
00257   using AbstractLinAlgPack::EtaVector;
00258   using AbstractLinAlgPack::transVtMtV;
00259   using AbstractLinAlgPack::num_bounded;
00260   using ConstrainedOptPack::MatrixExtractInvCholFactor;
00261 
00262   // /////////////////////////
00263   // Map to QPKWIK input
00264 
00265   // Validate that rHL is of the proper type.
00266   const MatrixExtractInvCholFactor &cG
00267     = dyn_cast<const MatrixExtractInvCholFactor>(G);
00268 
00269   // Determine the number of sparse bounds on variables and inequalities.
00270   // By default set for the dense case
00271   const value_type inf = this->infinite_bound();
00272   const size_type
00273     nd              = d->dim(),
00274     m_in            = E  ? b->dim()                  : 0,
00275     m_eq            = F  ? f->dim()                  : 0,
00276     nvarbounds      = dL ? num_bounded(*dL,*dU,inf)  : 0,
00277     ninequbounds    = E  ? num_bounded(*eL,*eU,inf)  : 0,
00278     nequalities     = F  ? f->dim()                  : 0;
00279 
00280   // Determine if this is a QP with a structure different from the
00281   // one just solved.
00282   
00283   const bool same_qp_struct = (  N_ == nd && M1_ == nvarbounds && M2_ == ninequbounds && M3_ == nequalities );
00284 
00286   // Set the input parameters to be sent to QPKWIKNEW
00287 
00288   // N
00289   N_ = nd;
00290 
00291   // M1
00292   M1_ = nvarbounds;
00293 
00294   // M2
00295   M2_ = ninequbounds;
00296 
00297   // M3
00298   M3_ = nequalities;
00299 
00300   // GRAD
00301   GRAD_ = VectorDenseEncap(g)();
00302 
00303   // UINV_AUG
00304   //
00305   // UINV_AUG = [ sqrt(bigM)  0  ]
00306   //            [ 0           L' ]
00307   //
00308   UINV_AUG_.resize(N_+1,N_+1);
00309   cG.extract_inv_chol( &nonconst_tri_ele( UINV_AUG_(2,N_+1,2,N_+1), BLAS_Cpp::upper ) );
00310   UINV_AUG_(1,1) = 1.0 / ::sqrt( NUMPARAM_[2] );
00311   UINV_AUG_.col(1)(2,N_+1) = 0.0;
00312   UINV_AUG_.row(1)(2,N_+1) = 0.0;
00313 
00314   // LDUINV_AUG
00315   LDUINV_AUG_ = UINV_AUG_.rows();
00316 
00317   // IBND, BL , BU, A, LDA, YPY
00318 
00319   IBND_INV_.resize( nd + m_in);
00320   std::fill( IBND_INV_.begin(), IBND_INV_.end(), 0 ); // Initialize the zero
00321   IBND_.resize( my_max( 1, M1_ + M2_ ) );
00322   BL_.resize( my_max( 1, M1_ + M2_ ) );
00323   BU_.resize( my_max( 1, M1_ + M2_ + M3_ ) );
00324   LDA_ = my_max( 1, M2_ + M3_ );
00325   A_.resize( LDA_, (  M2_ + M3_ > 0 ? N_ : 1 ) );
00326   YPY_.resize( my_max( 1, M1_ + M2_ ) );
00327   if(M1_)
00328     YPY_(1,M1_) = 0.0; // Must be for this QP interface
00329 
00330   // Initialize variable bound constraints
00331   if( dL ) {
00332     VectorDenseEncap dL_de(*dL);
00333     VectorDenseEncap dU_de(*dU);
00334     // read iterators
00335     AbstractLinAlgPack::sparse_bounds_itr
00336       dLU_itr( dL_de().begin(), dL_de().end()
00337           ,dU_de().begin(), dU_de().end()
00338           ,inf );
00339     // written iterators
00340     IBND_t::iterator
00341       IBND_itr = IBND_.begin(),
00342       IBND_end = IBND_.begin() + M1_;
00343     DVector::iterator
00344       BL_itr = BL_.begin(),
00345       BU_itr = BU_.begin(),
00346       YPY_itr = YPY_.begin();
00347     // Loop
00348     for( size_type ibnd_i = 1; IBND_itr != IBND_end; ++ibnd_i, ++dLU_itr ) {
00349       IBND_INV_[dLU_itr.index()-1] = ibnd_i;
00350       *IBND_itr++ = dLU_itr.index();
00351       *BL_itr++ = dLU_itr.lbound();
00352       *BU_itr++ = dLU_itr.ubound();
00353       *YPY_itr++  = 0.0; // Must be zero with this QP interface
00354     }
00355   }
00356 
00357   // Initialize inequality constraints
00358   
00359   if(M2_) {
00360     VectorDenseEncap eL_de(*eL);
00361     VectorDenseEncap eU_de(*eU);
00362     VectorDenseEncap b_de(*b);
00363     AbstractLinAlgPack::sparse_bounds_itr
00364       eLU_itr( eL_de().begin(), eL_de().end()
00365           ,eU_de().begin(), eU_de().end()
00366           ,inf );
00367     if( M2_ < m_in ) {
00368       // Initialize BL, BU, YPY and A for sparse bounds on general inequalities
00369       // written iterators
00370       DVector::iterator
00371         BL_itr    = BL_.begin() + M1_,
00372         BU_itr    = BU_.begin() + M1_,
00373         YPY_itr   = YPY_.begin() + M1_;
00374       IBND_t::iterator
00375         ibnds_itr = IBND_.begin() + M1_;
00376       // loop
00377       for(size_type i = 1; i <= M2_; ++i, ++eLU_itr, ++ibnds_itr ) {
00378         TEUCHOS_TEST_FOR_EXCEPT( !( !eLU_itr.at_end() ) );
00379         const size_type k      = eLU_itr.index();
00380         *BL_itr++              = eLU_itr.lbound();
00381         *BU_itr++              = eLU_itr.ubound();
00382         *YPY_itr++             = b_de()(k);
00383         *ibnds_itr             = k;  // Only for my record, not used by QPKWIK
00384         IBND_INV_[nd+k-1]      = M1_ + i;
00385         // Add the corresponding row of op(E) to A
00386         // y == A.row(i)'
00387         // y' = e_k' * op(E) => y = op(E')*e_k
00388         DVectorSlice y = A_.row(i);
00389         EtaVector e_k(k,eL_de().dim());
00390         V_MtV( &y( 1, N_ ), *E, BLAS_Cpp::trans_not(trans_E), e_k() ); // op(E')*e_k
00391       }
00392     }
00393     else {
00394       //
00395       // Initialize BL, BU, YPY and A for dense bounds on general inequalities
00396       //
00397       // Initialize BL(M1+1:M1+M2), BU(M1+1:M1+M2)
00398       // and IBND(M1+1:M1+M2) = identity (only for my record, not used by QPKWIK)
00399       DVector::iterator
00400         BL_itr    = BL_.begin() + M1_,
00401         BU_itr    = BU_.begin() + M1_;
00402       IBND_t::iterator
00403         ibnds_itr = IBND_.begin() + M1_;
00404       for(size_type i = 1; i <= m_in; ++i ) {
00405         if( !eLU_itr.at_end() && eLU_itr.index() == i ) {
00406           *BL_itr++ = eLU_itr.lbound();
00407           *BU_itr++ = eLU_itr.ubound();
00408           ++eLU_itr;
00409         }
00410         else {
00411           *BL_itr++ = -inf;
00412           *BU_itr++ = +inf;
00413         }
00414         *ibnds_itr++     = i;
00415         IBND_INV_[nd+i-1]= M1_ + i;
00416       }
00417       // A(1:M2,1:N) = op(E)
00418       assign( &A_(1,M2_,1,N_), *E, trans_E );
00419       // YPY
00420       YPY_(M1_+1,M1_+M2_) = b_de();
00421     }
00422   }
00423 
00424   // Initialize equalities
00425 
00426   if(M3_) {
00427     V_StV( &BU_( M1_ + M2_ + 1, M1_ + M2_ + M3_ ), -1.0, VectorDenseEncap(*f)() );
00428     assign( &A_( M2_ + 1, M2_ + M3_, 1, N_ ), *F, trans_F );
00429   }
00430 
00431   // IYPY
00432   IYPY_ = 1; // ???
00433 
00434   // WARM
00435   WARM_ = 0; // Cold start by default
00436 
00437   // MAX_ITER
00438   MAX_ITER_ = static_cast<f_int>(max_qp_iter_frac() * N_);
00439 
00440   // INF
00441   INF_ = ( same_qp_struct ? 1 : 0 );
00442   
00443   // Initilize output, internal state and workspace quantities.
00444   if(!same_qp_struct) {
00445     X_.resize(N_);
00446     NACTSTORE_ = 0;
00447     IACTSTORE_.resize(N_+1);
00448     IACT_.resize(N_+1);
00449     UR_.resize(N_+1);
00450     ISTATE_.resize( QPKWIKNEW_CppDecl::qpkwiknew_listate(N_,M1_,M2_,M3_) );
00451     LRW_ = QPKWIKNEW_CppDecl::qpkwiknew_lrw(N_,M1_,M2_,M3_);
00452     RW_.resize(LRW_);
00453   }
00454 
00455   // /////////////////////////////////////////////
00456   // Setup a warm start form the input arguments
00457   //
00458   // Interestingly enough, QPKWIK sorts all of the
00459   // constraints according to scaled multiplier values
00460   // and mixes equality with inequality constriants.
00461   // It seems to me that you should start with equality
00462   // constraints first.
00463 
00464   WARM_      = 0;
00465   NACTSTORE_ = 0;
00466 
00467   if( m_eq ) {
00468     // Add equality constraints first since we know these will
00469     // be part of the active set.
00470     for( size_type j = 1; j <= m_eq; ++j ) {
00471       IACTSTORE_[NACTSTORE_] = 2*M1_ + 2*M2_ + j;
00472       ++NACTSTORE_;
00473     }
00474   }
00475   if( ( nu && nu->nz() ) || ( mu && mu->nz() ) ) {
00476     // Add inequality constraints
00477     const size_type
00478       nu_nz = nu ? nu->nz() : 0,
00479       mu_nz = mu ? mu->nz() : 0;
00480     // Combine all the multipliers for the bound and general inequality
00481     // constraints and sort them from the largest to the smallest.  Hopefully
00482     // the constraints with the larger multiplier values will not be dropped
00483     // from the active set.
00484     SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
00485     typedef SpVector::element_type ele_t;
00486     if(nu && nu_nz) {
00487       VectorDenseEncap nu_de(*nu);
00488       DVectorSlice::const_iterator
00489         nu_itr = nu_de().begin(),
00490         nu_end = nu_de().end();
00491       index_type i = 1;
00492       while( nu_itr != nu_end ) {
00493         for( ; *nu_itr == 0.0; ++nu_itr, ++i );
00494         gamma.add_element(ele_t(i,*nu_itr));
00495         ++nu_itr; ++i;
00496       }
00497     }
00498     if(mu && mu_nz) {
00499       VectorDenseEncap mu_de(*mu);
00500       DVectorSlice::const_iterator
00501         mu_itr = mu_de().begin(),
00502         mu_end = mu_de().end();
00503       index_type i = 1;
00504       while( mu_itr != mu_end ) {
00505         for( ; *mu_itr == 0.0; ++mu_itr, ++i );
00506         gamma.add_element(ele_t(i+nd,*mu_itr));
00507         ++mu_itr; ++i;
00508       }
00509     }
00510     std::sort( gamma.begin(), gamma.end()
00511       , AbstractLinAlgPack::SortByDescendingAbsValue() );
00512     // Now add the inequality constraints in decreasing order
00513     const SpVector::difference_type o = gamma.offset();
00514     for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
00515       const size_type  j   = itr->index() + o;
00516       const value_type val = itr->value();
00517       if( j <= nd ) { // Variable bound
00518         const size_type ibnd_i = IBND_INV_[j-1];
00519         TEUCHOS_TEST_FOR_EXCEPT( !( ibnd_i ) );
00520         IACTSTORE_[NACTSTORE_]
00521           = (val < 0.0
00522              ? ibnd_i               // lower bound (see IACT(*))
00523              : M1_ + M2_ + ibnd_i   // upper bound (see IACT(*))
00524             );
00525         ++NACTSTORE_;
00526       }
00527       else if( j <= nd + m_in ) { // General inequality constraint
00528         const size_type ibnd_i = IBND_INV_[j-1]; // offset into M1_ + ibnd_j
00529         TEUCHOS_TEST_FOR_EXCEPT( !( ibnd_i ) );
00530         IACTSTORE_[NACTSTORE_]
00531           = (val < 0.0
00532              ? ibnd_i               // lower bound (see IACT(*))
00533              : M1_ + M2_ + ibnd_i   // upper bound (see IACT(*))
00534             );
00535         ++NACTSTORE_;
00536       }
00537     }
00538   }
00539   if( NACTSTORE_ > 0 )
00540     WARM_ = 1;
00541 
00542   // /////////////////////////
00543   // Call QPKWIK
00544 
00545   if( out && olevel > PRINT_NONE ) {
00546     *out
00547       << "\nCalling QPKWIK to solve QP problem ...\n";
00548   }
00549 
00550   QPKWIKNEW_CppDecl::qpkwiknew(
00551     N_, M1_, M2_, M3_, &GRAD_(1), &UINV_AUG_(1,1), LDUINV_AUG_, &IBND_[0]
00552     ,&BL_(1), &BU_(1), &A_(1,1), LDA_, &YPY_(1), IYPY_, WARM_, NUMPARAM_, MAX_ITER_, &X_(1)
00553     ,&NACTSTORE_, &IACTSTORE_[0], &INF_, &NACT_, &IACT_[0], &UR_[0], &EXTRA_
00554     ,&ITER_, &NUM_ADDS_, &NUM_DROPS_, &ISTATE_[0], LRW_, &RW_[0]
00555     );
00556 
00557   // ////////////////////////
00558   // Map from QPKWIK output
00559 
00560   // eta
00561   *eta = EXTRA_;
00562   // d
00563   (VectorDenseMutableEncap(*d))() = X_();
00564   // nu (simple variable bounds) and mu (general inequalities)
00565   if(nu) *nu = 0.0;
00566   if(mu) *mu = 0.0;
00567   // ToDo: Create VectorDenseMutableEncap views for faster access!
00568   {for(size_type i = 1; i <= NACT_; ++i) {
00569     size_type j = IACT_[i-1];
00570     EConstraintType type = constraint_type(M1_,M2_,M3_,j);
00571     FortranTypes::f_int idc = constraint_index(M1_,M2_,M3_,&IBND_[0],type,j);
00572     switch(type) {
00573       case NU_L:
00574         nu->set_ele( idc , -UR_(i) );
00575         break;
00576       case GAMA_L:
00577         mu->set_ele( IBND_[ M1_ + idc - 1 ], -UR_(i) );
00578         break;
00579       case NU_U:
00580         nu->set_ele( idc, UR_(i)) ;
00581         break;
00582       case GAMA_U:
00583         mu->set_ele( IBND_[ M1_ + idc - 1 ], UR_(i) );
00584         break;
00585       case LAMBDA:
00586         lambda->set_ele( idc, UR_(i) );
00587         break;
00588     }
00589   }}
00590   // obj_d (This could be updated within QPKWIK in the future)
00591   if(obj_d) {
00592     // obj_d = g'*d + 1/2 * d' * G * g
00593     *obj_d = dot(g,*d) + 0.5 * transVtMtV(*d,G,BLAS_Cpp::no_trans,*d);
00594   }
00595   // Ed (This could be updated within QPKWIK in the future)
00596   if(Ed) {
00597     V_MtV( Ed, *E, trans_E, *d );
00598   }
00599   // Fd (This could be updated within QPKWIK in the future)
00600   if(Fd) {
00601     V_MtV( Fd, *F, trans_F, *d );
00602   }
00603   // Set the QP statistics
00604   QPSolverStats::ESolutionType solution_type;
00605   if( INF_ >= 0 ) {
00606     solution_type = QPSolverStats::OPTIMAL_SOLUTION;
00607   }
00608   else if( INF_ == -1 ) { // Infeasible constraints
00609     TEUCHOS_TEST_FOR_EXCEPTION(
00610       true, QPSolverRelaxed::Infeasible
00611       ,"QPSolverRelaxedQPKWIK::solve_qp(...) : Error, QP is infeasible" );
00612   }
00613   else if( INF_ == -2 ) { // LRW too small
00614     TEUCHOS_TEST_FOR_EXCEPT( !( INF_ != -2 ) );  // Local programming error?
00615   }
00616   else if( INF_ == -3 ) { // Max iterations exceeded
00617     solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
00618   }
00619   else {
00620     TEUCHOS_TEST_FOR_EXCEPT(true); // Unknown return value!
00621   }
00622   qp_stats_.set_stats(
00623     solution_type, QPSolverStats::CONVEX
00624     ,ITER_, NUM_ADDS_, NUM_DROPS_
00625     ,WARM_==1, *eta > 0.0 );
00626 
00627   return qp_stats_.solution_type();
00628 }
00629 
00630 
00631 } // end namespace ConstrainedOptPack
00632 
00633 
00634 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_QPKWIK
00635 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends