Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 "TQpResidual.h"
57
58
59////////////////////////////////////////////////////////////////////////////////
60/// Constructor
61
63{
64 fNx = 0;
65 fMy = 0;
66 fMz = 0;
67
68 fNxup = 0.0;
69 fNxlo = 0.0;
70 fMcup = 0.0;
71 fMclo = 0.0;
72 fResidualNorm = 0.0;
73 fDualityGap = 0.0;
74}
75
76
77////////////////////////////////////////////////////////////////////////////////
78/// Constructor
79
82{
83 fNx = nx;
84 fMy = my;
85 fMz = mz;
86
87 if (ixlo.GetNrows() > 0) fXloIndex.Use(ixlo.GetNrows(),ixlo.GetMatrixArray());
88 if (ixup.GetNrows() > 0) fXupIndex.Use(ixup.GetNrows(),ixup.GetMatrixArray());
89 if (iclo.GetNrows() > 0) fCloIndex.Use(iclo.GetNrows(),iclo.GetMatrixArray());
90 if (icup.GetNrows() > 0) fCupIndex.Use(icup.GetNrows(),icup.GetMatrixArray());
91 fNxlo = ixlo.NonZeros();
92 fNxup = ixup.NonZeros();
93 fMclo = iclo.NonZeros();
94 fMcup = icup.NonZeros();
95
99
101 if (fMclo > 0) {
104 }
105 if (fMcup > 0) {
108 }
109 if (fNxlo > 0) {
112 }
113 if (fNxup > 0) {
116 }
117
118 fResidualNorm = 0.0;
119 fDualityGap = 0.0;
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Copy constructor
125
130
131
132////////////////////////////////////////////////////////////////////////////////
133/// Calculate residuals, their norms, and duality complementarity gap,
134/// given a problem and variable set.
135
137{
139
140 fRQ.ResizeTo(prob->fG); fRQ = prob->fG;
141 prob->Qmult(1.0,fRQ,1.0,vars->fX);
142
143 // calculate x^T (g+Qx) - contribution to the duality gap
144 Double_t gap = fRQ*vars->fX;
145
146 prob->ATransmult(1.0,fRQ,-1.0,vars->fY);
147 prob->CTransmult(1.0,fRQ,-1.0,vars->fZ);
148 if (fNxlo > 0) Add(fRQ,-1.0,vars->fGamma);
149 if (fNxup > 0) Add(fRQ, 1.0,vars->fPhi);
150
151 Double_t norm = 0.0;
154
155 fRA.ResizeTo(prob->fBa); fRA = prob->fBa;
156 prob->Amult(-1.0,fRA,1.0,vars->fX);
157
158 // contribution -d^T y to duality gap
159 gap -= prob->fBa*vars->fY;
160
163
164 fRC.ResizeTo(vars->fS); fRC = vars->fS;
165 prob->Cmult(-1.0,fRC,1.0,vars->fX);
166
169
170 fRz.ResizeTo(vars->fZ); fRz = vars->fZ;
171
172 if (fMclo > 0) {
173 Add(fRz,-1.0,vars->fLambda);
174
175 fRt.ResizeTo(vars->fS); fRt = vars->fS;
176 Add(fRt,-1.0,prob->GetSlowerBound());
178 Add(fRt,-1.0,vars->fT);
179
180 gap -= prob->fCloBound*vars->fLambda;
181
184 }
185
186 if (fMcup > 0) {
187 Add(fRz,1.0,vars->fPi);
188
189 fRu.ResizeTo(vars->fS); fRu = vars->fS;
190 Add(fRu,-1.0,prob->GetSupperBound() );
192 Add(fRu,1.0,vars->fU);
193
194 gap += prob->fCupBound*vars->fPi;
195
198 }
199
202
203 if (fNxlo > 0) {
204 fRv.ResizeTo(vars->fX); fRv = vars->fX;
205 Add(fRv,-1.0,prob->GetXlowerBound());
207 Add(fRv,-1.0,vars->fV);
208
209 gap -= prob->fXloBound*vars->fGamma;
210
213 }
214
215 if (fNxup > 0) {
216 fRw.ResizeTo(vars->fX); fRw = vars->fX;
217 Add(fRw,-1.0,prob->GetXupperBound());
219 Add(fRw,1.0,vars->fW);
220
221 gap += prob->fXupBound*vars->fPhi;
222
225 }
226
227 fDualityGap = gap;
229}
230
231
232////////////////////////////////////////////////////////////////////////////////
233/// Modify the "complementarity" component of the residuals, by adding the pairwise
234/// products of the complementary variables plus a constant alpha to this term.
235
237{
238 if (fMclo > 0) AddElemMult(fRlambda,1.0,vars->fT,vars->fLambda);
239 if (fMcup > 0) AddElemMult(fRpi ,1.0,vars->fU,vars->fPi);
240 if (fNxlo > 0) AddElemMult(fRgamma ,1.0,vars->fV,vars->fGamma);
241 if (fNxup > 0) AddElemMult(fRphi ,1.0,vars->fW,vars->fPhi);
242
243 if (alpha != 0.0) {
244 if (fMclo > 0) fRlambda.AddSomeConstant(alpha,fCloIndex);
245 if (fMcup > 0) fRpi .AddSomeConstant(alpha,fCupIndex);
246 if (fNxlo > 0) fRgamma .AddSomeConstant(alpha,fXloIndex);
247 if (fNxup > 0) fRphi .AddSomeConstant(alpha,fXupIndex);
248 }
249}
250
251
252////////////////////////////////////////////////////////////////////////////////
253/// Set the "complementarity" component of the residuals to the pairwise products of
254/// the complementary variables plus a constant alpha .
255
257{
258 this->Clear_r3();
259 this->Add_r3_xz_alpha(vars,alpha);
260}
261
262
263////////////////////////////////////////////////////////////////////////////////
264/// set the complementarity component of the residuals to 0.
265
267{
268 if (fMclo > 0) fRlambda.Zero();
269 if (fMcup > 0) fRpi .Zero();
270 if (fNxlo > 0) fRgamma .Zero();
271 if (fNxup > 0) fRphi .Zero();
272}
273
274
275////////////////////////////////////////////////////////////////////////////////
276/// set the noncomplementarity components of the residual (the terms arising from
277/// the linear equalities in the KKT conditions) to 0.
278
280{
281 fRQ.Zero();
282 fRA.Zero();
283 fRC.Zero();
284 fRz.Zero();
285 if (fNxlo > 0) fRv.Zero();
286 if (fNxup > 0) fRw.Zero();
287 if (fMclo > 0) fRt.Zero();
288 if (fMcup > 0) fRu.Zero();
289}
290
291
292////////////////////////////////////////////////////////////////////////////////
293/// Perform the projection operation required by Gondzio algorithm: replace each
294/// component r3_i of the complementarity component of the residuals by r3p_i-r3_i,
295/// where r3p_i is the projection of r3_i onto the box [rmin, rmax]. Then if the
296/// resulting value is less than -rmax, replace it by -rmax.
297
317
318
319////////////////////////////////////////////////////////////////////////////////
320/// Check if vector elements as selected through array indices are non-zero
321
323{
324 if (fNxlo > 0 &&
325 (!fRv .MatchesNonZeroPattern(fXloIndex) ||
327 return kFALSE;
328 }
329
330 if (fNxup > 0 &&
331 (!fRw .MatchesNonZeroPattern(fXupIndex) ||
333 return kFALSE;
334 }
335 if (fMclo > 0 &&
336 (!fRt .MatchesNonZeroPattern(fCloIndex) ||
338 return kFALSE;
339 }
340
341 if (fMcup > 0 &&
342 (!fRu .MatchesNonZeroPattern(fCupIndex) ||
344 return kFALSE;
345 }
346
347 return kTRUE;
348}
349
350
351////////////////////////////////////////////////////////////////////////////////
352/// Replace each component r3_i of the complementarity component of the residuals
353/// by r3p_i-r3_i, where r3p_i is the projection of r3_i onto the box [rmin, rmax].
354/// Then if the resulting value is less than -rmax, replace it by -rmax.
355
357{
358 Double_t * ep = v.GetMatrixArray();
359 const Double_t * const fep = ep+v.GetNrows();
360
361 while (ep < fep) {
362 if (*ep < rmin)
363 *ep = rmin - *ep;
364 else if (*ep > rmax)
365 *ep = rmax - *ep;
366 else
367 *ep = 0.0;
368 if (*ep < -rmax) *ep = -rmax;
369 ep++;
370 }
371}
372
373
374////////////////////////////////////////////////////////////////////////////////
375/// Assignment operator
376
378{
379 if (this != &source) {
381
382 fNx = source.fNx;
383 fMy = source.fMy;
384 fMz = source.fMz;
385
386 fNxup = source.fNxup;
387 fNxlo = source.fNxlo;
388 fMcup = source.fMcup;
389 fMclo = source.fMclo;
390
391 fXupIndex.ResizeTo(source.fXupIndex); fXupIndex = source.fXupIndex;
392 fXloIndex.ResizeTo(source.fXloIndex); fXloIndex = source.fXloIndex;
393 fCupIndex.ResizeTo(source.fCupIndex); fCupIndex = source.fCupIndex;
394 fCloIndex.ResizeTo(source.fCloIndex); fCupIndex = source.fCupIndex;
395
396 fRQ .ResizeTo(source.fRQ); fRQ = source.fRQ;
397 fRA .ResizeTo(source.fRA); fRA = source.fRA;
398 fRC .ResizeTo(source.fRC); fRC = source.fRC;
399 fRz .ResizeTo(source.fRz); fRz = source.fRz;
400 fRv .ResizeTo(source.fRv); fRv = source.fRv;
401 fRw .ResizeTo(source.fRw); fRw = source.fRw;
402 fRt .ResizeTo(source.fRt); fRt = source.fRt;
403 fRu .ResizeTo(source.fRu); fRu = source.fRu;
404 fRgamma .ResizeTo(source.fRgamma); fRgamma = source.fRgamma;
405 fRphi .ResizeTo(source.fRphi); fRphi = source.fRphi;
406 fRlambda.ResizeTo(source.fRlambda); fRlambda = source.fRlambda;
407 fRpi .ResizeTo(source.fRpi); fRpi = source.fRpi;
408
409 // LM: copy also these data members
410 fResidualNorm = source.fResidualNorm;
411 fDualityGap = source.fDualityGap;
412 }
413 return *this;
414}
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Mother of all ROOT objects.
Definition TObject.h:41
TObject & operator=(const TObject &rhs) noexcept
TObject assignment operator.
Definition TObject.h:299
Data for the general QP formulation.
Definition TQpDataBase.h:61
Data for the dense QP formulation.
Definition TQpDataDens.h:63
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.
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:452
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition TVectorT.cxx:602
void AddSomeConstant(Element val, const TVectorT< Element > &select)
Add to vector elements as selected through array select the value val.
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:293
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition TVectorT.cxx:348
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition TVectorT.cxx:543
TVectorT< Element > & AddElemMult(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementMult(source1,source2) .
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.