ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPSolverRelaxedQPOPT.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_QPOPT
43 
44 #include <assert.h>
45 
46 #include <algorithm>
47 #include <ostream>
48 #include <iomanip>
49 
50 #include "ConstrainedOptPack_QPSolverRelaxedQPOPT.hpp"
51 #include "ConstrainedOptPack_QPOPT_CppDecl.hpp"
52 #include "AbstractLinAlgPack_MatrixOp.hpp"
53 #include "AbstractLinAlgPack_SpVectorClass.hpp"
54 #include "AbstractLinAlgPack_EtaVector.hpp"
55 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
56 #include "AbstractLinAlgPack_VectorMutableDense.hpp"
57 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
58 #include "DenseLinAlgPack_DMatrixOut.hpp"
59 #include "DenseLinAlgPack_DVectorOut.hpp"
60 
61 // //////////////////////////////////////////////////////////
62 // Local implementation functions.
63 
64 namespace {
65 
66 typedef FortranTypes::f_int f_int;
67 typedef FortranTypes::f_dbl_prec f_dbl_prec;
68 typedef FortranTypes::f_logical f_logical;
69 
70 // Compute:
71 //
72 // HESS * x = [ G 0 ] * [ X(1,N-1) ] = [ G * X(1,N-1) ]
73 // [ 0 M ] [ X(N) ] [ M * X(N) ]
74 //
75 // The matrix vector product is implemented through the MatrixOp interface.
76 //
77 inline
78 void qphess_server_relax( const f_int& N, const f_int& LDH
79  , const f_int& JTHCOL, const f_dbl_prec* H, const f_dbl_prec* X, f_dbl_prec* HX
80  , f_int* IW, const f_int& LENIW, f_dbl_prec* W, const f_int& LENW )
81 {
82  using DenseLinAlgPack::DVectorSlice;
83  using AbstractLinAlgPack::SpVector;
84  using LinAlgOpPack::V_MtV;
85  using ConstrainedOptPack::QPSolverRelaxedQPOPT;
86 
87  // Here we have used some casting to pass on information about the qp solver
88  // that called QPSOL.
89  const QPSolverRelaxedQPOPT* qp_solver = reinterpret_cast<const QPSolverRelaxedQPOPT*>(H);
90  const AbstractLinAlgPack::MatrixOp &G = *qp_solver->G();
91 
92  const DVectorSlice x(const_cast<f_dbl_prec*>(X),N); // Just a view!
93  DVectorSlice hx(HX,N); // Just a view!
94 
95  if( JTHCOL == 0 ) {
96  // We are performing a dense mat-vec
97  // hx(1,N-1) = G * x(1,N-1)
98  AbstractLinAlgPack::VectorMutableDense x_n(x(1,N-1),Teuchos::null);
99  AbstractLinAlgPack::VectorMutableDense hx_n(hx(1,N-1),Teuchos::null);
100  V_MtV( &hx_n, G, BLAS_Cpp::no_trans, x_n );
101  // hx(N) = bigM * x(N)
102  hx(N) = qp_solver->use_as_bigM() * x(N);
103  }
104  else {
105  // we are extracting the JTHCOL column of G so use sparse x
106  if(JTHCOL == N) {
107  // 0
108  hx(1,N-1) = 0.0;
109  // bigM
110  hx(N) = qp_solver->use_as_bigM();
111  }
112  else {
113  // G(:,JTHCOL)
114  AbstractLinAlgPack::EtaVector e_j(JTHCOL,N-1);
115  AbstractLinAlgPack::VectorMutableDense hx_n(hx(1,N-1),Teuchos::null);
116  V_MtV( &hx_n, G, BLAS_Cpp::no_trans, e_j() );
117  // 0
118  hx(N) = 0.0;
119  }
120  }
121 }
122 
123 } // end namespace
124 
125 // ///////////////////////////////////////////////////////////////////////////
126 // Fortran declarations.
127 
128 extern "C" {
129 
130 // These are declarations for the subroutines that preform the communication
131 // between C++ and Fortran. There is no use in putting them in a
132 // namespace since the namespace name will not be used by the linker since
133 // we are using extern "C".
134 
135 //
136 FORTRAN_FUNC_DECL_UL_( void, QPHESS_SERVER_RELAX2, qphess_server_relax2 ) ( const f_int& N, const f_int& LDH
137  , const f_int& JTHCOL, const f_dbl_prec* H, const f_dbl_prec* X, f_dbl_prec* HX
138  , f_int* IW, const f_int& LENIW, f_dbl_prec* W, const f_int& LENW )
139 {
140  qphess_server_relax( N, LDH, JTHCOL, H, X, HX, IW, LENIW, W, LENW );
141 }
142 
143 } // end extern "C"
144 
145 namespace LinAlgOpPack {
148 }
149 
150 // ///////////////////////////////////////
151 // QPSolverRelaxedQPOPT
152 
153 namespace ConstrainedOptPack {
154 
155 QPSolverRelaxedQPOPT::QPSolverRelaxedQPOPT(
156  value_type max_qp_iter_frac
157  )
158  :max_qp_iter_frac_(max_qp_iter_frac)
159  ,ITMAX_(100)
160  ,FEATOL_(1.0e-10)
161  ,LDH_(1)
162 {}
163 
164 QPSolverRelaxedQPOPT::~QPSolverRelaxedQPOPT()
165 {
166  release_memory();
167 }
168 
169 // Overridden from QPSolverRelaxed
170 
171 void QPSolverRelaxedQPOPT::release_memory()
172 {
173  // ToDo: Resize all arrays to zero!
175 }
176 
177 // Overridden protected members
178 
179 QPSolverRelaxedQPOPT::f_int QPSolverRelaxedQPOPT::liwork(f_int N, f_int NCLIN) const
180 { return 2* N + 3; }
181 
182 QPSolverRelaxedQPOPT::f_int QPSolverRelaxedQPOPT::lrwork(f_int N, f_int NCLIN) const
183 { return NCLIN > 0 ? 2*N*N + 8*N + 5*NCLIN : N*N + 8 *N; }
184 
185 QPSolverRelaxedQPOPTSOL::EInform QPSolverRelaxedQPOPT::call_qp_solver(bool warm_start)
186 {
187 
188  // Set the rest of the QPOPT inputs that could not be set in the constructor.
189 
190  BIGBND_ = this->infinite_bound();
191  ITMAX_ = max_qp_iter_frac() * N_;
192  LDA_ = ( A_.rows() > 0 ? A_.rows() : 1 );
193  H_ = reinterpret_cast<f_dbl_prec*>(this); // used to implement QPHESS
194  AX_.resize( NCLIN_ > 0 ? NCLIN_ : 1 );
195 
196  // Set option parameters
197  {
198  namespace ns = QPOPT_CppDecl;
199  namespace ft = FortranTypes;
200  ns::reset_defaults();
201  ns::set_logical_option( ns::WARM_START , warm_start ? ft::F_FALSE : ft::F_TRUE );
202  ns::set_real_option( ns::FEASIBILITY_TOLERANCE , FEATOL_ );
203  ns::set_real_option( ns::INFINITE_BOUND_SIZE , BIGBND_ );
204  ns::set_int_option( ns::ITERATION_LIMIT , ITMAX_ );
205  ns::set_int_option( ns::PRINT_FILE , 0 );
206  ns::set_int_option( ns::SUMMARY_FILE , 0 );
207  ns::set_int_option( ns::PRINT_LEVEL , 0 );
208  ns::set_int_option( ns::PROBLEM_TYPE , ns::QP2 );
209  }
210 
211  QPOPT_CppDecl::qpopt(
212  N_, NCLIN_, LDA_, LDH_, NCLIN_ ? &A_(1,1) : NULL, &BL_(1), &BU_(1)
213  , &CVEC_(1), H_
214  , FORTRAN_NAME_UL_(QPHESS_SERVER_RELAX2,qphess_server_relax2)
215  , &ISTATE_[0], &X_(1), INFORM_, ITER_, OBJ_, &AX_(1)
216  , &CLAMDA_(1), &IWORK_[0], LIWORK_, &WORK_[0], LWORK_ );
217 
218  EInform return_inform;
219  typedef QPSolverRelaxedQPOPTSOL bc;
220  switch(INFORM_) {
221  case STRONG_LOCAL_MIN:
222  return_inform = bc::STRONG_LOCAL_MIN;
223  break;
224  case WEAK_LOCAL_MIN:
225  return_inform = bc::WEAK_LOCAL_MIN;
226  break;
227  case UNBOUNDED:
229  true, Unbounded
230  ,"QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
231  " QPOPT returned that the QP is unbounded!" );
232  case INFEASIBLE:
234  true, Infeasible
235  ,"QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
236  " QPOPT returned that the QP is infeasible!" );
237  case ITMAX_EXCEEDED:
238  return_inform = bc::MAX_ITER_EXCEEDED;
239  break;
240  case MAX_DOF_TOO_SMALL:
242  true, std::runtime_error,
243  "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
244  " QPOPT says that the max dof is too small" );
245  case INVALID_INPUT:
247  true, InvalidInput,
248  "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
249  " QPOPT returned that the input is invalid" );
250  case PROB_TYPE_NOT_REGOG:
252  true, std::logic_error,
253  "QPSolverRelaxedQPOPT::call_qp_solver() : Error,"
254  " QPOPT says that the problem type is not recognized" );
255  break;
256  default:
257  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen
258  }
259 
260  return return_inform;
261 }
262 
263 } // end namespace ConstrainedOptPack
264 
265 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_QPOPT
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)