Logo ROOT   6.13/01
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 // 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 
60 ClassImp(TQpResidual);
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 
83 TQpResidual::TQpResidual(Int_t nx,Int_t my,Int_t mz,TVectorD &ixlo,TVectorD &ixup,
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 
99  fRQ.ResizeTo(fNx);
100  fRA.ResizeTo(fMy);
101  fRC.ResizeTo(fMz);
102 
103  fRz.ResizeTo(fMz);
104  if (fMclo > 0) {
105  fRt.ResizeTo(fMz);
107  }
108  if (fMcup > 0) {
109  fRu.ResizeTo(fMz);
110  fRpi.ResizeTo(fMz);
111  }
112  if (fNxlo > 0) {
113  fRv.ResizeTo(fNx);
115  }
116  if (fNxup > 0) {
117  fRw.ResizeTo(fNx);
118  fRphi.ResizeTo(fNx);
119  }
120 
121  fResidualNorm = 0.0;
122  fDualityGap = 0.0;
123 }
124 
125 
126 ////////////////////////////////////////////////////////////////////////////////
127 /// Copy constructor
128 
129 TQpResidual::TQpResidual(const TQpResidual &another) : TObject(another)
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 
239 void TQpResidual::Add_r3_xz_alpha(TQpVar *vars,Double_t alpha)
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 
259 void TQpResidual::Set_r3_xz_alpha(TQpVar *vars,Double_t alpha)
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 
301 void TQpResidual::Project_r3(Double_t rmin,Double_t rmax)
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 
359 void TQpResidual::GondzioProjection(TVectorD &v,Double_t rmin,Double_t rmax)
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 
394  fXupIndex.ResizeTo(source.fXupIndex); fXupIndex = source.fXupIndex;
395  fXloIndex.ResizeTo(source.fXloIndex); fXloIndex = source.fXloIndex;
396  fCupIndex.ResizeTo(source.fCupIndex); fCupIndex = source.fCupIndex;
397  fCloIndex.ResizeTo(source.fCloIndex); fCupIndex = source.fCupIndex;
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
413  fResidualNorm = source.fResidualNorm;
414  fDualityGap = source.fDualityGap;
415  }
416  return *this;
417 }
TVectorD fRt
Definition: TQpResidual.h:94
TVectorD fXupBound
Definition: TQpDataBase.h:79
virtual void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fQ*x)
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
void CalcResids(TQpDataBase *problem, TQpVar *vars)
Calculate residuals, their norms, and duality complementarity gap, given a problem and variable set...
TVectorD fG
Definition: TQpDataBase.h:77
TVectorD fRphi
Definition: TQpResidual.h:97
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixT.cxx:2933
TVectorD fRgamma
Definition: TQpResidual.h:96
Bool_t ValidNonZeroPattern()
Check if vector elements as selected through array indices are non-zero.
virtual TVectorD & GetXlowerBound()
Definition: TQpDataBase.h:110
TVectorD fRlambda
Definition: TQpResidual.h:98
TVectorD fRpi
Definition: TQpResidual.h:99
Int_t GetNrows() const
Definition: TVectorT.h:75
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 fXupIndex
Definition: TQpResidual.h:80
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, where r3p_i is the projection of r3_i onto the box [rmin, rmax].
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
TVectorD fX
Definition: TQpVar.h:91
Double_t fMclo
Definition: TQpResidual.h:77
virtual void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC*x)
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:347
TVectorD fXloIndex
Definition: TQpResidual.h:81
virtual void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA*x)
TVectorD fCloBound
Definition: TQpDataBase.h:85
Double_t fDualityGap
Definition: TQpResidual.h:66
TQpResidual & operator=(const TQpResidual &source)
Assignment operator.
TVectorD fRA
Definition: TQpResidual.h:89
TVectorD fW
Definition: TQpVar.h:99
virtual void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC^T*x)
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...
Element * GetMatrixArray()
Definition: TVectorT.h:78
TVectorD fCloIndex
Definition: TQpResidual.h:83
TVectorD fRz
Definition: TQpResidual.h:91
TVectorD fRv
Definition: TQpResidual.h:92
SVector< double, 2 > v
Definition: Dict.h:5
void Clear_r3()
set the complementarity component of the residuals to 0.
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:451
TVectorD fV
Definition: TQpVar.h:96
TVectorD fXloBound
Definition: TQpDataBase.h:81
Double_t fMcup
Definition: TQpResidual.h:76
virtual void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA^T*x)
TVectorD fRu
Definition: TQpResidual.h:95
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:618
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:542
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...
TVectorD fBa
Definition: TQpDataBase.h:78
Double_t fNxup
Definition: TQpResidual.h:74
TVectorD fY
Definition: TQpVar.h:93
TVectorD fLambda
Definition: TQpVar.h:103
TVectorD fU
Definition: TQpVar.h:105
Double_t fNxlo
Definition: TQpResidual.h:75
TVectorD fPhi
Definition: TQpVar.h:97
TVectorD fGamma
Definition: TQpVar.h:100
void Clear_r1r2()
set the noncomplementarity components of the residual (the terms arising from the linear equalities i...
TQpResidual()
Constructor.
Definition: TQpResidual.cxx:65
virtual TVectorD & GetXupperBound()
Definition: TQpDataBase.h:108
Definition: TQpVar.h:59
TVectorD fCupBound
Definition: TQpDataBase.h:83
TVectorD fS
Definition: TQpVar.h:92
TVectorD fPi
Definition: TQpVar.h:106
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:1842
TVectorD fT
Definition: TQpVar.h:102
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1236
TVectorD fCupIndex
Definition: TQpResidual.h:82
TVectorD fRQ
Definition: TQpResidual.h:88
void AddSomeConstant(Element val, const TVectorT< Element > &select)
Add to vector elements as selected through array select the value val.
Definition: TVectorT.cxx:1282
virtual TVectorD & GetSlowerBound()
Definition: TQpDataBase.h:114
Double_t fResidualNorm
Definition: TQpResidual.h:65
TVectorD fZ
Definition: TQpVar.h:94
virtual TVectorD & GetSupperBound()
Definition: TQpDataBase.h:112
TVectorD fRC
Definition: TQpResidual.h:90
TVectorD fRw
Definition: TQpResidual.h:93