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

Generated on Wed May 12 21:52:29 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7