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

Generated on Thu Sep 18 12:34:16 2008 for ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization by doxygen 1.3.9.1