Logo ROOT   6.13/01
Reference Guide
TQpDataDens.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 #include "Riostream.h"
44 #include "TQpDataDens.h"
45 
46 //////////////////////////////////////////////////////////////////////////
47 // //
48 // TQpDataDens //
49 // //
50 // Data for the dense QP formulation //
51 // //
52 //////////////////////////////////////////////////////////////////////////
53 
54 ClassImp(TQpDataDens);
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// Constructor
58 
59 TQpDataDens::TQpDataDens(Int_t nx,Int_t my,Int_t mz)
60 : TQpDataBase(nx,my,mz)
61 {
62  fQ.ResizeTo(fNx,fNx);
63  fA.ResizeTo(fMy,fNx);
64  fC.ResizeTo(fMz,fNx);
65 }
66 
67 
68 ////////////////////////////////////////////////////////////////////////////////
69 /// Constructor
70 
72  TVectorD &xlow_in,TVectorD &ixlow_in,
73  TVectorD &xupp_in,TVectorD &ixupp_in,
74  TMatrixD &A_in, TVectorD &bA_in,
75  TMatrixD &C_in,
76  TVectorD &clow_in,TVectorD &iclow_in,
77  TVectorD &cupp_in,TVectorD &icupp_in)
78 {
79  fG .ResizeTo(c_in) ; fG = c_in;
80  fBa .ResizeTo(bA_in) ; fBa = bA_in;
81  fXloBound.ResizeTo(xlow_in) ; fXloBound = xlow_in;
82  fXloIndex.ResizeTo(ixlow_in); fXloIndex = ixlow_in;
83  fXupBound.ResizeTo(xupp_in) ; fXupBound = xupp_in;
84  fXupIndex.ResizeTo(ixupp_in); fXupIndex = ixupp_in;
85  fCloBound.ResizeTo(clow_in) ; fCloBound = clow_in;
86  fCloIndex.ResizeTo(iclow_in); fCloIndex = iclow_in;
87  fCupBound.ResizeTo(cupp_in) ; fCupBound = cupp_in;
88  fCupIndex.ResizeTo(icupp_in); fCupIndex = icupp_in;
89 
90  fNx = fG.GetNrows();
91  fQ.Use(Q_in);
92 
93  if (A_in.GetNrows() > 0) {
94  fA.Use(A_in);
95  fMy = fA.GetNrows();
96  } else
97  fMy = 0;
98 
99  if (C_in.GetNrows() > 0) {
100  fC.Use(C_in);
101  fMz = fC.GetNrows();
102  } else
103  fMz = 0;
104 }
105 
106 
107 ////////////////////////////////////////////////////////////////////////////////
108 /// Copy constructor
109 
111 {
112  *this = another;
113 }
114 
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 /// calculate y = beta*y + alpha*(fQ*x)
118 
119 void TQpDataDens::Qmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x )
120 {
121  y *= beta;
122  if (fQ.GetNoElements() > 0)
123  y += alpha*(fQ*x);
124 }
125 
126 
127 ////////////////////////////////////////////////////////////////////////////////
128 /// calculate y = beta*y + alpha*(fA*x)
129 
130 void TQpDataDens::Amult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
131 {
132  y *= beta;
133  if (fA.GetNoElements() > 0)
134  y += alpha*(fA*x);
135 }
136 
137 
138 ////////////////////////////////////////////////////////////////////////////////
139 /// calculate y = beta*y + alpha*(fC*x)
140 
141 void TQpDataDens::Cmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
142 {
143  y *= beta;
144  if (fC.GetNoElements() > 0)
145  y += alpha*(fC*x);
146 }
147 
148 
149 ////////////////////////////////////////////////////////////////////////////////
150 /// calculate y = beta*y + alpha*(fA^T*x)
151 
152 void TQpDataDens::ATransmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
153 {
154  y *= beta;
155  if (fA.GetNoElements() > 0)
156  y += alpha*(TMatrixD(TMatrixD::kTransposed,fA)*x);
157 }
158 
159 
160 ////////////////////////////////////////////////////////////////////////////////
161 /// calculate y = beta*y + alpha*(fC^T*x)
162 
163 void TQpDataDens::CTransmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
164 {
165  y *= beta;
166  if (fC.GetNoElements() > 0)
167  y += alpha*(TMatrixD(TMatrixD::kTransposed,fC)*x);
168 }
169 
170 
171 ////////////////////////////////////////////////////////////////////////////////
172 /// Return the largest component of several vectors in the data class
173 
175 {
176  Double_t norm = 0.0;
177 
178  Double_t componentNorm = fG.NormInf();
179  if (componentNorm > norm) norm = componentNorm;
180 
181  TMatrixDSym fQ_abs(fQ);
182  componentNorm = (fQ_abs.Abs()).Max();
183  if (componentNorm > norm) norm = componentNorm;
184 
185  componentNorm = fBa.NormInf();
186  if (componentNorm > norm) norm = componentNorm;
187 
188  TMatrixD fA_abs(fQ);
189  componentNorm = (fA_abs.Abs()).Max();
190  if (componentNorm > norm) norm = componentNorm;
191 
192  TMatrixD fC_abs(fQ);
193  componentNorm = (fC_abs.Abs()).Max();
194  if (componentNorm > norm) norm = componentNorm;
195 
197  componentNorm = fXloBound.NormInf();
198  if (componentNorm > norm) norm = componentNorm;
199 
201  componentNorm = fXupBound.NormInf();
202  if (componentNorm > norm) norm = componentNorm;
203 
205  componentNorm = fCloBound.NormInf();
206  if (componentNorm > norm) norm = componentNorm;
207 
209  componentNorm = fCupBound.NormInf();
210  if (componentNorm > norm) norm = componentNorm;
211 
212  return norm;
213 }
214 
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Print all class members
218 
219 void TQpDataDens::Print(Option_t * /*opt*/) const
220 {
221  fQ.Print("Q");
222  fG.Print("c");
223 
224  fXloBound.Print("xlow");
225  fXloIndex.Print("ixlow");
226 
227  fXupBound.Print("xupp");
228  fXupIndex.Print("ixupp");
229 
230  fA.Print("A");
231  fBa.Print("b");
232  fC.Print("C");
233 
234  fCloBound.Print("clow");
235  fCloIndex.Print("iclow");
236 
237  fCupBound.Print("cupp");
238  fCupIndex.Print("icupp");
239 }
240 
241 
242 ////////////////////////////////////////////////////////////////////////////////
243 /// Insert the Hessian Q into the matrix M at index (row,col) for the fundamental
244 /// linear system
245 
246 void TQpDataDens::PutQIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
247 {
248  m.SetSub(row,col,fQ);
249 }
250 
251 
252 ////////////////////////////////////////////////////////////////////////////////
253 /// Insert the constraint matrix A into the matrix M at index (row,col) for the fundamental
254 /// linear system
255 
256 void TQpDataDens::PutAIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
257 {
258  m.SetSub(row,col,fA);
259 }
260 
261 
262 ////////////////////////////////////////////////////////////////////////////////
263 /// Insert the constraint matrix C into the matrix M at index (row,col) for the fundamental
264 /// linear system
265 
266 void TQpDataDens::PutCIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
267 {
268  m.SetSub(row,col,fC);
269 }
270 
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Return in vector dq the diagonal of matrix fQ (Quadratic part of Objective function)
274 
276 {
277  const Int_t n = TMath::Min(fQ.GetNrows(),fQ.GetNcols());
278  dq.ResizeTo(n);
279  dq = TMatrixDDiag(fQ);
280 }
281 
282 
283 ////////////////////////////////////////////////////////////////////////////////
284 /// Return value of the objective function
285 
287 {
288  TVectorD tmp(fG);
289  this->Qmult(1.0,tmp,0.5,vars->fX);
290 
291  return tmp*vars->fX;
292 }
293 
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// Choose randomly a QP problem
297 
299 {
300  Double_t ix = 3074.20374;
301 
302  TVectorD xdual(fNx);
304  TVectorD sprime(fMz);
306 
307  fQ.RandomizePD(0.0,1.0,ix);
308  fA.Randomize(-10.0,10.0,ix);
309  fC.Randomize(-10.0,10.0,ix);
310  y .Randomize(-10.0,10.0,ix);
311 
312  fG = xdual;
313  fG -= fQ*x;
314 
317 
318  fBa = fA*x;
319  s = fC*x;
320 
321  // Now compute the real q = s-sprime
322  const TVectorD q = s-sprime;
323 
324  // Adjust fCloBound and fCupBound appropriately
325  Add(fCloBound,1.0,q);
326  Add(fCupBound,1.0,q);
327 
330 }
331 
332 
333 ////////////////////////////////////////////////////////////////////////////////
334 /// Assignment operator
335 
337 {
338  if (this != &source) {
339  TQpDataBase::operator=(source);
340  fQ.ResizeTo(source.fQ); fQ = source.fQ;
341  fA.ResizeTo(source.fA); fA = source.fA;
342  fC.ResizeTo(source.fC); fC = source.fC;
343  }
344  return *this;
345 }
TVectorD fXupBound
Definition: TQpDataBase.h:79
virtual void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fQ*x)
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
TQpDataBase & operator=(const TQpDataBase &source)
Assignment operator.
TVectorD fXloIndex
Definition: TQpDataBase.h:82
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
TVectorD fG
Definition: TQpDataBase.h:77
virtual void PutQIntoAt(TMatrixDBase &M, Int_t row, Int_t col)
Insert the Hessian Q into the matrix M at index (row,col) for the fundamental linear system...
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixT.cxx:2933
Int_t GetNcols() const
Definition: TMatrixTBase.h:125
virtual void Print(Option_t *opt="") const
Print all class members.
Int_t GetNrows() const
Definition: TVectorT.h:75
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
TMatrixD fC
Definition: TQpDataDens.h:70
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1361
TVectorD fX
Definition: TQpVar.h:91
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1210
you should not use this method at all Int_t y
Definition: TRolke.cxx:630
TVectorD fCupIndex
Definition: TQpDataBase.h:84
double beta(double x, double y)
Calculates the beta function.
virtual void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC*x)
virtual void PutCIntoAt(TMatrixDBase &M, Int_t row, Int_t col)
Insert the constraint matrix C into the matrix M at index (row,col) for the fundamental linear system...
void Randomize(Element alpha, Element beta, Double_t &seed)
randomize vector elements value
Definition: TVectorT.cxx:1303
TMatrixDSym fQ
Definition: TQpDataDens.h:68
Int_t GetNoElements() const
Definition: TMatrixTBase.h:126
virtual void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA*x)
TVectorD fCloBound
Definition: TQpDataBase.h:85
virtual Double_t ObjectiveValue(TQpVar *vars)
Return value of the objective function.
TMatrixD fA
Definition: TQpDataDens.h:69
virtual void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC^T*x)
virtual void PutAIntoAt(TMatrixDBase &M, Int_t row, Int_t col)
Insert the constraint matrix A into the matrix M at index (row,col) for the fundamental linear system...
virtual void DataRandom(TVectorD &x, TVectorD &y, TVectorD &z, TVectorD &s)
Choose randomly a QP problem.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t Int_t m
Definition: TRolke.cxx:637
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
TQpDataDens & operator=(const TQpDataDens &source)
Assignment operator.
virtual void GetDiagonalOfQ(TVectorD &dQ)
Return in vector dq the diagonal of matrix fQ (Quadratic part of Objective function) ...
virtual TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source)=0
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
Randomize matrix element values.
TVectorD fXloBound
Definition: TQpDataBase.h:81
TMatrixTDiag< Double_t > TMatrixDDiag
TVectorD fXupIndex
Definition: TQpDataBase.h:80
virtual void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA^T*x)
Linear Algebra Package.
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:542
TVectorD fBa
Definition: TQpDataBase.h:78
virtual Double_t DataNorm()
Return the largest component of several vectors in the data class.
virtual TMatrixTSym< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
static void RandomlyChooseBoundedVariables(TVectorD &x, TVectorD &dualx, TVectorD &blx, TVectorD &ixlow, TVectorD &bux, TVectorD &ixupp, Double_t &ix, Double_t percentLowerOnly, Double_t percentUpperOnly, Double_t percentBound)
Randomly choose x and its boundaries.
void Print(Option_t *name="") const
Print the matrix as a table of elements.
TMatrixT< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Element *data)
Use the array data to fill the matrix ([row_lwb..row_upb] x [col_lwb..col_upb])
Definition: TMatrixT.cxx:1053
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Definition: TQpVar.h:59
TVectorD fCupBound
Definition: TQpDataBase.h:83
TVectorD fCloIndex
Definition: TQpDataBase.h:86
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1236
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.