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