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
59
60////////////////////////////////////////////////////////////////////////////////
61/// Constructor
62
64{
65 fNx = 0;
66 fMy = 0;
67 fMz = 0;
68
69 fNxup = 0.0;
70 fNxlo = 0.0;
71 fMcup = 0.0;
72 fMclo = 0.0;
73 fResidualNorm = 0.0;
74 fDualityGap = 0.0;
75}
76
77
78////////////////////////////////////////////////////////////////////////////////
79/// Constructor
80
82 TVectorD &iclo,TVectorD &icup)
83{
84 fNx = nx;
85 fMy = my;
86 fMz = mz;
87
88 if (ixlo.GetNrows() > 0) fXloIndex.Use(ixlo.GetNrows(),ixlo.GetMatrixArray());
89 if (ixup.GetNrows() > 0) fXupIndex.Use(ixup.GetNrows(),ixup.GetMatrixArray());
90 if (iclo.GetNrows() > 0) fCloIndex.Use(iclo.GetNrows(),iclo.GetMatrixArray());
91 if (icup.GetNrows() > 0) fCupIndex.Use(icup.GetNrows(),icup.GetMatrixArray());
92 fNxlo = ixlo.NonZeros();
93 fNxup = ixup.NonZeros();
94 fMclo = iclo.NonZeros();
95 fMcup = icup.NonZeros();
96
100
102 if (fMclo > 0) {
105 }
106 if (fMcup > 0) {
109 }
110 if (fNxlo > 0) {
113 }
114 if (fNxup > 0) {
117 }
118
119 fResidualNorm = 0.0;
120 fDualityGap = 0.0;
121}
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Copy constructor
126
128{
129 *this = another;
130}
131
132
133////////////////////////////////////////////////////////////////////////////////
134/// Calculate residuals, their norms, and duality complementarity gap,
135/// given a problem and variable set.
136
138{
139 TQpDataDens *prob = (TQpDataDens *) prob_in;
140
141 fRQ.ResizeTo(prob->fG); fRQ = prob->fG;
142 prob->Qmult(1.0,fRQ,1.0,vars->fX);
143
144 // calculate x^T (g+Qx) - contribution to the duality gap
145 Double_t gap = fRQ*vars->fX;
146
147 prob->ATransmult(1.0,fRQ,-1.0,vars->fY);
148 prob->CTransmult(1.0,fRQ,-1.0,vars->fZ);
149 if (fNxlo > 0) Add(fRQ,-1.0,vars->fGamma);
150 if (fNxup > 0) Add(fRQ, 1.0,vars->fPhi);
151
152 Double_t norm = 0.0;
153 Double_t componentNorm = fRQ.NormInf();
154 if (componentNorm > norm) norm = componentNorm;
155
156 fRA.ResizeTo(prob->fBa); fRA = prob->fBa;
157 prob->Amult(-1.0,fRA,1.0,vars->fX);
158
159 // contribution -d^T y to duality gap
160 gap -= prob->fBa*vars->fY;
161
162 componentNorm = fRA.NormInf();
163 if( componentNorm > norm ) norm = componentNorm;
164
165 fRC.ResizeTo(vars->fS); fRC = vars->fS;
166 prob->Cmult(-1.0,fRC,1.0,vars->fX);
167
168 componentNorm = fRC.NormInf();
169 if( componentNorm > norm ) norm = componentNorm;
170
171 fRz.ResizeTo(vars->fZ); fRz = vars->fZ;
172
173 if (fMclo > 0) {
174 Add(fRz,-1.0,vars->fLambda);
175
176 fRt.ResizeTo(vars->fS); fRt = vars->fS;
177 Add(fRt,-1.0,prob->GetSlowerBound());
179 Add(fRt,-1.0,vars->fT);
180
181 gap -= prob->fCloBound*vars->fLambda;
182
183 componentNorm = fRt.NormInf();
184 if( componentNorm > norm ) norm = componentNorm;
185 }
186
187 if (fMcup > 0) {
188 Add(fRz,1.0,vars->fPi);
189
190 fRu.ResizeTo(vars->fS); fRu = vars->fS;
191 Add(fRu,-1.0,prob->GetSupperBound() );
193 Add(fRu,1.0,vars->fU);
194
195 gap += prob->fCupBound*vars->fPi;
196
197 componentNorm = fRu.NormInf();
198 if( componentNorm > norm ) norm = componentNorm;
199 }
200
201 componentNorm = fRz.NormInf();
202 if( componentNorm > norm ) norm = componentNorm;
203
204 if (fNxlo > 0) {
205 fRv.ResizeTo(vars->fX); fRv = vars->fX;
206 Add(fRv,-1.0,prob->GetXlowerBound());
208 Add(fRv,-1.0,vars->fV);
209
210 gap -= prob->fXloBound*vars->fGamma;
211
212 componentNorm = fRv.NormInf();
213 if( componentNorm > norm ) norm = componentNorm;
214 }
215
216 if (fNxup > 0) {
217 fRw.ResizeTo(vars->fX); fRw = vars->fX;
218 Add(fRw,-1.0,prob->GetXupperBound());
220 Add(fRw,1.0,vars->fW);
221
222 gap += prob->fXupBound*vars->fPhi;
223
224 componentNorm = fRw.NormInf();
225 if (componentNorm > norm) norm = componentNorm;
226 }
227
228 fDualityGap = gap;
229 fResidualNorm = norm;
230}
231
232
233////////////////////////////////////////////////////////////////////////////////
234/// Modify the "complementarity" component of the residuals, by adding the pairwise
235/// products of the complementary variables plus a constant alpha to this term.
236
238{
239 if (fMclo > 0) AddElemMult(fRlambda,1.0,vars->fT,vars->fLambda);
240 if (fMcup > 0) AddElemMult(fRpi ,1.0,vars->fU,vars->fPi);
241 if (fNxlo > 0) AddElemMult(fRgamma ,1.0,vars->fV,vars->fGamma);
242 if (fNxup > 0) AddElemMult(fRphi ,1.0,vars->fW,vars->fPhi);
243
244 if (alpha != 0.0) {
245 if (fMclo > 0) fRlambda.AddSomeConstant(alpha,fCloIndex);
246 if (fMcup > 0) fRpi .AddSomeConstant(alpha,fCupIndex);
247 if (fNxlo > 0) fRgamma .AddSomeConstant(alpha,fXloIndex);
248 if (fNxup > 0) fRphi .AddSomeConstant(alpha,fXupIndex);
249 }
250}
251
252
253////////////////////////////////////////////////////////////////////////////////
254/// Set the "complementarity" component of the residuals to the pairwise products of
255/// the complementary variables plus a constant alpha .
256
258{
259 this->Clear_r3();
260 this->Add_r3_xz_alpha(vars,alpha);
261}
262
263
264////////////////////////////////////////////////////////////////////////////////
265/// set the complementarity component of the residuals to 0.
266
268{
269 if (fMclo > 0) fRlambda.Zero();
270 if (fMcup > 0) fRpi .Zero();
271 if (fNxlo > 0) fRgamma .Zero();
272 if (fNxup > 0) fRphi .Zero();
273}
274
275
276////////////////////////////////////////////////////////////////////////////////
277/// set the noncomplementarity components of the residual (the terms arising from
278/// the linear equalities in the KKT conditions) to 0.
279
281{
282 fRQ.Zero();
283 fRA.Zero();
284 fRC.Zero();
285 fRz.Zero();
286 if (fNxlo > 0) fRv.Zero();
287 if (fNxup > 0) fRw.Zero();
288 if (fMclo > 0) fRt.Zero();
289 if (fMcup > 0) fRu.Zero();
290}
291
292
293////////////////////////////////////////////////////////////////////////////////
294/// Perform the projection operation required by Gondzio algorithm: replace each
295/// component r3_i of the complementarity component of the residuals by r3p_i-r3_i,
296/// where r3p_i is the projection of r3_i onto the box [rmin, rmax]. Then if the
297/// resulting value is less than -rmax, replace it by -rmax.
298
300{
301 if (fMclo > 0) {
302 GondzioProjection(fRlambda,rmin,rmax);
304 }
305 if (fMcup > 0) {
306 GondzioProjection(fRpi,rmin,rmax);
308 }
309 if (fNxlo > 0) {
310 GondzioProjection(fRgamma,rmin,rmax);
312 }
313 if (fNxup > 0) {
314 GondzioProjection(fRphi,rmin,rmax);
316 }
317}
318
319
320////////////////////////////////////////////////////////////////////////////////
321/// Check if vector elements as selected through array indices are non-zero
322
324{
325 if (fNxlo > 0 &&
326 (!fRv .MatchesNonZeroPattern(fXloIndex) ||
328 return kFALSE;
329 }
330
331 if (fNxup > 0 &&
332 (!fRw .MatchesNonZeroPattern(fXupIndex) ||
334 return kFALSE;
335 }
336 if (fMclo > 0 &&
337 (!fRt .MatchesNonZeroPattern(fCloIndex) ||
339 return kFALSE;
340 }
341
342 if (fMcup > 0 &&
343 (!fRu .MatchesNonZeroPattern(fCupIndex) ||
345 return kFALSE;
346 }
347
348 return kTRUE;
349}
350
351
352////////////////////////////////////////////////////////////////////////////////
353/// Replace each component r3_i of the complementarity component of the residuals
354/// by r3p_i-r3_i, where r3p_i is the projection of r3_i onto the box [rmin, rmax].
355/// Then if the resulting value is less than -rmax, replace it by -rmax.
356
358{
359 Double_t * ep = v.GetMatrixArray();
360 const Double_t * const fep = ep+v.GetNrows();
361
362 while (ep < fep) {
363 if (*ep < rmin)
364 *ep = rmin - *ep;
365 else if (*ep > rmax)
366 *ep = rmax - *ep;
367 else
368 *ep = 0.0;
369 if (*ep < -rmax) *ep = -rmax;
370 ep++;
371 }
372}
373
374
375////////////////////////////////////////////////////////////////////////////////
376/// Assignment operator
377
379{
380 if (this != &source) {
381 TObject::operator=(source);
382
383 fNx = source.fNx;
384 fMy = source.fMy;
385 fMz = source.fMz;
386
387 fNxup = source.fNxup;
388 fNxlo = source.fNxlo;
389 fMcup = source.fMcup;
390 fMclo = source.fMclo;
391
396
397 fRQ .ResizeTo(source.fRQ); fRQ = source.fRQ;
398 fRA .ResizeTo(source.fRA); fRA = source.fRA;
399 fRC .ResizeTo(source.fRC); fRC = source.fRC;
400 fRz .ResizeTo(source.fRz); fRz = source.fRz;
401 fRv .ResizeTo(source.fRv); fRv = source.fRv;
402 fRw .ResizeTo(source.fRw); fRw = source.fRw;
403 fRt .ResizeTo(source.fRt); fRt = source.fRt;
404 fRu .ResizeTo(source.fRu); fRu = source.fRu;
405 fRgamma .ResizeTo(source.fRgamma); fRgamma = source.fRgamma;
406 fRphi .ResizeTo(source.fRphi); fRphi = source.fRphi;
407 fRlambda.ResizeTo(source.fRlambda); fRlambda = source.fRlambda;
408 fRpi .ResizeTo(source.fRpi); fRpi = source.fRpi;
409
410 // LM: copy also these data members
412 fDualityGap = source.fDualityGap;
413 }
414 return *this;
415}
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
#define ClassImp(name)
Definition Rtypes.h:377
Mother of all ROOT objects.
Definition TObject.h:41
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition TObject.h:296
Data for the general QP formulation.
Definition TQpDataBase.h:61
TVectorD fXloBound
Definition TQpDataBase.h:81
virtual TVectorD & GetSlowerBound()
TVectorD fG
Definition TQpDataBase.h:77
TVectorD fXupBound
Definition TQpDataBase.h:79
virtual TVectorD & GetXlowerBound()
TVectorD fCupBound
Definition TQpDataBase.h:83
virtual TVectorD & GetSupperBound()
TVectorD fCloBound
Definition TQpDataBase.h:85
TVectorD fBa
Definition TQpDataBase.h:78
virtual TVectorD & GetXupperBound()
Data for the dense QP formulation.
Definition TQpDataDens.h:63
void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fQ*x)
void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fC*x)
void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fA*x)
void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fA^T*x)
void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fC^T*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.
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:453
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition TVectorT.cxx:603
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:294
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:349
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition TVectorT.cxx:544
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition TVectorT.cxx:620
Element * GetMatrixArray()
Definition TVectorT.h:78
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.