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

Generated on Thu Sep 18 12:34:16 2008 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by doxygen 1.3.9.1