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