Logo ROOT  
Reference Guide
TQpResidual.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 TQpResidual
46///
47/// The Residuals class calculates and stores the quantities that appear
48/// on the right-hand side of the linear systems that arise at each
49/// interior-point iteration. These residuals can be partitioned into
50/// two fundamental categories: the components arising from the linear
51/// equations in the KKT conditions, and the components arising from the
52/// complementarity conditions.
53///
54////////////////////////////////////////////////////////////////////////////////
55
56#include "Riostream.h"
57#include "TQpResidual.h"
58#include "TMatrixD.h"
59
61
62////////////////////////////////////////////////////////////////////////////////
63/// Constructor
64
66{
67 fNx = 0;
68 fMy = 0;
69 fMz = 0;
70
71 fNxup = 0.0;
72 fNxlo = 0.0;
73 fMcup = 0.0;
74 fMclo = 0.0;
75 fResidualNorm = 0.0;
76 fDualityGap = 0.0;
77}
78
79
80////////////////////////////////////////////////////////////////////////////////
81/// Constructor
82
84 TVectorD &iclo,TVectorD &icup)
85{
86 fNx = nx;
87 fMy = my;
88 fMz = mz;
89
90 if (ixlo.GetNrows() > 0) fXloIndex.Use(ixlo.GetNrows(),ixlo.GetMatrixArray());
91 if (ixup.GetNrows() > 0) fXupIndex.Use(ixup.GetNrows(),ixup.GetMatrixArray());
92 if (iclo.GetNrows() > 0) fCloIndex.Use(iclo.GetNrows(),iclo.GetMatrixArray());
93 if (icup.GetNrows() > 0) fCupIndex.Use(icup.GetNrows(),icup.GetMatrixArray());
94 fNxlo = ixlo.NonZeros();
95 fNxup = ixup.NonZeros();
96 fMclo = iclo.NonZeros();
97 fMcup = icup.NonZeros();
98
102
104 if (fMclo > 0) {
107 }
108 if (fMcup > 0) {
111 }
112 if (fNxlo > 0) {
115 }
116 if (fNxup > 0) {
119 }
120
121 fResidualNorm = 0.0;
122 fDualityGap = 0.0;
123}
124
125
126////////////////////////////////////////////////////////////////////////////////
127/// Copy constructor
128
130{
131 *this = another;
132}
133
134
135////////////////////////////////////////////////////////////////////////////////
136/// Calculate residuals, their norms, and duality complementarity gap,
137/// given a problem and variable set.
138
140{
141 TQpDataDens *prob = (TQpDataDens *) prob_in;
142
143 fRQ.ResizeTo(prob->fG); fRQ = prob->fG;
144 prob->Qmult(1.0,fRQ,1.0,vars->fX);
145
146 // calculate x^T (g+Qx) - contribution to the duality gap
147 Double_t gap = fRQ*vars->fX;
148
149 prob->ATransmult(1.0,fRQ,-1.0,vars->fY);
150 prob->CTransmult(1.0,fRQ,-1.0,vars->fZ);
151 if (fNxlo > 0) Add(fRQ,-1.0,vars->fGamma);
152 if (fNxup > 0) Add(fRQ, 1.0,vars->fPhi);
153
154 Double_t norm = 0.0;
155 Double_t componentNorm = fRQ.NormInf();
156 if (componentNorm > norm) norm = componentNorm;
157
158 fRA.ResizeTo(prob->fBa); fRA = prob->fBa;
159 prob->Amult(-1.0,fRA,1.0,vars->fX);
160
161 // contribution -d^T y to duality gap
162 gap -= prob->fBa*vars->fY;
163
164 componentNorm = fRA.NormInf();
165 if( componentNorm > norm ) norm = componentNorm;
166
167 fRC.ResizeTo(vars->fS); fRC = vars->fS;
168 prob->Cmult(-1.0,fRC,1.0,vars->fX);
169
170 componentNorm = fRC.NormInf();
171 if( componentNorm > norm ) norm = componentNorm;
172
173 fRz.ResizeTo(vars->fZ); fRz = vars->fZ;
174
175 if (fMclo > 0) {
176 Add(fRz,-1.0,vars->fLambda);
177
178 fRt.ResizeTo(vars->fS); fRt = vars->fS;
179 Add(fRt,-1.0,prob->GetSlowerBound());
181 Add(fRt,-1.0,vars->fT);
182
183 gap -= prob->fCloBound*vars->fLambda;
184
185 componentNorm = fRt.NormInf();
186 if( componentNorm > norm ) norm = componentNorm;
187 }
188
189 if (fMcup > 0) {
190 Add(fRz,1.0,vars->fPi);
191
192 fRu.ResizeTo(vars->fS); fRu = vars->fS;
193 Add(fRu,-1.0,prob->GetSupperBound() );
195 Add(fRu,1.0,vars->fU);
196
197 gap += prob->fCupBound*vars->fPi;
198
199 componentNorm = fRu.NormInf();
200 if( componentNorm > norm ) norm = componentNorm;
201 }
202
203 componentNorm = fRz.NormInf();
204 if( componentNorm > norm ) norm = componentNorm;
205
206 if (fNxlo > 0) {
207 fRv.ResizeTo(vars->fX); fRv = vars->fX;
208 Add(fRv,-1.0,prob->GetXlowerBound());
210 Add(fRv,-1.0,vars->fV);
211
212 gap -= prob->fXloBound*vars->fGamma;
213
214 componentNorm = fRv.NormInf();
215 if( componentNorm > norm ) norm = componentNorm;
216 }
217
218 if (fNxup > 0) {
219 fRw.ResizeTo(vars->fX); fRw = vars->fX;
220 Add(fRw,-1.0,prob->GetXupperBound());
222 Add(fRw,1.0,vars->fW);
223
224 gap += prob->fXupBound*vars->fPhi;
225
226 componentNorm = fRw.NormInf();
227 if (componentNorm > norm) norm = componentNorm;
228 }
229
230 fDualityGap = gap;
231 fResidualNorm = norm;
232}
233
234
235////////////////////////////////////////////////////////////////////////////////
236/// Modify the "complementarity" component of the residuals, by adding the pairwise
237/// products of the complementary variables plus a constant alpha to this term.
238
240{
241 if (fMclo > 0) AddElemMult(fRlambda,1.0,vars->fT,vars->fLambda);
242 if (fMcup > 0) AddElemMult(fRpi ,1.0,vars->fU,vars->fPi);
243 if (fNxlo > 0) AddElemMult(fRgamma ,1.0,vars->fV,vars->fGamma);
244 if (fNxup > 0) AddElemMult(fRphi ,1.0,vars->fW,vars->fPhi);
245
246 if (alpha != 0.0) {
247 if (fMclo > 0) fRlambda.AddSomeConstant(alpha,fCloIndex);
248 if (fMcup > 0) fRpi .AddSomeConstant(alpha,fCupIndex);
249 if (fNxlo > 0) fRgamma .AddSomeConstant(alpha,fXloIndex);
250 if (fNxup > 0) fRphi .AddSomeConstant(alpha,fXupIndex);
251 }
252}
253
254
255////////////////////////////////////////////////////////////////////////////////
256/// Set the "complementarity" component of the residuals to the pairwise products of
257/// the complementary variables plus a constant alpha .
258
260{
261 this->Clear_r3();
262 this->Add_r3_xz_alpha(vars,alpha);
263}
264
265
266////////////////////////////////////////////////////////////////////////////////
267/// set the complementarity component of the residuals to 0.
268
270{
271 if (fMclo > 0) fRlambda.Zero();
272 if (fMcup > 0) fRpi .Zero();
273 if (fNxlo > 0) fRgamma .Zero();
274 if (fNxup > 0) fRphi .Zero();
275}
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// set the noncomplementarity components of the residual (the terms arising from
280/// the linear equalities in the KKT conditions) to 0.
281
283{
284 fRQ.Zero();
285 fRA.Zero();
286 fRC.Zero();
287 fRz.Zero();
288 if (fNxlo > 0) fRv.Zero();
289 if (fNxup > 0) fRw.Zero();
290 if (fMclo > 0) fRt.Zero();
291 if (fMcup > 0) fRu.Zero();
292}
293
294
295////////////////////////////////////////////////////////////////////////////////
296/// Perform the projection operation required by Gondzio algorithm: replace each
297/// component r3_i of the complementarity component of the residuals by r3p_i-r3_i,
298/// where r3p_i is the projection of r3_i onto the box [rmin, rmax]. Then if the
299/// resulting value is less than -rmax, replace it by -rmax.
300
302{
303 if (fMclo > 0) {
304 GondzioProjection(fRlambda,rmin,rmax);
306 }
307 if (fMcup > 0) {
308 GondzioProjection(fRpi,rmin,rmax);
310 }
311 if (fNxlo > 0) {
312 GondzioProjection(fRgamma,rmin,rmax);
314 }
315 if (fNxup > 0) {
316 GondzioProjection(fRphi,rmin,rmax);
318 }
319}
320
321
322////////////////////////////////////////////////////////////////////////////////
323/// Check if vector elements as selected through array indices are non-zero
324
326{
327 if (fNxlo > 0 &&
328 (!fRv .MatchesNonZeroPattern(fXloIndex) ||
330 return kFALSE;
331 }
332
333 if (fNxup > 0 &&
334 (!fRw .MatchesNonZeroPattern(fXupIndex) ||
336 return kFALSE;
337 }
338 if (fMclo > 0 &&
339 (!fRt .MatchesNonZeroPattern(fCloIndex) ||
341 return kFALSE;
342 }
343
344 if (fMcup > 0 &&
345 (!fRu .MatchesNonZeroPattern(fCupIndex) ||
347 return kFALSE;
348 }
349
350 return kTRUE;
351}
352
353
354////////////////////////////////////////////////////////////////////////////////
355/// Replace each component r3_i of the complementarity component of the residuals
356/// by r3p_i-r3_i, where r3p_i is the projection of r3_i onto the box [rmin, rmax].
357/// Then if the resulting value is less than -rmax, replace it by -rmax.
358
360{
361 Double_t * ep = v.GetMatrixArray();
362 const Double_t * const fep = ep+v.GetNrows();
363
364 while (ep < fep) {
365 if (*ep < rmin)
366 *ep = rmin - *ep;
367 else if (*ep > rmax)
368 *ep = rmax - *ep;
369 else
370 *ep = 0.0;
371 if (*ep < -rmax) *ep = -rmax;
372 ep++;
373 }
374}
375
376
377////////////////////////////////////////////////////////////////////////////////
378/// Assignment operator
379
381{
382 if (this != &source) {
383 TObject::operator=(source);
384
385 fNx = source.fNx;
386 fMy = source.fMy;
387 fMz = source.fMz;
388
389 fNxup = source.fNxup;
390 fNxlo = source.fNxlo;
391 fMcup = source.fMcup;
392 fMclo = source.fMclo;
393
398
399 fRQ .ResizeTo(source.fRQ); fRQ = source.fRQ;
400 fRA .ResizeTo(source.fRA); fRA = source.fRA;
401 fRC .ResizeTo(source.fRC); fRC = source.fRC;
402 fRz .ResizeTo(source.fRz); fRz = source.fRz;
403 fRv .ResizeTo(source.fRv); fRv = source.fRv;
404 fRw .ResizeTo(source.fRw); fRw = source.fRw;
405 fRt .ResizeTo(source.fRt); fRt = source.fRt;
406 fRu .ResizeTo(source.fRu); fRu = source.fRu;
407 fRgamma .ResizeTo(source.fRgamma); fRgamma = source.fRgamma;
408 fRphi .ResizeTo(source.fRphi); fRphi = source.fRphi;
409 fRlambda.ResizeTo(source.fRlambda); fRlambda = source.fRlambda;
410 fRpi .ResizeTo(source.fRpi); fRpi = source.fRpi;
411
412 // LM: copy also these data members
414 fDualityGap = source.fDualityGap;
415 }
416 return *this;
417}
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
TVectorT< Element > & AddElemMult(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementMult(source1,source2) .
Definition: TVectorT.cxx:1845
Mother of all ROOT objects.
Definition: TObject.h:37
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:283
Data for the general QP formulation.
Definition: TQpDataBase.h:61
TVectorD fXloBound
Definition: TQpDataBase.h:81
virtual TVectorD & GetSlowerBound()
Definition: TQpDataBase.h:114
TVectorD fG
Definition: TQpDataBase.h:77
TVectorD fXupBound
Definition: TQpDataBase.h:79
virtual TVectorD & GetXlowerBound()
Definition: TQpDataBase.h:110
TVectorD fCupBound
Definition: TQpDataBase.h:83
virtual TVectorD & GetSupperBound()
Definition: TQpDataBase.h:112
TVectorD fCloBound
Definition: TQpDataBase.h:85
TVectorD fBa
Definition: TQpDataBase.h:78
virtual TVectorD & GetXupperBound()
Definition: TQpDataBase.h:108
Data for the dense QP formulation.
Definition: TQpDataDens.h:63
virtual void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC^T*x)
virtual void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA*x)
virtual void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC*x)
virtual void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA^T*x)
virtual void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fQ*x)
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition: TQpResidual.h:62
void Clear_r1r2()
set the noncomplementarity components of the residual (the terms arising from the linear equalities i...
Double_t fDualityGap
Definition: TQpResidual.h:66
Double_t fNxup
Definition: TQpResidual.h:74
void Clear_r3()
set the complementarity component of the residuals to 0.
Double_t fResidualNorm
Definition: TQpResidual.h:65
TVectorD fRgamma
Definition: TQpResidual.h:96
Double_t fMcup
Definition: TQpResidual.h:76
TVectorD fXupIndex
Definition: TQpResidual.h:80
TQpResidual()
Constructor.
Definition: TQpResidual.cxx:65
TVectorD fRu
Definition: TQpResidual.h:95
Double_t fMclo
Definition: TQpResidual.h:77
TVectorD fRpi
Definition: TQpResidual.h:99
TVectorD fXloIndex
Definition: TQpResidual.h:81
TVectorD fRC
Definition: TQpResidual.h:90
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...
Bool_t ValidNonZeroPattern()
Check if vector elements as selected through array indices are non-zero.
TVectorD fCloIndex
Definition: TQpResidual.h:83
TVectorD fRQ
Definition: TQpResidual.h:88
TVectorD fRt
Definition: TQpResidual.h:94
TVectorD fRw
Definition: TQpResidual.h:93
void CalcResids(TQpDataBase *problem, TQpVar *vars)
Calculate residuals, their norms, and duality complementarity gap, given a problem and variable set.
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...
TVectorD fRlambda
Definition: TQpResidual.h:98
TVectorD fRA
Definition: TQpResidual.h:89
TVectorD fRv
Definition: TQpResidual.h:92
Double_t fNxlo
Definition: TQpResidual.h:75
TVectorD fRphi
Definition: TQpResidual.h:97
TQpResidual & operator=(const TQpResidual &source)
Assignment operator.
void Project_r3(Double_t rmin, Double_t rmax)
Perform the projection operation required by Gondzio algorithm: replace each component r3_i of the co...
TVectorD fCupIndex
Definition: TQpResidual.h:82
static void GondzioProjection(TVectorD &v, Double_t rmin, Double_t rmax)
Replace each component r3_i of the complementarity component of the residuals by r3p_i-r3_i,...
TVectorD fRz
Definition: TQpResidual.h:91
Class containing the variables for the general QP formulation.
Definition: TQpVar.h:60
TVectorD fU
Definition: TQpVar.h:105
TVectorD fX
Definition: TQpVar.h:91
TVectorD fLambda
Definition: TQpVar.h:103
TVectorD fGamma
Definition: TQpVar.h:100
TVectorD fT
Definition: TQpVar.h:102
TVectorD fPi
Definition: TQpVar.h:106
TVectorD fS
Definition: TQpVar.h:92
TVectorD fPhi
Definition: TQpVar.h:97
TVectorD fV
Definition: TQpVar.h:96
TVectorD fW
Definition: TQpVar.h:99
TVectorD fY
Definition: TQpVar.h:93
TVectorD fZ
Definition: TQpVar.h:94
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:454
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:604
void AddSomeConstant(Element val, const TVectorT< Element > &select)
Add to vector elements as selected through array select the value val.
Definition: TVectorT.cxx:1285
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1239
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:295
Int_t GetNrows() const
Definition: TVectorT.h:75
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:350
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:545
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:621
Element * GetMatrixArray()
Definition: TVectorT.h:78
void Add(RHist< DIMENSIONS, PRECISION, STAT_TO... > &to, const RHist< DIMENSIONS, PRECISION, STAT_FROM... > &from)
Add two histograms.
Definition: RHist.hxx:323