Logo ROOT  
Reference Guide
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
45Class containing the variables for the general QP formulation
46
47In terms of in our abstract problem formulation, these variables are
48the vectors x, y, z and s.
49*/
50
51#include "Riostream.h"
52#include "TQpVar.h"
53#include "TMatrixD.h"
54
56
57////////////////////////////////////////////////////////////////////////////////
58/// Default constructor
59
61{
62 fNx = 0;
63 fMy = 0;
64 fMz = 0;
65 fNxup = 0;
66 fNxlo = 0;
67 fMcup = 0;
68 fMclo = 0;
70}
71
72
73////////////////////////////////////////////////////////////////////////////////
74/// Constructor
75
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
106
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
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);
156 }
157 if (fMcup > 0) {
158 fU.ResizeTo(fMz);
160 }
161 if (fNxlo > 0) {
162 fV.ResizeTo(fNx);
164 }
165
166 if (fNxup > 0) {
167 fW.ResizeTo(fNx);
169 }
170
171 fS.ResizeTo(fMz);
172 fX.ResizeTo(fNx);
173 fY.ResizeTo(fMy);
174 fZ.ResizeTo(fMz);
176}
177
178
179////////////////////////////////////////////////////////////////////////////////
180/// Copy constructor
181
182TQpVar::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
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
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
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);
225 }
226 return mu;
227}
228
229
230////////////////////////////////////////////////////////////////////////////////
231/// Perform a "saxpy" operation on all data vectors : x += alpha*y
232
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
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
306{
307 Double_t maxStep = 1.0;
308
309 if (fMclo > 0 ) {
310 R__ASSERT(fT .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));
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));
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));
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
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
377{
378 Bool_t interior = kTRUE;
379 if (fMclo > 0)
380 interior = interior &&
382
383 if (fMcup > 0)
384 interior = interior &&
386
387 if (fNxlo > 0)
388 interior = interior &&
390
391 if (fNxup > 0)
392 interior = interior &&
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
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
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
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
536{
537 fS.Zero();
538 fX.Zero();
539 fY.Zero();
540 fZ.Zero();
541
542 if (fNxlo > 0) {
543 fV = alpha;
545 fGamma = beta;
547 }
548
549 if (fNxup > 0) {
550 fW = alpha;
552 fPhi = beta;
554 }
555
556 if (fMclo > 0 ) {
557 fT = alpha;
559 fLambda = beta;
561 }
562
563 if (fMcup > 0) {
564 fU = alpha;
566 fPi = beta;
568 }
569}
570
571
572////////////////////////////////////////////////////////////////////////////////
573/// The amount by which the current variables violate the non-negativity constraints.
574
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
617{
618 if (fNxlo > 0) {
621 }
622 if (fNxup > 0) {
625 }
626 if (fMclo > 0) {
629 }
630 if (fMcup > 0) {
633 }
634}
635
636
637////////////////////////////////////////////////////////////////////////////////
638/// Print class members
639
640void 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
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
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
742{
743 if (fNxlo > 0 &&
744 ( !fV .MatchesNonZeroPattern(fXloIndex) ||
746 return kFALSE;
747 }
748
749 if (fNxup > 0 &&
750 ( !fW .MatchesNonZeroPattern(fXupIndex) ||
752 return kFALSE;
753 }
754 if (fMclo > 0 &&
755 ( !fT .MatchesNonZeroPattern(fCloIndex) ||
757 return kFALSE;
758 }
759
760 if (fMcup > 0 &&
761 ( !fU .MatchesNonZeroPattern(fCupIndex) ||
763 return kFALSE;
764 }
765
766 return kTRUE;
767}
768
769
770////////////////////////////////////////////////////////////////////////////////
771/// Assignment operator
772
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
809 }
810 return *this;
811}
#define b(i)
Definition: RSha256.hxx:100
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
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 .
Mother of all ROOT objects.
Definition: TObject.h:37
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:283
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
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
Int_t fNxlo
Definition: TQpVar.h:67
virtual Bool_t ValidNonZeroPattern()
Check that the variables conform to the non-zero indices.
Definition: TQpVar.cxx:741
virtual Double_t Violation()
The amount by which the current variables violate the non-negativity constraints.
Definition: TQpVar.cxx:575
TVectorD fLambda
Definition: TQpVar.h:103
Int_t fNComplementaryVariables
Definition: TQpVar.h:88
TQpVar & operator=(const TQpVar &source)
Assignment operator.
Definition: TQpVar.cxx:773
Int_t fMz
Definition: TQpVar.h:65
TQpVar()
Default constructor.
Definition: TQpVar.cxx:60
TVectorD fGamma
Definition: TQpVar.h:100
Int_t fMy
Definition: TQpVar.h:64
Int_t fNx
Definition: TQpVar.h:63
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:535
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:349
TVectorD fT
Definition: TQpVar.h:102
TVectorD fPi
Definition: TQpVar.h:106
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:212
virtual Double_t GetMu()
compute complementarity gap, obtained by taking the inner product of the complementary vectors and di...
Definition: TQpVar.cxx:193
TVectorD fS
Definition: TQpVar.h:92
Int_t fMcup
Definition: TQpVar.h:68
TVectorD fXloIndex
Definition: TQpVar.h:72
virtual Double_t Norm1()
Return the sum of the vector-norm1's.
Definition: TQpVar.cxx:677
TVectorD fPhi
Definition: TQpVar.h:97
virtual Double_t NormInf()
Return the sum of the vector-normInf's.
Definition: TQpVar.cxx:701
TVectorD fV
Definition: TQpVar.h:96
TVectorD fW
Definition: TQpVar.h:99
Int_t fNxup
Definition: TQpVar.h:66
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:616
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:466
virtual Bool_t IsInteriorPoint()
Is the current position an interior point ?
Definition: TQpVar.cxx:376
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:447
TVectorD fXupIndex
Definition: TQpVar.h:73
Int_t fMclo
Definition: TQpVar.h:69
virtual void Print(Option_t *option="") const
Print class members.
Definition: TQpVar.cxx:640
TVectorD fY
Definition: TQpVar.h:93
TVectorD fCloIndex
Definition: TQpVar.h:75
TVectorD fZ
Definition: TQpVar.h:94
virtual void Negate()
Perform a "negate" operation on all data vectors : x = -x.
Definition: TQpVar.cxx:273
TVectorD fCupIndex
Definition: TQpVar.h:74
virtual void Saxpy(TQpVar *b, Double_t alpha)
Perform a "saxpy" operation on all data vectors : x += alpha*y.
Definition: TQpVar.cxx:233
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:454
Element Min() const
return minimum vector element value
Definition: TVectorT.cxx:655
Element Norm1() const
Compute the 1-norm of the vector SUM{ |v[i]| }.
Definition: TVectorT.cxx:568
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:604
void AddSomeConstant(Element val, const TVectorT< Element > &select)
Add to vector elements as selected through array select the value val.
Definition: TVectorT.cxx:1285
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1239
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:295
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:350
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:545
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1364
Bool_t SomePositive(const TVectorT< Element > &select)
Check if vector elements as selected through array select are all positive.
Definition: TVectorT.cxx:1262
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:621
Element * GetMatrixArray()
Definition: TVectorT.h:78
double beta(double x, double y)
Calculates the beta function.
const Int_t n
Definition: legend1.C:16
void Add(RHist< DIMENSIONS, PRECISION, STAT_TO... > &to, const RHist< DIMENSIONS, PRECISION, STAT_FROM... > &from)
Add two histograms.
Definition: RHist.hxx:323