Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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/// \class TMehrotraSolver
46///
47/// Derived class of TQpSolverBase implementing the original Mehrotra
48/// predictor-corrector algorithm
49///
50////////////////////////////////////////////////////////////////////////////////
51
52#include <iostream>
53#include "TMath.h"
54#include "TMehrotraSolver.h"
55
56
57////////////////////////////////////////////////////////////////////////////////
58/// Default constructor
59
61{
62 fPrintlevel = 0;
63 fTsig = 0.0;
64 fStep = nullptr;
65 fFactory = nullptr;
66}
67
68
69////////////////////////////////////////////////////////////////////////////////
70/// Constructor
71
73{
74 fFactory = of;
76
77 fPrintlevel = verbose;
78 fTsig = 3.0; // the usual value for the centering exponent (tau)
79}
80
81
82////////////////////////////////////////////////////////////////////////////////
83/// Copy constructor
84
89
90
91////////////////////////////////////////////////////////////////////////////////
92/// Solve the quadratic programming problem as formulated through prob, store
93/// the final solution in iterate->fX . Monitor the residuals during the iterations
94/// through resid . The status is returned as defined in TQpSolverBase::ETerminationCode .
95
97{
98 Int_t status_code;
99 Double_t alpha = 1;
100 Double_t sigma = 1;
101
102 fDnorm = prob->DataNorm();
103
104 // initialization of (x,y,z) and factorization routine.
106 this->Start(fFactory,iterate,prob,resid,fStep);
107
108 fIter = 0;
109 Double_t mu = iterate->GetMu();
110
111 Int_t done = 0;
112 do {
113 fIter++;
114
115 // evaluate residuals and update algorithm status:
116 resid->CalcResids(prob,iterate);
117
118 // termination test:
119 status_code = this->DoStatus(prob,iterate,resid,fIter,mu,0);
120 if (status_code != kNOT_FINISHED ) break;
121 if (fPrintlevel >= 10)
122 this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,0);
123
124 // *** Predictor step ***
125
126 resid->Set_r3_xz_alpha(iterate,0.0);
127
130 fStep->Negate();
131
132 alpha = iterate->StepBound(fStep);
133
134 // calculate centering parameter
135 Double_t muaff = iterate->MuStep(fStep,alpha);
137
138 // *** Corrector step ***
139
140 // form right hand side of linear system:
141 resid->Add_r3_xz_alpha(fStep,-sigma*mu );
142
144 fStep->Negate();
145
146 // We've finally decided on a step direction, now calculate the
147 // length using Mehrotra's heuristic.
148 alpha = this->FinalStepLength(iterate,fStep);
149
150 // alternatively, just use a crude step scaling factor.
151 //alpha = 0.995 * iterate->StepBound(fStep);
152
153 // actually take the step and calculate the new mu
154 iterate->Saxpy(fStep,alpha);
155 mu = iterate->GetMu();
156 } while(!done);
157
158 resid->CalcResids(prob,iterate);
159 if (fPrintlevel >= 10)
160 this->DoMonitor(prob,iterate,resid,alpha,sigma,fIter,mu,status_code,1);
161
162 return status_code;
163}
164
165
166////////////////////////////////////////////////////////////////////////////////
167/// Print information about the optimization process and monitor the convergence
168/// status of thye algorithm
169
170void TMehrotraSolver::DefMonitor(TQpDataBase * /* data */,TQpVar * /* vars */,
172 Double_t alpha,Double_t /* sigma */,Int_t i,Double_t mu,
173 Int_t status_code,Int_t level)
174{
175 switch (level) {
176 case 0 : case 1:
177 {
178 std::cout << std::endl << "Duality Gap: " << resids->GetDualityGap() << std::endl;
179 if (i > 1) {
180 std::cout << " alpha = " << alpha << std::endl;
181 }
182 std::cout << " *** Iteration " << i << " *** " << std::endl;
183 std::cout << " mu = " << mu << " relative residual norm = "
184 << resids->GetResidualNorm()/fDnorm << std::endl;
185
186 if (level == 1) {
187 // Termination has been detected by the status check; print
188 // appropriate message
189 switch (status_code) {
191 std::cout << std::endl << " *** SUCCESSFUL TERMINATION ***" << std::endl;
192 break;
194 std::cout << std::endl << " *** MAXIMUM ITERATIONS REACHED *** " << std::endl;
195 break;
196 case kINFEASIBLE:
197 std::cout << std::endl << " *** TERMINATION: PROBABLY INFEASIBLE *** " << std::endl;
198 break;
199 case kUNKNOWN:
200 std::cout << std::endl << " *** TERMINATION: STATUS UNKNOWN *** " << std::endl;
201 break;
202 }
203 }
204 } break; // end case 0: case 1:
205 } // end switch(level)
206}
207
208
209////////////////////////////////////////////////////////////////////////////////
210/// Deconstructor
211
216
217
218////////////////////////////////////////////////////////////////////////////////
219/// Assignment operator
220
222{
223 if (this != &source) {
225
226 fPrintlevel = source.fPrintlevel;
227 fTsig = source.fTsig;
228
229 if (fStep) delete fStep;
230
231 fStep = new TQpVar(*source.fStep);
232 fFactory = source.fFactory;
233 }
234 return *this;
235}
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
@ kNOT_FINISHED
@ kINFEASIBLE
@ kMAX_ITS_EXCEEDED
@ kUNKNOWN
@ kSUCCESSFUL_TERMINATION
Derived class of TQpSolverBase implementing the original Mehrotra predictor-corrector algorithm.
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) override
Print information about the optimization process and monitor the convergence status of thye algorithm...
TMehrotraSolver()
Default constructor.
TMehrotraSolver & operator=(const TMehrotraSolver &source)
Assignment operator.
~TMehrotraSolver() override
Deconstructor.
Int_t Solve(TQpDataBase *prob, TQpVar *iterate, TQpResidual *resid) override
Solve the quadratic programming problem as formulated through prob, store the final solution in itera...
TQpProbBase * fFactory
Data for the general QP formulation.
Definition TQpDataBase.h:61
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
default general problem formulation:
Definition TQpProbBase.h:89
virtual TQpLinSolverBase * MakeLinSys(const TQpDataBase *data)=0
virtual TQpVar * MakeVariables(const TQpDataBase *data)=0
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition TQpResidual.h:62
The Solver class contains methods for monitoring and checking the convergence status of the algorithm...
TQpLinSolverBase * fSys
Double_t fDnorm
virtual void Start(TQpProbBase *formulation, TQpVar *iterate, TQpDataBase *prob, TQpResidual *resid, TQpVar *step)
Implements a default starting-point heuristic.
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.
virtual Int_t DoStatus(TQpDataBase *data, TQpVar *vars, TQpResidual *resids, Int_t i, Double_t mu, Int_t level)
Tests for termination.
virtual Double_t FinalStepLength(TQpVar *iterate, TQpVar *step)
Implements a version of Mehrotra starting point heuristic, modified to ensure identical steps in the ...
TQpSolverBase & operator=(const TQpSolverBase &source)
Assignment operator.
Class containing the variables for the general QP formulation.
Definition TQpVar.h:60
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition TQpVar.cxx:271
const Double_t sigma
int iterate(rng_state_t *X)
Definition mixmax.icc:34
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732