Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TQpVar.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 /** \class TQpVar
44 
45 Class containing the variables for the general QP formulation
46 
47 In terms of in our abstract problem formulation, these variables are
48 the vectors x, y, z and s.
49 */
50 
51 #include "Riostream.h"
52 #include "TQpVar.h"
53 #include "TMatrixD.h"
54 
55 ClassImp(TQpVar);
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// Default constructor
59 
60 TQpVar::TQpVar()
61 {
62  fNx = 0;
63  fMy = 0;
64  fMz = 0;
65  fNxup = 0;
66  fNxlo = 0;
67  fMcup = 0;
68  fMclo = 0;
69  fNComplementaryVariables = 0;
70 }
71 
72 
73 ////////////////////////////////////////////////////////////////////////////////
74 /// Constructor
75 
76 TQpVar::TQpVar(TVectorD &x_in,TVectorD &s_in,TVectorD &y_in,TVectorD &z_in,
77  TVectorD &v_in,TVectorD &gamma_in,TVectorD &w_in,TVectorD &phi_in,
78  TVectorD &t_in,TVectorD &lambda_in,TVectorD &u_in,TVectorD &pi_in,
79  TVectorD &ixlow_in,TVectorD &ixupp_in,TVectorD &iclow_in,TVectorD &icupp_in)
80 {
81  if (x_in .GetNrows() > 0) fX. Use(x_in .GetNrows(),x_in .GetMatrixArray());
82  if (s_in .GetNrows() > 0) fS. Use(s_in .GetNrows(),s_in .GetMatrixArray());
83  if (y_in .GetNrows() > 0) fY. Use(y_in .GetNrows(),y_in .GetMatrixArray());
84  if (z_in .GetNrows() > 0) fZ. Use(z_in .GetNrows(),z_in .GetMatrixArray());
85  if (v_in .GetNrows() > 0) fV. Use(v_in .GetNrows(),v_in .GetMatrixArray());
86  if (phi_in .GetNrows() > 0) fPhi. Use(phi_in .GetNrows(),phi_in .GetMatrixArray());
87  if (w_in .GetNrows() > 0) fW. Use(w_in .GetNrows(),w_in .GetMatrixArray());
88  if (gamma_in .GetNrows() > 0) fGamma. Use(gamma_in .GetNrows(),gamma_in .GetMatrixArray());
89  if (t_in .GetNrows() > 0) fT. Use(t_in .GetNrows(),t_in .GetMatrixArray());
90  if (lambda_in.GetNrows() > 0) fLambda. Use(lambda_in.GetNrows(),lambda_in.GetMatrixArray());
91  if (u_in .GetNrows() > 0) fU. Use(u_in .GetNrows(),u_in .GetMatrixArray());
92  if (pi_in .GetNrows() > 0) fPi. Use(pi_in .GetNrows(),pi_in .GetMatrixArray());
93  if (ixlow_in .GetNrows() > 0) fXloIndex.Use(ixlow_in .GetNrows(),ixlow_in .GetMatrixArray());
94  if (ixupp_in .GetNrows() > 0) fXupIndex.Use(ixupp_in .GetNrows(),ixupp_in .GetMatrixArray());
95  if (iclow_in .GetNrows() > 0) fCloIndex.Use(iclow_in .GetNrows(),iclow_in .GetMatrixArray());
96  if (icupp_in .GetNrows() > 0) fCupIndex.Use(icupp_in .GetNrows(),icupp_in .GetMatrixArray());
97 
98  fNx = fX.GetNrows();
99  fMy = fY.GetNrows();
100  fMz = fZ.GetNrows();
101 
102  R__ASSERT(fNx == fXloIndex.GetNrows() || 0 == fXloIndex.GetNrows());
103  R__ASSERT(fNx == fXloIndex.GetNrows() || 0 == fXloIndex.GetNrows());
104  R__ASSERT(fMz == fCloIndex.GetNrows() || 0 == fCloIndex.GetNrows());
105  R__ASSERT(fMz == fCupIndex.GetNrows() || 0 == fCupIndex.GetNrows());
106 
107  fNxlo = fXloIndex.NonZeros();
108  fNxup = fXupIndex.NonZeros();
109  fMclo = fCloIndex.NonZeros();
110  fMcup = fCupIndex.NonZeros();
111  fNComplementaryVariables = fMclo+fMcup+fNxlo+fNxup;
112 
113  R__ASSERT(fMz == fS.GetNrows());
114  R__ASSERT(fNx == fV .GetNrows() || (0 == fV .GetNrows() && fNxlo == 0));
115  R__ASSERT(fNx == fGamma .GetNrows() || (0 == fGamma .GetNrows() && fNxlo == 0));
116 
117  R__ASSERT(fNx == fW .GetNrows() || (0 == fW .GetNrows() && fNxup == 0));
118  R__ASSERT(fNx == fPhi .GetNrows() || (0 == fPhi .GetNrows() && fNxup == 0));
119 
120  R__ASSERT(fMz == fT .GetNrows() || (0 == fT .GetNrows() && fMclo == 0));
121  R__ASSERT(fMz == fLambda.GetNrows() || (0 == fLambda.GetNrows() && fMclo == 0));
122 
123  R__ASSERT(fMz == fU .GetNrows() || (0 == fU .GetNrows() && fMcup == 0));
124  R__ASSERT(fMz == fPi .GetNrows() || (0 == fPi .GetNrows() && fMcup == 0));
125 }
126 
127 
128 ////////////////////////////////////////////////////////////////////////////////
129 /// Constructor
130 
131 TQpVar::TQpVar(Int_t nx,Int_t my,Int_t mz,TVectorD &ixlow,TVectorD &ixupp,
132  TVectorD &iclow,TVectorD &icupp)
133 {
134  R__ASSERT(nx == ixlow.GetNrows() || 0 == ixlow.GetNrows());
135  R__ASSERT(nx == ixlow.GetNrows() || 0 == ixlow.GetNrows());
136  R__ASSERT(mz == iclow.GetNrows() || 0 == iclow.GetNrows());
137  R__ASSERT(mz == icupp.GetNrows() || 0 == icupp.GetNrows());
138 
139  fNxlo = ixlow.NonZeros();
140  fNxup = ixupp.NonZeros();
141  fMclo = iclow.NonZeros();
142  fMcup = icupp.NonZeros();
143 
144  if (ixlow.GetNrows() > 0) fXloIndex.Use(ixlow.GetNrows(),ixlow.GetMatrixArray());
145  if (ixupp.GetNrows() > 0) fXupIndex.Use(ixupp.GetNrows(),ixupp.GetMatrixArray());
146  if (iclow.GetNrows() > 0) fCloIndex.Use(iclow.GetNrows(),iclow.GetMatrixArray());
147  if (icupp.GetNrows() > 0) fCupIndex.Use(icupp.GetNrows(),icupp.GetMatrixArray());
148 
149  fNx = nx;
150  fMy = my;
151  fMz = mz;
152 
153  if (fMclo > 0) {
154  fT.ResizeTo(fMz);
155  fLambda.ResizeTo(fMz);
156  }
157  if (fMcup > 0) {
158  fU.ResizeTo(fMz);
159  fPi.ResizeTo(fMz);
160  }
161  if (fNxlo > 0) {
162  fV.ResizeTo(fNx);
163  fGamma.ResizeTo(fNx);
164  }
165 
166  if (fNxup > 0) {
167  fW.ResizeTo(fNx);
168  fPhi.ResizeTo(fNx);
169  }
170 
171  fS.ResizeTo(fMz);
172  fX.ResizeTo(fNx);
173  fY.ResizeTo(fMy);
174  fZ.ResizeTo(fMz);
175  fNComplementaryVariables = fMclo+fMcup+fNxlo+fNxup;
176 }
177 
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 /// Copy constructor
181 
182 TQpVar::TQpVar(const TQpVar &another) : TObject(another)
183 {
184  *this = another;
185 }
186 
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// compute complementarity gap, obtained by taking the inner product of the
190 /// complementary vectors and dividing by the total number of components
191 /// computes mu = (t'lambda +u'pi + v'gamma + w'phi)/(mclow+mcupp+nxlow+nxupp)
192 
193 Double_t TQpVar::GetMu()
194 {
195  Double_t mu = 0.0;
196  if (fNComplementaryVariables > 0 ) {
197  if (fMclo > 0) mu += fT*fLambda;
198  if (fMcup > 0) mu += fU*fPi;
199  if (fNxlo > 0) mu += fV*fGamma;
200  if (fNxup > 0) mu += fW*fPhi;
201 
202  mu /= fNComplementaryVariables;
203  }
204  return mu;
205 }
206 
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 /// Compute the complementarity gap resulting from a step of length "alpha" along
210 /// direction "step"
211 
212 Double_t TQpVar::MuStep(TQpVar *step,Double_t alpha)
213 {
214  Double_t mu = 0.0;
215  if (fNComplementaryVariables > 0) {
216  if (fMclo > 0)
217  mu += (fT+alpha*step->fT)*(fLambda+alpha*step->fLambda);
218  if (fMcup > 0)
219  mu += (fU+alpha*step->fU)*(fPi+alpha*step->fPi);
220  if (fNxlo > 0)
221  mu += (fV+alpha*step->fV)*(fGamma+alpha*step->fGamma);
222  if (fNxup > 0)
223  mu += (fW+alpha*step->fW)*(fPhi+alpha*step->fPhi);
224  mu /= fNComplementaryVariables;
225  }
226  return mu;
227 }
228 
229 
230 ////////////////////////////////////////////////////////////////////////////////
231 /// Perform a "saxpy" operation on all data vectors : x += alpha*y
232 
233 void TQpVar::Saxpy(TQpVar *b,Double_t alpha)
234 {
235  Add(fX,alpha,b->fX);
236  Add(fY,alpha,b->fY);
237  Add(fZ,alpha,b->fZ);
238  Add(fS,alpha,b->fS);
239  if (fMclo > 0) {
240  R__ASSERT((b->fT) .MatchesNonZeroPattern(fCloIndex) &&
241  (b->fLambda).MatchesNonZeroPattern(fCloIndex));
242 
243  Add(fT ,alpha,b->fT);
244  Add(fLambda,alpha,b->fLambda);
245  }
246  if (fMcup > 0) {
247  R__ASSERT((b->fU) .MatchesNonZeroPattern(fCupIndex) &&
248  (b->fPi).MatchesNonZeroPattern(fCupIndex));
249 
250  Add(fU ,alpha,b->fU);
251  Add(fPi,alpha,b->fPi);
252  }
253  if (fNxlo > 0) {
254  R__ASSERT((b->fV) .MatchesNonZeroPattern(fXloIndex) &&
255  (b->fGamma).MatchesNonZeroPattern(fXloIndex));
256 
257  Add(fV ,alpha,b->fV);
258  Add(fGamma,alpha,b->fGamma);
259  }
260  if (fNxup > 0) {
261  R__ASSERT((b->fW) .MatchesNonZeroPattern(fXupIndex) &&
262  (b->fPhi).MatchesNonZeroPattern(fXupIndex));
263 
264  Add(fW ,alpha,b->fW);
265  Add(fPhi,alpha,b->fPhi);
266  }
267 }
268 
269 
270 ////////////////////////////////////////////////////////////////////////////////
271 /// Perform a "negate" operation on all data vectors : x = -x
272 
273 void TQpVar::Negate()
274 {
275  fS *= -1.;
276  fX *= -1.;
277  fY *= -1.;
278  fZ *= -1.;
279  if (fMclo > 0) {
280  fT *= -1.;
281  fLambda *= -1.;
282  }
283  if (fMcup > 0) {
284  fU *= -1.;
285  fPi *= -1.;
286  }
287  if (fNxlo > 0) {
288  fV *= -1.;
289  fGamma *= -1.;
290  }
291  if (fNxup > 0) {
292  fW *= -1.;
293  fPhi *= -1.;
294  }
295 }
296 
297 
298 ////////////////////////////////////////////////////////////////////////////////
299 /// calculate the largest alpha in (0,1] such that the/ nonnegative variables stay
300 /// nonnegative in the given search direction. In the general QP problem formulation
301 /// this is the largest value of alpha such that
302 /// (t,u,v,w,lambda,pi,phi,gamma) + alpha * (b->t,b->u,b->v,b->w,b->lambda,b->pi,
303 /// b->phi,b->gamma) >= 0.
304 
305 Double_t TQpVar::StepBound(TQpVar *b)
306 {
307  Double_t maxStep = 1.0;
308 
309  if (fMclo > 0 ) {
310  R__ASSERT(fT .SomePositive(fCloIndex));
311  R__ASSERT(fLambda.SomePositive(fCloIndex));
312 
313  maxStep = this->StepBound(fT, b->fT, maxStep);
314  maxStep = this->StepBound(fLambda,b->fLambda,maxStep);
315  }
316 
317  if (fMcup > 0 ) {
318  R__ASSERT(fU .SomePositive(fCupIndex));
319  R__ASSERT(fPi.SomePositive(fCupIndex));
320 
321  maxStep = this->StepBound(fU, b->fU, maxStep);
322  maxStep = this->StepBound(fPi,b->fPi,maxStep);
323  }
324 
325  if (fNxlo > 0 ) {
326  R__ASSERT(fV .SomePositive(fXloIndex));
327  R__ASSERT(fGamma.SomePositive(fXloIndex));
328 
329  maxStep = this->StepBound(fV, b->fV, maxStep);
330  maxStep = this->StepBound(fGamma,b->fGamma,maxStep);
331  }
332 
333  if (fNxup > 0 ) {
334  R__ASSERT(fW .SomePositive(fXupIndex));
335  R__ASSERT(fPhi.SomePositive(fXupIndex));
336 
337  maxStep = this->StepBound(fW, b->fW, maxStep);
338  maxStep = this->StepBound(fPhi,b->fPhi,maxStep);
339  }
340 
341  return maxStep;
342 }
343 
344 
345 ////////////////////////////////////////////////////////////////////////////////
346 /// Find the maximum stepsize of v in direction dir
347 /// before violating the nonnegativity constraints
348 
349 Double_t TQpVar::StepBound(TVectorD &v,TVectorD &dir,Double_t maxStep)
350 {
351  if (!AreCompatible(v,dir)) {
352  ::Error("StepBound(TVectorD &,TVectorD &,Double_t)","vector's not compatible");
353  return kFALSE;
354  }
355 
356  const Int_t n = v.GetNrows();
357  const Double_t * const pD = dir.GetMatrixArray();
358  const Double_t * const pV = v.GetMatrixArray();
359 
360  Double_t bound = maxStep;
361  for (Int_t i = 0; i < n; i++) {
362  Double_t tmp = pD[i];
363  if ( pV[i] >= 0 && tmp < 0 ) {
364  tmp = -pV[i]/tmp;
365  if (tmp < bound)
366  bound = tmp;
367  }
368  }
369  return bound;
370 }
371 
372 
373 ////////////////////////////////////////////////////////////////////////////////
374 /// Is the current position an interior point ?
375 
376 Bool_t TQpVar::IsInteriorPoint()
377 {
378  Bool_t interior = kTRUE;
379  if (fMclo > 0)
380  interior = interior &&
381  fT.SomePositive(fCloIndex) && fLambda.SomePositive(fCloIndex);
382 
383  if (fMcup > 0)
384  interior = interior &&
385  fU.SomePositive(fCupIndex) && fPi.SomePositive(fCupIndex);
386 
387  if (fNxlo > 0)
388  interior = interior &&
389  fV.SomePositive(fXloIndex) && fGamma.SomePositive(fXloIndex);
390 
391  if (fNxup > 0)
392  interior = interior &&
393  fW.SomePositive(fXupIndex) && fPhi.SomePositive(fXupIndex);
394 
395  return interior;
396 }
397 
398 
399 ////////////////////////////////////////////////////////////////////////////////
400 /// Performs the same function as StepBound, and supplies additional information about
401 /// which component of the nonnegative variables is responsible for restricting alpha.
402 /// In terms of the abstract formulation, the components have the following meanings :
403 ///
404 /// primalValue : the value of the blocking component of the primal variables (u,t,v,w).
405 /// primalStep : the corresponding value of the blocking component of the primal step
406 /// variables (b->u,b->t,b->v,b->w)
407 /// dualValue : the value of the blocking component of the dual variables/
408 /// (lambda,pi,phi,gamma).
409 /// dualStep : the corresponding value of the blocking component of the dual step
410 /// variables (b->lambda,b->pi,b->phi,b->gamma)
411 /// firstOrSecond : 1 if the primal step is blocking,
412 /// 2 if the dual step is block,
413 /// 0 if no step is blocking.
414 
415 Double_t TQpVar::FindBlocking(TQpVar *step,
416  Double_t &primalValue,
417  Double_t &primalStep,
418  Double_t &dualValue,
419  Double_t &dualStep,
420  Int_t &fIrstOrSecond)
421 {
422  fIrstOrSecond = 0;
423  Double_t alpha = 1.0;
424  if (fMclo > 0)
425  alpha = FindBlocking(fT,step->fT,fLambda,step->fLambda,alpha,
426  primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
427 
428  if (fMcup > 0)
429  alpha = FindBlocking(fU,step->fU,fPi,step->fPi,alpha,
430  primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
431 
432  if (fNxlo > 0)
433  alpha = FindBlocking(fV,step->fV,fGamma,step->fGamma,alpha,
434  primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
435 
436  if (fNxup > 0)
437  alpha = FindBlocking(fW,step->fW,fPhi,step->fPhi,alpha,
438  primalValue,primalStep,dualValue,dualStep,fIrstOrSecond);
439 
440  return alpha;
441 }
442 
443 
444 ////////////////////////////////////////////////////////////////////////////////
445 /// See other FindBlocking function
446 
447 Double_t TQpVar::FindBlocking(TVectorD &w,TVectorD &wstep,TVectorD &u,TVectorD &ustep,
448  Double_t maxStep,Double_t &w_elt,Double_t &wstep_elt,Double_t &u_elt,
449  Double_t &ustep_elt,int& fIrst_or_second)
450 {
451  return FindBlockingSub(w.GetNrows(),
452  w.GetMatrixArray(), 1,
453  wstep.GetMatrixArray(),1,
454  u.GetMatrixArray(), 1,
455  ustep.GetMatrixArray(),1,
456  maxStep,
457  w_elt,wstep_elt,
458  u_elt,ustep_elt,
459  fIrst_or_second);
460 }
461 
462 
463 ////////////////////////////////////////////////////////////////////////////////
464 /// See FindBlocking function
465 
466 Double_t TQpVar::FindBlockingSub(Int_t n,
467  Double_t *w, Int_t incw,
468  Double_t *wstep,Int_t incwstep,
469  Double_t *u, Int_t incu,
470  Double_t *ustep,Int_t incustep,
471  Double_t maxStep,
472  Double_t &w_elt,Double_t &wstep_elt,
473  Double_t &u_elt,Double_t &ustep_elt,
474  Int_t &fIrst_or_second)
475 {
476  Double_t bound = maxStep;
477 
478  Int_t i = n-1;
479  Int_t lastBlocking = -1;
480 
481  // Search backward so that we fInd the blocking constraint of lowest
482  // index. We do this to make things consistent with MPI's MPI_MINLOC,
483  // which returns the processor with smallest rank where a min occurs.
484  //
485  // Still, going backward is ugly!
486  Double_t *pw = w +(n-1)*incw;
487  Double_t *pwstep = wstep+(n-1)*incwstep;
488  Double_t *pu = u +(n-1)*incu;
489  Double_t *pustep = ustep+(n-1)*incustep;
490 
491  while (i >= 0) {
492  Double_t tmp = *pwstep;
493  if (*pw > 0 && tmp < 0) {
494  tmp = -*pw/tmp;
495  if( tmp <= bound ) {
496  bound = tmp;
497  lastBlocking = i;
498  fIrst_or_second = 1;
499  }
500  }
501  tmp = *pustep;
502  if (*pu > 0 && tmp < 0) {
503  tmp = -*pu/tmp;
504  if( tmp <= bound ) {
505  bound = tmp;
506  lastBlocking = i;
507  fIrst_or_second = 2;
508  }
509  }
510 
511  i--;
512  if (i >= 0) {
513  // It is safe to decrement the pointers
514  pw -= incw;
515  pwstep -= incwstep;
516  pu -= incu;
517  pustep -= incustep;
518  }
519  }
520 
521  if (lastBlocking > -1) {
522  // fIll out the elements
523  w_elt = w[lastBlocking];
524  wstep_elt = wstep[lastBlocking];
525  u_elt = u[lastBlocking];
526  ustep_elt = ustep[lastBlocking];
527  }
528  return bound;
529 }
530 
531 
532 ////////////////////////////////////////////////////////////////////////////////
533 /// Sets components of (u,t,v,w) to alpha and of (lambda,pi,phi,gamma) to beta
534 
535 void TQpVar::InteriorPoint(Double_t alpha,Double_t beta)
536 {
537  fS.Zero();
538  fX.Zero();
539  fY.Zero();
540  fZ.Zero();
541 
542  if (fNxlo > 0) {
543  fV = alpha;
544  fV.SelectNonZeros(fXloIndex);
545  fGamma = beta;
546  fGamma.SelectNonZeros(fXloIndex);
547  }
548 
549  if (fNxup > 0) {
550  fW = alpha;
551  fW.SelectNonZeros(fXupIndex);
552  fPhi = beta;
553  fPhi.SelectNonZeros(fXupIndex);
554  }
555 
556  if (fMclo > 0 ) {
557  fT = alpha;
558  fT.SelectNonZeros(fCloIndex);
559  fLambda = beta;
560  fLambda.SelectNonZeros(fCloIndex);
561  }
562 
563  if (fMcup > 0) {
564  fU = alpha;
565  fU.SelectNonZeros(fCupIndex);
566  fPi = beta;
567  fPi.SelectNonZeros(fCupIndex);
568  }
569 }
570 
571 
572 ////////////////////////////////////////////////////////////////////////////////
573 /// The amount by which the current variables violate the non-negativity constraints.
574 
575 Double_t TQpVar::Violation()
576 {
577  Double_t viol = 0.0;
578  Double_t cmin;
579 
580  if (fNxlo > 0) {
581  cmin = fV.Min();
582  if (cmin < viol) viol = cmin;
583 
584  cmin = fGamma.Min();
585  if (cmin < viol) viol = cmin;
586  }
587  if (fNxup > 0) {
588  cmin = fW.Min();
589  if (cmin < viol) viol = cmin;
590 
591  cmin = fPhi.Min();
592  if (cmin < viol) viol = cmin;
593  }
594  if (fMclo > 0) {
595  cmin = fT.Min();
596  if (cmin < viol) viol = cmin;
597 
598  cmin = fLambda.Min();
599  if (cmin < viol) viol = cmin;
600  }
601  if (fMcup > 0) {
602  cmin = fU.Min();
603  if (cmin < viol) viol = cmin;
604 
605  cmin = fPi.Min();
606  if (cmin < viol) viol = cmin;
607  }
608 
609  return -viol;
610 }
611 
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 /// Add alpha to components of (u,t,v,w) and beta to components of (lambda,pi,phi,gamma)
615 
616 void TQpVar::ShiftBoundVariables(Double_t alpha,Double_t beta)
617 {
618  if (fNxlo > 0) {
619  fV .AddSomeConstant(alpha,fXloIndex);
620  fGamma.AddSomeConstant(beta, fXloIndex);
621  }
622  if (fNxup > 0) {
623  fW .AddSomeConstant(alpha,fXupIndex);
624  fPhi.AddSomeConstant(beta, fXupIndex);
625  }
626  if (fMclo > 0) {
627  fT .AddSomeConstant(alpha,fCloIndex);
628  fLambda.AddSomeConstant(beta, fCloIndex);
629  }
630  if (fMcup > 0) {
631  fU .AddSomeConstant(alpha,fCupIndex);
632  fPi.AddSomeConstant(beta, fCupIndex);
633  }
634 }
635 
636 
637 ////////////////////////////////////////////////////////////////////////////////
638 /// Print class members
639 
640 void TQpVar::Print(Option_t * /*option*/) const
641 {
642  std::cout << "fNx : " << fNx << std::endl;
643  std::cout << "fMy : " << fMy << std::endl;
644  std::cout << "fMz : " << fMz << std::endl;
645  std::cout << "fNxup: " << fNxup << std::endl;
646  std::cout << "fNxlo: " << fNxlo << std::endl;
647  std::cout << "fMcup: " << fMcup << std::endl;
648  std::cout << "fMclo: " << fMclo << std::endl;
649 
650  fXloIndex.Print("fXloIndex");
651  fXupIndex.Print("fXupIndex");
652  fCupIndex.Print("fCupIndex");
653  fCloIndex.Print("fCloIndex");
654 
655  fX.Print("fX");
656  fS.Print("fS");
657  fY.Print("fY");
658  fZ.Print("fZ");
659 
660  fV.Print("fV");
661  fPhi.Print("fPhi");
662 
663  fW.Print("fW");
664  fGamma.Print("fGamma");
665 
666  fT.Print("fT");
667  fLambda.Print("fLambda");
668 
669  fU.Print("fU");
670  fPi.Print("fPi");
671 }
672 
673 
674 ////////////////////////////////////////////////////////////////////////////////
675 /// Return the sum of the vector-norm1's
676 
677 Double_t TQpVar::Norm1()
678 {
679  Double_t norm = 0.0;
680  norm += fX.Norm1();
681  norm += fS.Norm1();
682  norm += fY.Norm1();
683  norm += fZ.Norm1();
684 
685  norm += fV.Norm1();
686  norm += fPhi.Norm1();
687  norm += fW.Norm1();
688  norm += fGamma.Norm1();
689  norm += fT.Norm1();
690  norm += fLambda.Norm1();
691  norm += fU.Norm1();
692  norm += fPi.Norm1();
693 
694  return norm;
695 }
696 
697 
698 ////////////////////////////////////////////////////////////////////////////////
699 /// Return the sum of the vector-normInf's
700 
701 Double_t TQpVar::NormInf()
702 {
703  Double_t norm = 0.0;
704 
705  Double_t tmp = fX.NormInf();
706  if (tmp > norm) norm = tmp;
707  tmp = fS.NormInf();
708  if (tmp > norm) norm = tmp;
709  tmp = fY.NormInf();
710  if (tmp > norm) norm = tmp;
711  tmp = fZ.NormInf();
712  if (tmp > norm) norm = tmp;
713 
714  tmp = fV.NormInf();
715  if (tmp > norm) norm = tmp;
716  tmp = fPhi.NormInf();
717  if (tmp > norm) norm = tmp;
718 
719  tmp = fW.NormInf();
720  if (tmp > norm) norm = tmp;
721  tmp = fGamma.NormInf();
722  if (tmp > norm) norm = tmp;
723 
724  tmp = fT.NormInf();
725  if (tmp > norm) norm = tmp;
726  tmp = fLambda.NormInf();
727  if (tmp > norm) norm = tmp;
728 
729  tmp = fU.NormInf();
730  if (tmp > norm) norm = tmp;
731  tmp = fPi.NormInf();
732  if (tmp > norm) norm = tmp;
733 
734  return norm;
735 }
736 
737 
738 ////////////////////////////////////////////////////////////////////////////////
739 /// Check that the variables conform to the non-zero indices
740 
741 Bool_t TQpVar::ValidNonZeroPattern()
742 {
743  if (fNxlo > 0 &&
744  ( !fV .MatchesNonZeroPattern(fXloIndex) ||
745  !fGamma.MatchesNonZeroPattern(fXloIndex) ) ) {
746  return kFALSE;
747  }
748 
749  if (fNxup > 0 &&
750  ( !fW .MatchesNonZeroPattern(fXupIndex) ||
751  !fPhi.MatchesNonZeroPattern(fXupIndex) ) ) {
752  return kFALSE;
753  }
754  if (fMclo > 0 &&
755  ( !fT .MatchesNonZeroPattern(fCloIndex) ||
756  !fLambda.MatchesNonZeroPattern(fCloIndex) ) ) {
757  return kFALSE;
758  }
759 
760  if (fMcup > 0 &&
761  ( !fU .MatchesNonZeroPattern(fCupIndex) ||
762  !fPi.MatchesNonZeroPattern(fCupIndex) ) ) {
763  return kFALSE;
764  }
765 
766  return kTRUE;
767 }
768 
769 
770 ////////////////////////////////////////////////////////////////////////////////
771 /// Assignment operator
772 
773 TQpVar &TQpVar::operator=(const TQpVar &source)
774 {
775  if (this != &source) {
776  TObject::operator=(source);
777  fNx = source.fNx;
778  fMy = source.fMy;
779  fMz = source.fMz;
780  fNxup = source.fNxup;
781  fNxlo = source.fNxlo;
782  fMcup = source.fMcup;
783  fMclo = source.fMclo;
784 
785  fXloIndex = source.fXloIndex;
786  fXupIndex = source.fXupIndex;
787  fCupIndex = source.fCupIndex;
788  fCloIndex = source.fCloIndex;
789 
790  fX .ResizeTo(source.fX); fX = source.fX;
791  fS .ResizeTo(source.fS); fS = source.fS;
792  fY .ResizeTo(source.fY); fY = source.fY;
793  fZ .ResizeTo(source.fZ); fZ = source.fZ;
794 
795  fV .ResizeTo(source.fV); fV = source.fV;
796  fPhi .ResizeTo(source.fPhi); fPhi = source.fPhi;
797 
798  fW .ResizeTo(source.fW); fW = source.fW;
799  fGamma .ResizeTo(source.fGamma) ; fGamma = source.fGamma;
800 
801  fT .ResizeTo(source.fT); fT = source.fT;
802  fLambda.ResizeTo(source.fLambda); fLambda = source.fLambda;
803 
804  fU .ResizeTo(source.fU); fU = source.fU;
805  fPi .ResizeTo(source.fPi); fPi = source.fPi;
806 
807  // LM: copy also this data member
808  fNComplementaryVariables = source.fNComplementaryVariables;
809  }
810  return *this;
811 }