ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSolverRelaxedQPOPT.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 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_QPOPT
00043 
00044 #include <assert.h>
00045 
00046 #include <algorithm>
00047 #include <ostream>
00048 #include <iomanip>
00049 
00050 #include "ConstrainedOptPack_QPSolverRelaxedQPOPT.hpp"
00051 #include "ConstrainedOptPack_QPOPT_CppDecl.hpp"
00052 #include "AbstractLinAlgPack_MatrixOp.hpp"
00053 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00054 #include "AbstractLinAlgPack_EtaVector.hpp"
00055 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00056 #include "AbstractLinAlgPack_VectorMutableDense.hpp"
00057 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00058 #include "DenseLinAlgPack_DMatrixOut.hpp"
00059 #include "DenseLinAlgPack_DVectorOut.hpp"
00060 
00061 // //////////////////////////////////////////////////////////
00062 // Local implementation functions.
00063 
00064 namespace {
00065 
00066 typedef FortranTypes::f_int       f_int;
00067 typedef FortranTypes::f_dbl_prec  f_dbl_prec;
00068 typedef FortranTypes::f_logical   f_logical;
00069 
00070 // Compute:
00071 //
00072 // HESS * x = [ G  0 ] * [ X(1,N-1) ] = [ G * X(1,N-1) ]
00073 //            [ 0  M ]   [   X(N)   ]   [ M * X(N)     ]
00074 //
00075 // The matrix vector product is implemented through the MatrixOp interface.
00076 //
00077 inline
00078 void qphess_server_relax( const f_int& N, const f_int& LDH
00079   , const f_int& JTHCOL, const f_dbl_prec* H, const f_dbl_prec* X, f_dbl_prec* HX
00080   , f_int* IW, const f_int& LENIW, f_dbl_prec* W, const f_int& LENW )
00081 {
00082   using DenseLinAlgPack::DVectorSlice;
00083   using AbstractLinAlgPack::SpVector;
00084   using LinAlgOpPack::V_MtV;
00085   using ConstrainedOptPack::QPSolverRelaxedQPOPT;
00086 
00087   // Here we have used some casting to pass on information about the qp solver
00088   // that called QPSOL.
00089   const QPSolverRelaxedQPOPT* qp_solver = reinterpret_cast<const QPSolverRelaxedQPOPT*>(H);
00090   const AbstractLinAlgPack::MatrixOp &G = *qp_solver->G();
00091 
00092   const DVectorSlice x(const_cast<f_dbl_prec*>(X),N);   // Just a view!
00093   DVectorSlice hx(HX,N);                                // Just a view!
00094 
00095   if( JTHCOL == 0 ) {
00096     // We are performing a dense mat-vec
00097     // hx(1,N-1) = G * x(1,N-1)
00098     AbstractLinAlgPack::VectorMutableDense x_n(x(1,N-1),Teuchos::null);
00099     AbstractLinAlgPack::VectorMutableDense hx_n(hx(1,N-1),Teuchos::null);
00100     V_MtV( &hx_n, G, BLAS_Cpp::no_trans, x_n );
00101     // hx(N) = bigM * x(N)
00102     hx(N) = qp_solver->use_as_bigM() * x(N);
00103   }
00104   else {
00105     // we are extracting the JTHCOL column of G so use sparse x
00106     if(JTHCOL == N) {
00107       // 0
00108       hx(1,N-1) = 0.0;
00109       // bigM
00110       hx(N) = qp_solver->use_as_bigM();
00111     }
00112     else {
00113       // G(:,JTHCOL)
00114       AbstractLinAlgPack::EtaVector e_j(JTHCOL,N-1);
00115       AbstractLinAlgPack::VectorMutableDense hx_n(hx(1,N-1),Teuchos::null);
00116       V_MtV( &hx_n, G, BLAS_Cpp::no_trans, e_j() );
00117       // 0
00118       hx(N) = 0.0;
00119     }
00120   }
00121 }
00122 
00123 } // end namespace
00124 
00125 // ///////////////////////////////////////////////////////////////////////////
00126 // Fortran declarations.
00127 
00128 extern "C" {
00129 
00130 // These are declarations for the subroutines that preform the communication
00131 // between C++ and Fortran.  There is no use in putting them in a
00132 // namespace since the namespace name will not be used by the linker since
00133 // we are using extern "C".
00134 
00135 //
00136 FORTRAN_FUNC_DECL_UL_( void, QPHESS_SERVER_RELAX2, qphess_server_relax2 ) ( const f_int& N, const f_int& LDH
00137   , const f_int& JTHCOL, const f_dbl_prec* H, const f_dbl_prec* X, f_dbl_prec* HX
00138   , f_int* IW, const f_int& LENIW, f_dbl_prec* W, const f_int& LENW )
00139 {
00140   qphess_server_relax( N, LDH, JTHCOL, H, X, HX, IW, LENIW, W, LENW );
00141 }
00142 
00143 } // end extern "C"
00144 
00145 namespace LinAlgOpPack {
00146   using AbstractLinAlgPack::Mp_StM;
00147   using AbstractLinAlgPack::Vp_StMtV;
00148 }
00149 
00150 // ///////////////////////////////////////
00151 // QPSolverRelaxedQPOPT
00152 
00153 namespace ConstrainedOptPack {
00154 
00155 QPSolverRelaxedQPOPT::QPSolverRelaxedQPOPT(
00156   value_type        max_qp_iter_frac
00157   )
00158   :max_qp_iter_frac_(max_qp_iter_frac)
00159   ,ITMAX_(100)
00160   ,FEATOL_(1.0e-10)
00161   ,LDH_(1)
00162 {}
00163 
00164 QPSolverRelaxedQPOPT::~QPSolverRelaxedQPOPT()
00165 {
00166   release_memory();
00167 }
00168 
00169 // Overridden from QPSolverRelaxed
00170 
00171 void QPSolverRelaxedQPOPT::release_memory()
00172 {
00173   // ToDo: Resize all arrays to zero!
00174   QPSolverRelaxedQPOPTSOL::release_memory();
00175 }
00176 
00177 // Overridden protected members
00178 
00179 QPSolverRelaxedQPOPT::f_int QPSolverRelaxedQPOPT::liwork(f_int N, f_int NCLIN) const
00180 { return 2* N + 3; }
00181 
00182 QPSolverRelaxedQPOPT::f_int QPSolverRelaxedQPOPT::lrwork(f_int N, f_int NCLIN) const
00183 { return NCLIN > 0 ? 2*N*N + 8*N + 5*NCLIN : N*N + 8 *N; }
00184 
00185 QPSolverRelaxedQPOPTSOL::EInform QPSolverRelaxedQPOPT::call_qp_solver(bool warm_start)
00186 {
00187   
00188   // Set the rest of the QPOPT inputs that could not be set in the constructor.
00189 
00190   BIGBND_ = this->infinite_bound();
00191   ITMAX_ = max_qp_iter_frac() * N_;
00192   LDA_  = ( A_.rows() > 0 ? A_.rows() : 1 );
00193   H_    = reinterpret_cast<f_dbl_prec*>(this);  // used to implement QPHESS
00194   AX_.resize( NCLIN_ > 0 ? NCLIN_ : 1 );
00195 
00196   // Set option parameters
00197   {
00198     namespace ns = QPOPT_CppDecl;
00199     namespace ft = FortranTypes;
00200     ns::reset_defaults();
00201     ns::set_logical_option( ns::WARM_START            , warm_start ? ft::F_FALSE : ft::F_TRUE   );
00202     ns::set_real_option(    ns::FEASIBILITY_TOLERANCE , FEATOL_                                 );
00203     ns::set_real_option(    ns::INFINITE_BOUND_SIZE   , BIGBND_                                 );
00204     ns::set_int_option(     ns::ITERATION_LIMIT       , ITMAX_                                  );
00205     ns::set_int_option(     ns::PRINT_FILE            , 0                                       );
00206     ns::set_int_option(     ns::SUMMARY_FILE          , 0                                       );
00207     ns::set_int_option(     ns::PRINT_LEVEL           , 0                                       );
00208     ns::set_int_option(     ns::PROBLEM_TYPE          , ns::QP2                                 );
00209   }
00210 
00211   QPOPT_CppDecl::qpopt(
00212     N_, NCLIN_, LDA_, LDH_, NCLIN_ ? &A_(1,1) : NULL, &BL_(1), &BU_(1)
00213     , &CVEC_(1), H_
00214     , FORTRAN_NAME_UL_(QPHESS_SERVER_RELAX2,qphess_server_relax2)
00215     , &ISTATE_[0], &X_(1), INFORM_, ITER_, OBJ_, &AX_(1)
00216     , &CLAMDA_(1), &IWORK_[0], LIWORK_, &WORK_[0], LWORK_ );
00217 
00218   EInform return_inform;
00219   typedef QPSolverRelaxedQPOPTSOL bc;
00220   switch(INFORM_) {
00221     case STRONG_LOCAL_MIN:
00222       return_inform = bc::STRONG_LOCAL_MIN;
00223       break;
00224     case WEAK_LOCAL_MIN:
00225       return_inform = bc::WEAK_LOCAL_MIN;
00226       break;
00227     case UNBOUNDED:
00228       TEUCHOS_TEST_FOR_EXCEPTION(
00229         true, Unbounded
00230         ,"QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
00231         " QPOPT returned that the QP is unbounded!" );
00232     case INFEASIBLE:
00233       TEUCHOS_TEST_FOR_EXCEPTION(
00234         true, Infeasible
00235         ,"QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
00236         " QPOPT returned that the QP is infeasible!" );
00237     case ITMAX_EXCEEDED:
00238       return_inform = bc::MAX_ITER_EXCEEDED;
00239       break;
00240     case MAX_DOF_TOO_SMALL:
00241       TEUCHOS_TEST_FOR_EXCEPTION(
00242         true, std::runtime_error,
00243         "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
00244         " QPOPT says that the max dof is too small" );
00245     case INVALID_INPUT:
00246       TEUCHOS_TEST_FOR_EXCEPTION(
00247         true, InvalidInput,
00248         "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
00249         " QPOPT returned that the input is invalid" );
00250     case PROB_TYPE_NOT_REGOG:
00251       TEUCHOS_TEST_FOR_EXCEPTION(
00252         true, std::logic_error,
00253         "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
00254         " QPOPT says that the problem type is not recognized" );
00255       break;
00256     default:
00257       TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen
00258   }
00259   
00260   return return_inform;
00261 }
00262 
00263 } // end namespace ConstrainedOptPack
00264 
00265 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_QPOPT
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends