ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_QPSolverRelaxed.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 #include <assert.h>
00043 
00044 #include <limits>
00045 
00046 #include "ConstrainedOptPack_QPSolverRelaxed.hpp"
00047 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00048 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00049 #include "AbstractLinAlgPack_VectorMutable.hpp"
00050 #include "AbstractLinAlgPack_VectorOut.hpp"
00051 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00052 #include "ProfileHackPack_profile_hack.hpp"
00053 #include "Teuchos_Assert.hpp"
00054 
00055 namespace ConstrainedOptPack {
00056 
00057 // public
00058 
00059 QPSolverRelaxed::QPSolverRelaxed()
00060   :infinite_bound_(std::numeric_limits<value_type>::max())
00061 {}
00062 
00063 QPSolverStats::ESolutionType
00064 QPSolverRelaxed::solve_qp(
00065   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00066   ,const Vector& g, const MatrixSymOp& G
00067   ,value_type etaL
00068   ,const Vector& dL, const Vector& dU
00069   ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
00070   ,const Vector& eL, const Vector& eU
00071   ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
00072   ,value_type* obj_d
00073   ,value_type* eta, VectorMutable* d
00074   ,VectorMutable* nu
00075   ,VectorMutable* mu, VectorMutable* Ed
00076   ,VectorMutable* lambda, VectorMutable* Fd
00077   )
00078 {
00079   return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
00080     ,&E,trans_E,&b,&eL,&eU,&F,trans_F,&f
00081     ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00082 }
00083 
00084 QPSolverStats::ESolutionType
00085 QPSolverRelaxed::solve_qp(
00086   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00087   ,const Vector& g, const MatrixSymOp& G
00088   ,value_type etaL
00089   ,const Vector& dL, const Vector& dU
00090   ,const MatrixOp& E, BLAS_Cpp::Transp trans_E, const Vector& b
00091   ,const Vector& eL, const Vector& eU
00092   ,value_type* obj_d
00093   ,value_type* eta, VectorMutable* d
00094   ,VectorMutable* nu
00095   ,VectorMutable* mu, VectorMutable* Ed
00096   )
00097 {
00098   return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
00099     ,&E,trans_E,&b,&eL,&eU,NULL,BLAS_Cpp::no_trans,NULL
00100     ,obj_d,eta,d,nu,mu,Ed,NULL,NULL);
00101 }
00102 
00103 QPSolverStats::ESolutionType
00104 QPSolverRelaxed::solve_qp(
00105   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00106   ,const Vector& g, const MatrixSymOp& G
00107   ,value_type etaL
00108   ,const Vector& dL, const Vector& dU
00109   ,const MatrixOp& F, BLAS_Cpp::Transp trans_F, const Vector& f
00110   ,value_type* obj_d
00111   ,value_type* eta, VectorMutable* d
00112   ,VectorMutable* nu
00113   ,VectorMutable* lambda, VectorMutable* Fd
00114   )
00115 {
00116   return solve_qp(out,olevel,test_what,g,G,etaL,&dL,&dU
00117     ,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL,&F,trans_F,&f
00118     ,obj_d,eta,d,nu,NULL,NULL,lambda,Fd);
00119 }
00120 
00121 
00122 QPSolverStats::ESolutionType
00123 QPSolverRelaxed::solve_qp(
00124   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00125   ,const Vector& g, const MatrixSymOp& G
00126   ,const Vector& dL, const Vector& dU
00127   ,value_type* obj_d
00128   ,VectorMutable* d
00129   ,VectorMutable* nu
00130   )
00131 {
00132   value_type eta = 0, etaL = 0;
00133   return solve_qp(
00134     out,olevel,test_what,g,G,etaL,&dL,&dU
00135     ,NULL,BLAS_Cpp::no_trans,NULL,NULL,NULL
00136     ,NULL,BLAS_Cpp::no_trans,NULL
00137     ,obj_d,&eta,d,nu,NULL,NULL,NULL,NULL);
00138 }
00139 
00140 QPSolverStats::ESolutionType
00141 QPSolverRelaxed::solve_qp(
00142   std::ostream* out, EOutputLevel olevel, ERunTests test_what
00143   ,const Vector& g, const MatrixSymOp& G
00144   ,value_type etaL
00145   ,const Vector* dL, const Vector* dU
00146   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00147   ,const Vector* eL, const Vector* eU
00148   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00149   ,value_type* obj_d
00150   ,value_type* eta, VectorMutable* d
00151   ,VectorMutable* nu
00152   ,VectorMutable* mu, VectorMutable* Ed
00153   ,VectorMutable* lambda, VectorMutable* Fd
00154   )
00155 {
00156 #ifdef PROFILE_HACK_ENABLED
00157   ProfileHackPack::ProfileTiming profile_timing( "QPSolverRelaxed::solve_qp(...)" );
00158 #endif
00159   validate_input(
00160     infinite_bound(),g,G,etaL,dL,dU
00161     ,E,trans_E,b,eL,eU,F,trans_F,f
00162     ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00163   print_qp_input(
00164     infinite_bound(),out,olevel,g,G,etaL,dL,dU,E,trans_E,b,eL,eU
00165     ,F,trans_F,f,eta,d,nu,mu,lambda );
00166   QPSolverStats::ESolutionType
00167     solve_return = imp_solve_qp(
00168       out,olevel,test_what,g,G,etaL,dL,dU
00169       ,E,trans_E,b,eL,eU,F,trans_F,f
00170       ,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00171   print_qp_output(
00172     infinite_bound(),out,olevel,obj_d,eta,d,nu,mu,Ed,lambda,Fd);
00173   return solve_return;
00174 }
00175 
00176 void QPSolverRelaxed::validate_input(
00177   const value_type infinite_bound
00178   ,const Vector& g, const MatrixSymOp& G
00179   ,value_type etaL
00180   ,const Vector* dL, const Vector* dU
00181   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00182   ,const Vector* eL, const Vector* eU
00183   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00184   ,const value_type* obj_d
00185   ,const value_type* eta, const Vector* d
00186   ,const Vector* nu
00187   ,const Vector* mu, const Vector* Ed
00188   ,const Vector* lambda, const Vector* Fd
00189   )
00190 {
00191   // Validate output arguments
00192   TEUCHOS_TEST_FOR_EXCEPTION(
00193     !d, std::invalid_argument
00194     ,"QPSolverRelaxed::validate_input(...) : Error, "
00195     "If d!=NULL is not allowed." );
00196   TEUCHOS_TEST_FOR_EXCEPTION(
00197     ( E || F ) && !eta, std::invalid_argument
00198     ,"QPSolverRelaxed::validate_input(...) : Error, "
00199     "If eta!=NULL is not allowed if E!=NULL or F!=NULL." );
00200 
00201   // Validate the sets of constraints arguments
00202   TEUCHOS_TEST_FOR_EXCEPTION(
00203     dL && ( !dU || !nu ), std::invalid_argument
00204     ,"QPSolverRelaxed::validate_input(...) : Error, "
00205     "If dL!=NULL then dU!=NULL and nu!=NULL must also be true." );
00206   TEUCHOS_TEST_FOR_EXCEPTION(
00207     E && ( !b || !eL || !eU || !mu ), std::invalid_argument
00208     ,"QPSolverRelaxed::validate_input(...) : Error, "
00209     "If E!=NULL then b!=NULL, eL!=NULL, eU!=NULL and mu!=NULL must also "
00210     "be true." );
00211   TEUCHOS_TEST_FOR_EXCEPTION(
00212     F && ( !f || !lambda ), std::invalid_argument
00213     ,"QPSolverRelaxed::validate_input(...) : Error, "
00214     "If F!=NULL then f!=NULL and lambda!=NULL must also "
00215     "be true." );
00216 
00217   // ToDo: Replace the below code with checks of compatibility of the vector spaces!
00218 
00219   // Validate the sizes of the arguments
00220   const size_type
00221     nd = d->dim();
00222   TEUCHOS_TEST_FOR_EXCEPTION(
00223     g.dim() != nd, std::invalid_argument
00224     ,"QPSolverRelaxed::validate_input(...) : Error, "
00225     "g.dim() != d->dim()." );
00226   TEUCHOS_TEST_FOR_EXCEPTION(
00227     G.rows() != nd || G.cols() != nd, std::invalid_argument
00228     ,"QPSolverRelaxed::validate_input(...) : Error, "
00229     "G.rows() != d->dim() or G.cols() != d->dim()." );
00230   TEUCHOS_TEST_FOR_EXCEPTION(
00231     dL && dL->dim() != nd, std::invalid_argument
00232     ,"QPSolverRelaxed::validate_input(...) : Error, "
00233     "dL->dim() = " << dL->dim() << " != d->dim() = " << nd << "." );
00234   TEUCHOS_TEST_FOR_EXCEPTION(
00235     dU && dU->dim() != nd, std::invalid_argument
00236     ,"QPSolverRelaxed::validate_input(...) : Error, "
00237     "dU->dim() = " << dU->dim() << " != d->dim() = " << nd << "." );
00238   if( E ) {
00239     const size_type
00240       m_in = BLAS_Cpp::rows( E->rows(), E->cols(), trans_E );
00241     TEUCHOS_TEST_FOR_EXCEPTION(
00242       BLAS_Cpp::cols( E->rows(), E->cols(), trans_E ) != nd, std::invalid_argument
00243       ,"QPSolverRelaxed::validate_input(...) : Error, op(E).cols() != d->dim()." );
00244     TEUCHOS_TEST_FOR_EXCEPTION(
00245       b->dim() != m_in, std::invalid_argument
00246       ,"QPSolverRelaxed::validate_input(...) : Error, b->dim() != op(E).rows()." );
00247     TEUCHOS_TEST_FOR_EXCEPTION(
00248       eL->dim() != m_in, std::invalid_argument
00249       ,"QPSolverRelaxed::validate_input(...) : Error, eL->dim() != op(E).rows()." );
00250     TEUCHOS_TEST_FOR_EXCEPTION(
00251       eU->dim() != m_in, std::invalid_argument
00252       ,"QPSolverRelaxed::validate_input(...) : Error, eU->dim() != op(E).rows()." );
00253     TEUCHOS_TEST_FOR_EXCEPTION(
00254       Ed && Ed->dim() != m_in, std::invalid_argument
00255       ,"QPSolverRelaxed::validate_input(...) : Error, Ed->dim() != op(E).rows()." );
00256   }
00257   if( F ) {
00258     const size_type
00259       m_eq = BLAS_Cpp::rows( F->rows(), F->cols(), trans_F );
00260     TEUCHOS_TEST_FOR_EXCEPTION(
00261       BLAS_Cpp::cols( F->rows(), F->cols(), trans_F ) != nd, std::invalid_argument
00262       ,"QPSolverRelaxed::validate_input(...) : Error, op(F).cols() != d->dim()." );
00263     TEUCHOS_TEST_FOR_EXCEPTION(
00264       f->dim() != m_eq, std::invalid_argument
00265       ,"QPSolverRelaxed::validate_input(...) : Error, f->dim() != op(F).rows()." );
00266     TEUCHOS_TEST_FOR_EXCEPTION(
00267       lambda->dim() != m_eq, std::invalid_argument
00268       ,"QPSolverRelaxed::validate_input(...) : Error, lambda->dim() != op(F).rows()." );
00269     TEUCHOS_TEST_FOR_EXCEPTION(
00270       Fd && Fd->dim() != m_eq, std::invalid_argument
00271       ,"QPSolverRelaxed::validate_input(...) : Error, Fd->dim() != op(F).rows()." );
00272   }
00273 }
00274 
00275 void QPSolverRelaxed::print_qp_input( 
00276   const value_type infinite_bound
00277   ,std::ostream* out, EOutputLevel olevel
00278   ,const Vector& g, const MatrixSymOp& G
00279   ,value_type etaL
00280   ,const Vector* dL, const Vector* dU
00281   ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
00282   ,const Vector* eL, const Vector* eU
00283   ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
00284   ,value_type* eta, VectorMutable* d
00285   ,VectorMutable* nu
00286   ,VectorMutable* mu
00287   ,VectorMutable* lambda
00288   )
00289 {
00290   using AbstractLinAlgPack::num_bounded;
00291   if( out && (int)olevel >= (int)PRINT_ITER_STEPS ) {
00292     *out<< "\n*** Printing input to QPSolverRelaxed::solve_qp(...) ...\n";
00293     // g
00294     *out << "\n||g||inf = " << g.norm_inf() << std::endl;
00295     if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00296       *out<< "g =\n" << g;
00297     // G
00298     if( (int)olevel >= (int)PRINT_EVERY_THING )
00299       *out<< "\nG =\n" << G;
00300     // etaL
00301     *out << "\netaL = " << etaL << std::endl;
00302     // eta
00303     *out << "\neta  = " << *eta << std::endl;
00304     if(dL) {
00305       // dL, dU
00306       *out << "\ndL.dim()   = " << dL->dim();
00307       *out << "\ndU.dim()   = " << dU->dim();
00308       *out << "\nnum_bounded(dL,dU) = " << num_bounded(*dL,*dU,infinite_bound) << std::endl;
00309       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00310         *out << "dL =\n" << *dL;
00311       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00312         *out << "dU =\n" << *dU;
00313     }
00314     else {
00315       *out << "\ndL = -inf";
00316       *out << "\ndU = +inf";
00317     }
00318     // d
00319     *out << "\n||d||inf = " << d->norm_inf() << std::endl;
00320     if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00321       *out<< "d =\n" << *d;
00322     // nu
00323     if(nu) {
00324       *out<< "\nnu->nz() = " << nu->nz() << std::endl
00325         << "||nu||inf  = " << nu->norm_inf() << std::endl;
00326       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00327         *out<< "nu =\n" << *nu;
00328     }
00329     if(E) {
00330       // op(E)
00331       if( (int)olevel >= (int)PRINT_EVERY_THING )
00332         *out<< "\nE" << std::endl << *E
00333           << "trans_E = " << BLAS_Cpp::trans_to_string(trans_E) << std::endl;
00334       // b
00335       *out << "\n||b||inf = " << b->norm_inf() << std::endl;
00336       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00337       *out<< "b =\n" << *b;
00338       // eL, eU
00339       *out<< "\neL.dim()   = " << eL->dim();
00340       *out<< "\neU.dim()   = " << eU->dim();
00341       *out << "\nnum_bounded(eL,eU) = " << num_bounded(*eL,*eU,infinite_bound) << std::endl;
00342       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00343         *out<< "eL =\n" << *eL;
00344       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00345         *out<< "eU =\n" << *eU;
00346       // mu
00347       *out<< "\nmu.nz()   = " << mu->nz() << std::endl
00348         << "||mu||inf = " << mu->norm_inf() << std::endl;
00349       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00350         *out<< "mu =\n" << *mu;
00351     }
00352     if(F) {
00353       // op(F)
00354       if( (int)olevel >= (int)PRINT_EVERY_THING )
00355         *out<< "\nF" << std::endl << *F
00356           << "trans_F = " << BLAS_Cpp::trans_to_string(trans_F) << std::endl;
00357       // f
00358       *out<< "\n||f||inf = " << f->norm_inf() << std::endl;
00359       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00360         *out<< "f =\n" << *f;
00361       // lambda
00362       *out<< "\n||lambda||inf = " << lambda->norm_inf() << std::endl;
00363       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00364         *out<< "lambda =\n" << *lambda;
00365     }
00366     *out<< "\nEnd input to QPSolverRelaxed::solve_qp(...)\n";
00367   }
00368 }
00369 
00370 void QPSolverRelaxed::print_qp_output(
00371   const value_type infinite_bound
00372   ,std::ostream* out, EOutputLevel olevel
00373   ,const value_type* obj_d
00374   ,const value_type* eta, const Vector* d
00375   ,const Vector* nu
00376   ,const Vector* mu, const Vector* Ed
00377   ,const Vector* lambda, const Vector* Fd
00378   )
00379 {
00380   if( out && (int)olevel > (int)PRINT_ITER_STEPS ) {
00381     *out<< "\n*** Printing output from QPSolverRelaxed::solve_qp(...) ...\n";
00382     // obj_d
00383     if(obj_d)
00384       *out << "\nobj_d = " << *obj_d << std::endl;
00385     // eta
00386     *out << "\neta = " << *eta << std::endl;
00387     // d
00388     *out << "\n||d||inf = " << d->norm_inf() << std::endl;
00389     if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00390       *out<< "d =\n" << *d;
00391     // nu
00392     if(nu) {
00393       *out<< "\nnu.nz()   = " << nu->nz() << std::endl
00394         << "||nu||inf = " << nu->norm_inf() << std::endl;
00395       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00396         *out<< "nu =\n" << *nu;
00397     }
00398     // Ed
00399     if(Ed) {
00400       *out << "\n||Ed||inf = " << Ed->norm_inf() << std::endl;
00401       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00402         *out<< "Ed =\n" << *Ed;
00403     }
00404     // mu
00405     if(mu) {
00406       *out<< "\nmu.nz()   = " << mu->nz() << std::endl
00407         << "||mu||inf = " << mu->norm_inf() << std::endl;
00408       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00409         *out<< "mu =\n" << *mu;
00410     }
00411     // lambda
00412     if(lambda) {
00413       *out<< "\n||lambda||inf = " << lambda->norm_inf() << std::endl;
00414       if( (int)olevel >= (int)PRINT_ITER_ACT_SET )
00415         *out<< "lambda =\n" << *lambda;
00416     }
00417     // Fd
00418     if(Fd) {
00419       *out << "\n||Fd||inf = " << Fd->norm_inf() << std::endl;
00420       if( (int)olevel >= (int)PRINT_ITER_VECTORS )
00421         *out<< "Fd =\n" << *Fd;
00422     }
00423     *out<< "\nEnd output from QPSolverRelaxed::solve_qp(...)\n";
00424   }
00425 }
00426 
00427 } // end namespace ConstrainedOptPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends