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