Logo ROOT   6.12/07
Reference Guide
TMehrotraSolver.cxx
Go to the documentation of this file.
1 // @(#)root/quadp:$Id$
2 // Author: Eddy Offermann May 2004
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /*************************************************************************
13  * Parts of this file are copied from the OOQP distribution and *
14  * are subject to the following license: *
15  * *
16  * COPYRIGHT 2001 UNIVERSITY OF CHICAGO *
17  * *
18  * The copyright holder hereby grants you royalty-free rights to use, *
19  * reproduce, prepare derivative works, and to redistribute this software*
20  * to others, provided that any changes are clearly documented. This *
21  * software was authored by: *
22  * *
23  * E. MICHAEL GERTZ gertz@mcs.anl.gov *
24  * Mathematics and Computer Science Division *
25  * Argonne National Laboratory *
26  * 9700 S. Cass Avenue *
27  * Argonne, IL 60439-4844 *
28  * *
29  * STEPHEN J. WRIGHT swright@cs.wisc.edu *
30  * Computer Sciences Department *
31  * University of Wisconsin *
32  * 1210 West Dayton Street *
33  * Madison, WI 53706 FAX: (608)262-9777 *
34  * *
35  * Any questions or comments may be directed to one of the authors. *
36  * *
37  * ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES OF *
38  * ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT, AND *
39  * OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A CONTRACT *
40  * WITH THE DEPARTMENT OF ENERGY. *
41  *************************************************************************/
42 
43 //////////////////////////////////////////////////////////////////////////
44 // //
45 // TMehrotraSolver //
46 // //
47 // Derived class of TQpSolverBase implementing the original Mehrotra //
48 // predictor-corrector algorithm //
49 // //
50 //////////////////////////////////////////////////////////////////////////
51 
52 #include "Riostream.h"
53 #include "TMath.h"
54 #include "TMehrotraSolver.h"
55 
57 
58 ////////////////////////////////////////////////////////////////////////////////
59 /// Default constructor
60 
62 {
63  fPrintlevel = 0;
64  fTsig = 0.0;
65  fStep = 0;
66  fFactory = 0;
67 }
68 
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 /// Constructor
72 
74 {
75  fFactory = of;
76  fStep = fFactory->MakeVariables(prob);
77 
78  fPrintlevel = verbose;
79  fTsig = 3.0; // the usual value for the centering exponent (tau)
80 }
81 
82 
83 ////////////////////////////////////////////////////////////////////////////////
84 /// Copy constructor
85 
87 {
88  *this = another;
89 }
90 
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// Solve the quadratic programming problem as formulated through prob, store
94 /// the final solution in iterate->fX . Monitor the residuals during the iterations
95 /// through resid . The status is returned as defined in TQpSolverBase::ETerminationCode .
96 
98 {
99  Int_t status_code;
100  Double_t alpha = 1;
101  Double_t sigma = 1;
102 
103  fDnorm = prob->DataNorm();
104 
105  // initialization of (x,y,z) and factorization routine.
106  fSys = fFactory->MakeLinSys(prob);
107  this->Start(fFactory,iterate,prob,resid,fStep);
108 
109  fIter = 0;
110  Double_t mu = iterate->GetMu();
111 
112  Int_t done = 0;
113  do {
114  fIter++;
115 
116  // evaluate residuals and update algorithm status:
117  resid->CalcResids(prob,iterate);
118 
119  // termination test:
120  status_code = this->DoStatus(prob,iterate,resid,fIter,mu,0);
121  if (status_code != kNOT_FINISHED ) break;
122  if (fPrintlevel >= 10)
123  this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,0);
124 
125  // *** Predictor step ***
126 
127  resid->Set_r3_xz_alpha(iterate,0.0);
128 
129  fSys->Factor(prob,iterate);
130  fSys->Solve(prob,iterate,resid,fStep);
131  fStep->Negate();
132 
133  alpha = iterate->StepBound(fStep);
134 
135  // calculate centering parameter
136  Double_t muaff = iterate->MuStep(fStep,alpha);
137  sigma = TMath::Power(muaff/mu,fTsig);
138 
139  // *** Corrector step ***
140 
141  // form right hand side of linear system:
142  resid->Add_r3_xz_alpha(fStep,-sigma*mu );
143 
144  fSys->Solve(prob,iterate,resid,fStep);
145  fStep->Negate();
146 
147  // We've finally decided on a step direction, now calculate the
148  // length using Mehrotra's heuristic.
149  alpha = this->FinalStepLength(iterate,fStep);
150 
151  // alternatively, just use a crude step scaling factor.
152  //alpha = 0.995 * iterate->StepBound(fStep);
153 
154  // actually take the step and calculate the new mu
155  iterate->Saxpy(fStep,alpha);
156  mu = iterate->GetMu();
157  } while(!done);
158 
159  resid->CalcResids(prob,iterate);
160  if (fPrintlevel >= 10)
161  this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,1);
162 
163  return status_code;
164 }
165 
166 
167 ////////////////////////////////////////////////////////////////////////////////
168 /// Print information about the optimization process and monitor the convergence
169 /// status of thye algorithm
170 
171 void TMehrotraSolver::DefMonitor(TQpDataBase * /* data */,TQpVar * /* vars */,
172  TQpResidual *resids,
173  Double_t alpha,Double_t /* sigma */,Int_t i,Double_t mu,
174  Int_t status_code,Int_t level)
175 {
176  switch (level) {
177  case 0 : case 1:
178  {
179  std::cout << std::endl << "Duality Gap: " << resids->GetDualityGap() << std::endl;
180  if (i > 1) {
181  std::cout << " alpha = " << alpha << std::endl;
182  }
183  std::cout << " *** Iteration " << i << " *** " << std::endl;
184  std::cout << " mu = " << mu << " relative residual norm = "
185  << resids->GetResidualNorm()/fDnorm << std::endl;
186 
187  if (level == 1) {
188  // Termination has been detected by the status check; print
189  // appropriate message
190  switch (status_code) {
192  std::cout << std::endl << " *** SUCCESSFUL TERMINATION ***" << std::endl;
193  break;
194  case kMAX_ITS_EXCEEDED:
195  std::cout << std::endl << " *** MAXIMUM ITERATIONS REACHED *** " << std::endl;
196  break;
197  case kINFEASIBLE:
198  std::cout << std::endl << " *** TERMINATION: PROBABLY INFEASIBLE *** " << std::endl;
199  break;
200  case kUNKNOWN:
201  std::cout << std::endl << " *** TERMINATION: STATUS UNKNOWN *** " << std::endl;
202  break;
203  }
204  }
205  } break; // end case 0: case 1:
206  } // end switch(level)
207 }
208 
209 
210 ////////////////////////////////////////////////////////////////////////////////
211 /// Deconstructor
212 
214 {
215  delete fStep;
216 }
217 
218 
219 ////////////////////////////////////////////////////////////////////////////////
220 /// Assignment operator
221 
223 {
224  if (this != &source) {
225  TQpSolverBase::operator=(source);
226 
227  fPrintlevel = source.fPrintlevel;
228  fTsig = source.fTsig;
229 
230  if (fStep) delete fStep;
231 
232  fStep = new TQpVar(*source.fStep);
233  fFactory = source.fFactory;
234  }
235  return *this;
236 }
virtual Double_t DataNorm()=0
int iterate(rng_state_t *X)
Definition: mixmax.icc:34
virtual TQpVar * MakeVariables(const TQpDataBase *data)=0
virtual void DoMonitor(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Double_t alpha, Double_t sigma, Int_t i, Double_t mu, Int_t stop_code, Int_t level)
Monitor progress / convergence aat each interior-point iteration.
void CalcResids(TQpDataBase *problem, TQpVar *vars)
Calculate residuals, their norms, and duality complementarity gap, given a problem and variable set...
TMehrotraSolver & operator=(const TMehrotraSolver &source)
Assignment operator.
TMehrotraSolver()
Default constructor.
virtual Double_t MuStep(TQpVar *step, Double_t alpha)
Compute the complementarity gap resulting from a step of length "alpha" along direction "step"...
Definition: TQpVar.cxx:210
int Int_t
Definition: RtypesCore.h:41
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:627
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
virtual Double_t FinalStepLength(TQpVar *iterate, TQpVar *step)
Implements a version of Mehrotra starting point heuristic, modified to ensure identical steps in the ...
virtual Int_t DoStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Tests for termination.
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
virtual Double_t GetMu()
compute complementarity gap, obtained by taking the inner product of the complementary vectors and di...
Definition: TQpVar.cxx:191
static Double_t StepBound(TVectorD &v, TVectorD &dir, Double_t maxStep)
Find the maximum stepsize of v in direction dir before violating the nonnegativity constraints...
Definition: TQpVar.cxx:347
const Double_t sigma
virtual void Start(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Implements a default starting-point heuristic.
Double_t GetResidualNorm()
Definition: TQpResidual.h:108
void Add_r3_xz_alpha(TQpVar *vars, Double_t alpha)
Modify the "complementarity" component of the residuals, by adding the pairwise products of the compl...
Double_t fDnorm
Definition: TQpSolverBase.h:78
TQpProbBase * fFactory
Double_t GetDualityGap()
Definition: TQpResidual.h:109
TQpLinSolverBase * fSys
Definition: TQpSolverBase.h:76
virtual TQpLinSolverBase * MakeLinSys(const TQpDataBase *data)=0
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition: TQpVar.cxx:271
void Set_r3_xz_alpha(TQpVar *vars, Double_t alpha)
Set the "complementarity" component of the residuals to the pairwise products of the complementary va...
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
virtual void Saxpy(TQpVar *b, Double_t alpha)
Perform a "saxpy" operation on all data vectors : x += alpha*y.
Definition: TQpVar.cxx:231
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
virtual void DefMonitor(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Double_t alpha, Double_t sigma, Int_t i, Double_t mu, Int_t status_code, Int_t level)
Print information about the optimization process and monitor the convergence status of thye algorithm...
Definition: TQpVar.h:59
virtual ~TMehrotraSolver()
Deconstructor.
virtual Int_t Solve(TQpDataBase *prob, TQpVar *iterate, TQpResidual *resid)
Solve the quadratic programming problem as formulated through prob, store the final solution in itera...