ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSolverRelaxedLOQO.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 // Here we map from the QPSolverRelaxed QP formulation to the LOQO QP formulation.
00043 //
00044 // QPSolverRelaxed QP formulation:
00045 // ------------------------------
00046 //
00047 // min     g'*d + 1/2 * d'*G*d + (eta + 1/2*eta^2) * M
00048 // s.t.    dL   <= d                   <= dU
00049 //         etaL <= eta
00050 //         eL   <= op(E)*d - b*eta     <= eU
00051 //                 op(F)*d + (1-eta)*f  = 0
00052 //
00053 // LOQO QP formulation:
00054 // -------------------
00055 //
00056 // min      c'*x + 1/2 * x'*Q*x
00057 // s.t.     b <= A*x <= b + r
00058 //          l <= x <= u
00059 //
00060 // Mapping =>
00061 //
00062 // LOQO   QPSolverRelaxed 
00063 // ----   ---------------
00064 // x      [ d; eta ]
00065 // c      [ g; M ]
00066 // Q      [ G, 0; 0, M ]
00067 // A      [ op(E), -b; op(F), -f ]
00068 // b      [ eL; -f ]
00069 // r      [ eU-eL; 0 ]
00070 // l      [ dL, etaL ]
00071 // u      [ dU, +inf ]
00072 //
00073 // Above, in the LOQO formulation all singly bounded inequalities
00074 // must be formulated as b(j) <= A(j,:)*x with r(j) = inf.  This
00075 // will require some fudging since eL(j) == -inf may be true in some
00076 // cases.  Here we will have to exchange eL(j) and eU(j) and use
00077 // A(j,:) = -op(E)(j,:).
00078 //
00079 
00080 #ifdef CONSTRAINED_OPTIMIZATION_PACK_USE_LOQO
00081 
00082 #include <assert.h>
00083 
00084 #include <vector>
00085 
00086 #include "ConstrainedOptPack_QPSolverRelaxedLOQO.hpp"
00087 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixExtractInvCholFactor.hpp"
00088 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00089 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOp.hpp"
00090 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00091 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00092 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_EtaVector.hpp"
00093 #include "AbstractLinAlgPack_sparse_bounds.hpp"
00094 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00095 #include "Midynamic_cast_verbose.h"
00096 #include "MiWorkspacePack.h"
00097 
00098 extern "C" {
00099 #include "loqo.h"     // -I$(LOQODIR)
00100 #include "myalloc.h"  // -I$(LOQODIR)
00101 } // end extern "C"
00102 
00103 namespace LinAlgOpPack {
00104   using AbstractLinAlgPack::Vp_StV;
00105   using AbstractLinAlgPack::Mp_StM;
00106   using AbstractLinAlgPack::Vp_StMtV;
00107 }
00108 
00109 namespace ConstrainedOptPack {
00110 
00111 // ///////////////////////////////////////
00112 // Members for QPSolverRelaxedLOQO::InitLOQOHessianJacobian
00113 
00114 void QPSolverRelaxedLOQO::InitLOQOHessianJacobian::init_hess_jacob(
00115   const MatrixOp& G, const value_type bigM
00116   , const MatrixOp* E, BLAS_Cpp::Transp trans_E, const DVectorSlice* b
00117   , const int loqo_b_stat[], const size_type num_inequal
00118   , const MatrixOp* F, BLAS_Cpp::Transp trans_F, const DVectorSlice* f
00119   , void* _loqo_lp
00120   ) const
00121 {
00122 
00123   LOQO* loqo_lp = (LOQO*)_loqo_lp;
00124 
00125   const size_type
00126     nd    = G.rows(),
00127     m_in  = E ? b->size() : 0,
00128     m_eq  = F ? f->size() : 0;
00129 
00130   TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->n == nd + 1  ) );
00131   TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->m == num_inequal + m_eq  ) );
00132 
00133   // This default implementation assumes G, E and F are completely dense!
00134 
00135   //
00136   // Setup Q
00137   //
00138 
00139   loqo_lp->qnz = nd*nd + 1;
00140   MALLOC( loqo_lp->Q, loqo_lp->qnz, double );
00141   MALLOC( loqo_lp->iQ, loqo_lp->qnz, int );
00142   MALLOC( loqo_lp->kQ, nd+2, int );
00143   // Setup kQ[] and iQ[]
00144   {for( size_type j = 1; j <= nd; ++j ) {
00145     loqo_lp->kQ[j-1] = nd*(j-1);
00146     for( size_type i = 1; i <= nd; ++i )
00147       loqo_lp->iQ[ loqo_lp->kQ[j-1] + (i-1) ] = i-1; // zero based in LOQO
00148   }}
00149   loqo_lp->kQ[nd]              = nd*nd;
00150   loqo_lp->iQ[loqo_lp->kQ[nd]] = nd; // zero based in LOQO
00151   loqo_lp->kQ[nd+1]            = nd*nd + 1;
00152   // Setup Q[]
00153   {
00154     DMatrixSlice Q( loqo_lp->Q, nd*nd, nd, nd, nd );
00155     LinAlgOpPack::assign( &Q, G, BLAS_Cpp::no_trans );
00156     loqo_lp->Q[nd*nd] = bigM;
00157   }
00158 
00159   //
00160   // Setup A
00161   //
00162 
00163   loqo_lp->nz = (num_inequal+m_eq) * (nd+1);
00164   MALLOC( loqo_lp->A,  loqo_lp->nz, double );
00165   MALLOC( loqo_lp->iA, loqo_lp->nz, int );
00166   MALLOC( loqo_lp->kA, nd+2, int );
00167 
00168   if( num_inequal == m_in ) {
00169     // All the inequalities have finite bounds
00170     // Setup kA[] and iA[]
00171     {for( size_type j = 1; j <= nd+1; ++j ) {
00172       loqo_lp->kA[j-1] = (m_in+m_eq)*(j-1);
00173       for( size_type i = 1; i <= m_in+m_eq; ++i )
00174         loqo_lp->iA[ loqo_lp->kA[j-1] + (i-1) ] = i-1; // zero based in LOQO
00175     }}
00176     loqo_lp->kA[nd+1] = (m_in+m_eq)*(nd+1);
00177     // Setup A[]
00178     DMatrixSlice A( loqo_lp->A, loqo_lp->nz, loqo_lp->m, loqo_lp->m, nd+1 );
00179     if(E) {
00180       LinAlgOpPack::assign( &A(1,m_in,1,nd), *E, trans_E );  // A(1:m_in,1:nd) = op(E)
00181       LinAlgOpPack::V_StV( &A.col(nd+1)(1,m_in), -1.0, *b ); // A(1:m_in,nd+1) = -b
00182     }
00183     if(F) {
00184       LinAlgOpPack::assign( &A(m_in+1,m_in+m_eq,1,nd), *F, trans_F );  // A(m_in+1:m_in+m_eq,1:nd) = op(F)
00185       LinAlgOpPack::V_StV( &A.col(nd+1)(m_in+1,m_in+m_eq), -1.0, *f ); // A(m_in+1:m_in+m_eq,nd+1) = -f
00186     }
00187   }
00188   else {
00189     // At least one of the inequality constriants has
00190     // both infinite upper and lower bounds.
00191     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Finish this!
00192   }
00193   
00194   // Loop through and adjust A for absent lower bound and using upper bound
00195   if( num_inequal ) {
00196     DMatrixSlice A( loqo_lp->A, loqo_lp->nz, loqo_lp->m, loqo_lp->m, nd+1 );
00197     for(size_type k = 1; k <= num_inequal; ++k ) {
00198       const int j = loqo_b_stat[k-1];
00199       if( j < 0 )
00200         DenseLinAlgPack::Vt_S( &A.row(j), -1.0 );
00201     }
00202   }
00203 
00204 }
00205 
00206 // ///////////////////////////////////////
00207 // Members for QPSolverRelaxedLOQO
00208 
00209 QPSolverRelaxedLOQO::QPSolverRelaxedLOQO(
00210   const init_hess_jacob_ptr_t  init_hess_jacob
00211   ,value_type                  bigM
00212   ,value_type                  nonbinding_lag_mult
00213   )
00214   :init_hess_jacob_(init_hess_jacob)
00215   ,bigM_(bigM)
00216   ,nonbinding_lag_mult_(nonbinding_lag_mult)
00217 {
00218 //  bigM_ = 1.0; // Just test this!
00219   nonbinding_lag_mult_ = 1e-6;
00220 }
00221 
00222 QPSolverRelaxedLOQO::~QPSolverRelaxedLOQO()
00223 {
00224   this->release_memory();
00225 }
00226 
00227 // Overridden from QPSolverRelaxed
00228 
00229 QPSolverStats
00230 QPSolverRelaxedLOQO::get_qp_stats() const
00231 {
00232   return qp_stats_;
00233 }
00234 
00235 void QPSolverRelaxedLOQO::release_memory()
00236 {
00237   // Todo: resize to zero all the workspace!
00238 }
00239 
00240 QPSolverStats::ESolutionType
00241 QPSolverRelaxedLOQO::imp_solve_qp(
00242       std::ostream* out, EOutputLevel olevel, ERunTests test_what
00243     , const DVectorSlice& g, const MatrixOp& G
00244     , value_type etaL
00245     , const SpVectorSlice& dL, const SpVectorSlice& dU
00246     , const MatrixOp* E, BLAS_Cpp::Transp trans_E, const DVectorSlice* b
00247       , const SpVectorSlice* eL, const SpVectorSlice* eU
00248     , const MatrixOp* F, BLAS_Cpp::Transp trans_F, const DVectorSlice* f
00249     , value_type* obj_d
00250     , value_type* eta, DVectorSlice* d
00251     , SpVector* nu
00252     , SpVector* mu, DVectorSlice* Ed
00253     , DVectorSlice* lambda, DVectorSlice* Fd
00254   )
00255 {
00256   using Teuchos::Workspace;
00257   Teuchos::WorkspaceStore* wss = wsp::default_workspace_store.get();
00258 
00259   const value_type inf_bnd  = std::numeric_limits<value_type>::max();
00260 //  const value_type real_big = 1e+20;
00261   const value_type real_big = HUGE_VAL;
00262 
00263   const size_type
00264     nd   = g.size(),
00265     m_in = E ? b->size() : 0,
00266     m_eq = F ? f->size() : 0;
00267 
00268   //
00269   // Create a LOQO QP definition struct
00270   //
00271 
00272   LOQO *loqo_lp = openlp();
00273   TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp  ) );
00274 
00275   //
00276   // Setup loqo_r and loqo_b and count the number of actual
00277   // constraints.
00278   //
00279 
00280   // LOQO's b vector storage
00281   MALLOC( loqo_lp->b, m_in+m_eq, double ); // May not use all of this storage
00282   DVectorSlice loqo_b( loqo_lp->b, m_in+m_eq );
00283   // LOQO's r vector storage
00284   MALLOC( loqo_lp->r, m_in+m_eq, double ); // May not use all of this storage
00285   DVectorSlice loqo_r( loqo_lp->r, m_in+m_eq );
00286   // Gives status of b.
00287   //                  /  j : if eL(j) > -inf_bnd
00288   // loqo_b_stat(k) = |
00289   //                  \ -j : if eL(j) <= -inf_bnd && eU(j) < +inf_bnd
00290   //
00291   // , for k = 1...num_inequal
00292   //
00293   Workspace<int>               loqo_b_stat_ws(wss,m_in); // May not use all of this
00294   DenseLinAlgPack::VectorSliceTmpl<int>  loqo_b_stat(&loqo_b_stat_ws[0],loqo_b_stat_ws.size());
00295   std::fill( loqo_b_stat.begin(), loqo_b_stat.end(), 0 ); // Initialize to zero
00296 
00297   // Fill up loqo_b, loqo_r and loqo_b_stat
00298   size_type num_inequal = 0; // The actual number of bouned general inequalities
00299   if(E) {
00300     // Read iterators
00301     AbstractLinAlgPack::sparse_bounds_itr
00302       eLU_itr( eL->begin(), eL->end(), eL->offset()
00303            , eU->begin(), eU->end(), eU->offset(), inf_bnd );
00304     // written iterators
00305     DVectorSlice::iterator
00306       b_itr   = loqo_b.begin(),
00307       r_itr   = loqo_r.begin();
00308     DenseLinAlgPack::VectorSliceTmpl<int>::iterator
00309       b_stat_itr  = loqo_b_stat.begin();
00310     // loop
00311     for( int k = 1; !eLU_itr.at_end(); ++k, ++eLU_itr, ++b_itr, ++r_itr, ++b_stat_itr, ++num_inequal )
00312     {
00313       const size_type j = eLU_itr.indice();
00314       if(eLU_itr.lbound() > -inf_bnd) {
00315         *b_itr = eLU_itr.lbound();
00316         *r_itr = eLU_itr.ubound() >= inf_bnd ? real_big : eLU_itr.ubound() - eLU_itr.lbound();
00317         *b_stat_itr = j; // We need to make A(k,:) = [ +op(E)(j,:), -b(j) ]
00318       }
00319       else {
00320         TEUCHOS_TEST_FOR_EXCEPT( !( eLU_itr.ubound() < +inf_bnd ) );
00321         *b_itr = -eLU_itr.ubound();
00322         *r_itr = eLU_itr.lbound() <= -inf_bnd ? real_big : - eLU_itr.lbound() + eLU_itr.ubound();
00323         *b_stat_itr = -j; // We need to make A(k,:) = [ -op(E)(j,:), +b(j) ]
00324       }
00325     }
00326   }
00327   if(F) {
00328     LinAlgOpPack::V_StV( &loqo_b(num_inequal+1,num_inequal+m_eq), -1.0, *f );
00329     loqo_r(num_inequal+1,num_inequal+m_eq) = 0.0;
00330   }
00331 
00332   //
00333   // Setup the QP dimensions
00334   //
00335 
00336   loqo_lp->n = nd+1;
00337   loqo_lp->m = num_inequal + m_eq;
00338 
00339   //
00340   // Setup loqo_c, loqo_l and loqo_u
00341   //
00342 
00343   // LOQO's c vector storage
00344   MALLOC( loqo_lp->c, nd+1, double );
00345   DVectorSlice loqo_c( loqo_lp->c, nd+1 );
00346   loqo_c(1,nd) = g;
00347   loqo_c(nd+1) = bigM();
00348 
00349   // LOQO's l vector storage
00350   MALLOC( loqo_lp->l, nd+1, double );
00351   DVectorSlice loqo_l( loqo_lp->l, nd+1 );
00352   std::fill( loqo_l.begin(), loqo_l.end(), -real_big );
00353   {
00354     SpVectorSlice::const_iterator
00355       dL_itr = dL.begin(),
00356       dL_end = dL.end();
00357     for( ; dL_itr != dL_end; ++dL_itr )
00358       loqo_l( dL_itr->indice() + dL.offset() ) = dL_itr->value();
00359   }
00360   loqo_l(nd+1) = etaL;
00361 
00362   // LOQO's u vector storage
00363   MALLOC( loqo_lp->u, nd+1, double );
00364   DVectorSlice loqo_u( loqo_lp->u, nd+1 );
00365   std::fill( loqo_u.begin(), loqo_u.end(), +real_big );
00366   {
00367     SpVectorSlice::const_iterator
00368       dU_itr = dU.begin(),
00369       dU_end = dU.end();
00370     for( ; dU_itr != dU_end; ++dU_itr )
00371       loqo_u( dU_itr->indice() + dU.offset() ) = dU_itr->value();
00372   }
00373   loqo_u(nd+1) = +real_big;
00374   
00375   //
00376   // Setup the objective and constraint matrices (using strategy interface).
00377   //
00378 
00379   init_hess_jacob().init_hess_jacob(
00380     G,bigM(),E,trans_E,b,&loqo_b_stat[0],num_inequal,F,trans_F,f
00381     ,loqo_lp);
00382 
00383   //
00384   // Setup the starting point
00385   //
00386 
00387   MALLOC( loqo_lp->x, nd+1, double );
00388   DVectorSlice loqo_x( loqo_lp->x, nd+1 );
00389   loqo_x(1,nd) = *d;
00390   loqo_x(nd+1) = *eta;
00391 
00392   //
00393   // Set some control parameters
00394   //
00395   
00396 //  strcpy( loqo_lp->name, "loqo_qp" );
00397   loqo_lp->quadratic = 1;
00398   loqo_lp->convex    = 1;
00399   switch( olevel ) {
00400     case PRINT_NONE:
00401       loqo_lp->verbose = 0;
00402       break;
00403     case PRINT_BASIC_INFO:
00404       loqo_lp->verbose = 1;
00405       break;
00406     case PRINT_ITER_SUMMARY:
00407       loqo_lp->verbose = 2;
00408       break;
00409     case PRINT_ITER_STEPS:
00410       loqo_lp->verbose = 3;
00411       break;
00412     case PRINT_ITER_ACT_SET:
00413       loqo_lp->verbose = 4;
00414       break;
00415     case PRINT_ITER_VECTORS:
00416       loqo_lp->verbose = 5;
00417       break;
00418     case PRINT_EVERY_THING:
00419       loqo_lp->verbose = 6;
00420       break;
00421     default:
00422       TEUCHOS_TEST_FOR_EXCEPT(true);
00423   }
00424 
00425   //
00426   // Solve the QP
00427   //
00428 
00429   if( out && olevel >= PRINT_BASIC_INFO ) {
00430     *out << "\nSolving QP using LOQO ...\n";
00431     out->flush();
00432   }
00433   
00434   const int loqo_status = solvelp(loqo_lp);
00435 
00436   if( out && olevel >= PRINT_BASIC_INFO ) {
00437     *out << "\nLOQO returned status = " << loqo_status << "\n";
00438   }
00439 
00440   //
00441   // Map the solution to the output arguments
00442   //
00443 
00444   TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->x  ) );
00445   DVectorSlice loqo_x_sol( loqo_lp->x, nd+1 );
00446 
00447   // d
00448   *d    = loqo_x_sol(1,nd);
00449 
00450   // eta
00451   *eta  = loqo_x_sol(nd+1);
00452 
00453   // obj_d
00454   if(obj_d)
00455     *obj_d = loqo_lp->primal_obj - (*eta + 0.5 * (*eta)*(*eta)) * bigM();
00456 
00457   // nu
00458   if(nu) {
00459     nu->resize(nd,nd);
00460     TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->z  ) );
00461     TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->s  ) );
00462     const DVectorSlice
00463       loqo_z(loqo_lp->z,loqo_lp->n),   // Multipliers for l - x <= 0
00464       loqo_s(loqo_lp->s,loqo_lp->n);   // Multipliers for x - u <= 0
00465     DVectorSlice::const_iterator
00466       z_itr = loqo_z.begin(),
00467       s_itr = loqo_s.begin();
00468     typedef SpVector::element_type ele_t;
00469     for( size_type i = 1; i <= nd; ++i, ++z_itr, ++s_itr ) {
00470       if( *z_itr > *s_itr && *z_itr >= nonbinding_lag_mult() ) {
00471         // Lower bound is active
00472         nu->add_element(ele_t(i,-(*z_itr)));
00473       }
00474       else if( *s_itr > *z_itr && *s_itr >= nonbinding_lag_mult() ) {
00475         // Upper bound is active
00476         nu->add_element(ele_t(i,+(*s_itr)));
00477       }
00478     }
00479     // We could look at z(nd+1) and s(nd+1) for the value of kappa?
00480     nu->assume_sorted(true);
00481   }
00482 
00483   // mu
00484   if(mu) {
00485     mu->resize(m_in,num_inequal);
00486     DenseLinAlgPack::VectorSliceTmpl<int>::iterator
00487       b_stat_itr  = loqo_b_stat.begin();
00488     TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->v  ) );
00489     TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->q  ) );
00490     const DVectorSlice
00491       loqo_v(loqo_lp->v,loqo_lp->m),   // Multipliers for b <= A*x
00492       loqo_q(loqo_lp->q,loqo_lp->m);   // Multipliers for A*x <= b + r
00493     DVectorSlice::const_iterator
00494       v_itr = loqo_v.begin(),
00495       q_itr = loqo_q.begin();
00496     // loop
00497     typedef SpVector::element_type ele_t;
00498     for( size_type k = 1; k <= num_inequal; ++k, ++b_stat_itr, ++v_itr, ++q_itr ) {
00499       const int j = *b_stat_itr;
00500       if( *v_itr > *q_itr && *v_itr >= nonbinding_lag_mult() ) {
00501         // Lower bound is active
00502         if( j < 0 ) // We had to flip this since it was really and upper bound
00503           mu->add_element(ele_t(-j,+(*v_itr)));
00504         else // This really was a lower bound
00505           mu->add_element(ele_t(+j,-(*v_itr)));
00506       }
00507       else if( *q_itr > *v_itr && *q_itr >= nonbinding_lag_mult() ) {
00508         // Upper bound is active
00509         mu->add_element(ele_t(+j,+(*q_itr)));
00510       }
00511     }
00512   }
00513 
00514   // Ed
00515   if(Ed) {
00516     LinAlgOpPack::V_MtV( Ed, *E, trans_E, *d );
00517   }
00518 
00519   // lambda
00520   if(lambda) {
00521     TEUCHOS_TEST_FOR_EXCEPT( !(  loqo_lp->y  ) );
00522     const DVectorSlice
00523       loqo_y(loqo_lp->y,loqo_lp->m);         // Multipliers for equalities
00524     DVectorSlice::const_iterator
00525       y_itr = loqo_y.begin() + num_inequal;  // Get iterators to equalities
00526     DVectorSlice::iterator
00527       lambda_itr = lambda->begin();
00528     // loop
00529     for( size_type k = 1; k <= m_eq; ++k, ++y_itr, ++lambda_itr ) {
00530       *lambda_itr = -(*y_itr);
00531     }
00532   }
00533 
00534   // Fd
00535   if(Fd) {
00536     LinAlgOpPack::V_MtV( Fd, *F, trans_F, *d );
00537   }
00538 
00539   //
00540   // Setup the QP statistics
00541   //
00542 
00543   QPSolverStats::ESolutionType solution_type = QPSolverStats::OPTIMAL_SOLUTION; // Assume this?
00544   switch( loqo_status ) { // I had to find this out by trial and error!
00545       case 0:
00546       solution_type = QPSolverStats::OPTIMAL_SOLUTION;
00547       break;
00548     case 2:
00549       solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
00550       break;
00551     default:
00552       TEUCHOS_TEST_FOR_EXCEPT(true);
00553   }
00554 
00555   qp_stats_.set_stats(
00556     solution_type, QPSolverStats::CONVEX
00557     ,loqo_lp->iter, QPSolverStats::NOT_KNOWN, QPSolverStats::NOT_KNOWN
00558     ,false, *eta > 0.0 );
00559 
00560   //
00561   // Clean up dynamically allocated memory for LOQO
00562   //
00563 
00564   inv_clo();          // frees memory associated with matrix factorization
00565   closelp(loqo_lp);   // frees all allocated arrays with free(...).
00566 
00567   return qp_stats_.solution_type();
00568 
00569 }
00570 
00571 } // end namespace ConstrainedOptPack
00572 
00573 #endif // CONSTRAINED_OPTIMIZATION_PACK_USE_LOQO
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends