Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
60 ClassImp(TQpResidual);
61 
62 ////////////////////////////////////////////////////////////////////////////////
63 /// Constructor
64 
65 TQpResidual::TQpResidual()
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);
106  fRlambda.ResizeTo(fMz);
107  }
108  if (fMcup > 0) {
109  fRu.ResizeTo(fMz);
110  fRpi.ResizeTo(fMz);
111  }
112  if (fNxlo > 0) {
113  fRv.ResizeTo(fNx);
114  fRgamma.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 
139 void TQpResidual::CalcResids(TQpDataBase *prob_in,TQpVar *vars)
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());
180  fRt.SelectNonZeros(fCloIndex);
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() );
194  fRu.SelectNonZeros(fCupIndex);
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());
209  fRv.SelectNonZeros(fXloIndex);
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());
221  fRw.SelectNonZeros(fXupIndex);
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 
269 void TQpResidual::Clear_r3()
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 
282 void TQpResidual::Clear_r1r2()
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);
305  fRlambda.SelectNonZeros(fCloIndex);
306  }
307  if (fMcup > 0) {
308  GondzioProjection(fRpi,rmin,rmax);
309  fRpi.SelectNonZeros(fCupIndex);
310  }
311  if (fNxlo > 0) {
312  GondzioProjection(fRgamma,rmin,rmax);
313  fRgamma.SelectNonZeros(fXloIndex);
314  }
315  if (fNxup > 0) {
316  GondzioProjection(fRphi,rmin,rmax);
317  fRphi.SelectNonZeros(fXupIndex);
318  }
319 }
320 
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 /// Check if vector elements as selected through array indices are non-zero
324 
325 Bool_t TQpResidual::ValidNonZeroPattern()
326 {
327  if (fNxlo > 0 &&
328  (!fRv .MatchesNonZeroPattern(fXloIndex) ||
329  !fRgamma.MatchesNonZeroPattern(fXloIndex)) ) {
330  return kFALSE;
331  }
332 
333  if (fNxup > 0 &&
334  (!fRw .MatchesNonZeroPattern(fXupIndex) ||
335  !fRphi.MatchesNonZeroPattern(fXupIndex)) ) {
336  return kFALSE;
337  }
338  if (fMclo > 0 &&
339  (!fRt .MatchesNonZeroPattern(fCloIndex) ||
340  !fRlambda.MatchesNonZeroPattern(fCloIndex)) ) {
341  return kFALSE;
342  }
343 
344  if (fMcup > 0 &&
345  (!fRu .MatchesNonZeroPattern(fCupIndex) ||
346  !fRpi.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 
380 TQpResidual &TQpResidual::operator=(const TQpResidual &source)
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 }