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

Generated on Tue Jul 13 09:30:51 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7