Logo ROOT   6.10/09
Reference Guide
TQpDataSparse.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 "TQpDataSparse.h"
45 
46 //////////////////////////////////////////////////////////////////////////
47 // //
48 // TQpDataSparse //
49 // //
50 // Data for the sparse QP formulation //
51 // //
52 //////////////////////////////////////////////////////////////////////////
53 
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// Constructor
58 
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  TMatrixDSparse &A_in, TVectorD &bA_in,
75  TMatrixDSparse &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()) {
100  fC.Use(C_in);
101  fMz = fC.GetNrows();
102  } else
103  fMz = 0;
104  // fQ.Print();
105  // fA.Print();
106  // fC.Print();
107  // printf("fNx: %d\n",fNx);
108  // printf("fMy: %d\n",fMy);
109  // printf("fMz: %d\n",fMz);
110 }
111 
112 
113 ////////////////////////////////////////////////////////////////////////////////
114 /// Copy constructor
115 
117 {
118  *this = another;
119 }
120 
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 /// Allocate space for the appropriate number of non-zeros in the matrices
124 
126 {
127  fQ.SetSparseIndex(nnzQ);
128  fA.SetSparseIndex(nnzA);
129  fC.SetSparseIndex(nnzC);
130 }
131 
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// calculate y = beta*y + alpha*(fQ*x)
135 
137 {
138  y *= beta;
139  if (fQ.GetNoElements() > 0)
140  y += alpha*(fQ*x);
141 }
142 
143 
144 ////////////////////////////////////////////////////////////////////////////////
145 /// calculate y = beta*y + alpha*(fA*x)
146 
148 {
149  y *= beta;
150  if (fA.GetNoElements() > 0)
151  y += alpha*(fA*x);
152 }
153 
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 /// calculate y = beta*y + alpha*(fC*x)
157 
159 {
160  y *= beta;
161  if (fC.GetNoElements() > 0)
162  y += alpha*(fC*x);
163 }
164 
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 /// calculate y = beta*y + alpha*(fA^T*x)
168 
170 {
171  y *= beta;
172  if (fA.GetNoElements() > 0)
174 }
175 
176 
177 ////////////////////////////////////////////////////////////////////////////////
178 /// calculate y = beta*y + alpha*(fC^T*x)
179 
181 {
182  y *= beta;
183  if (fC.GetNoElements() > 0)
185 }
186 
187 
188 ////////////////////////////////////////////////////////////////////////////////
189 /// Return the largest component of several vectors in the data class
190 
192 {
193  Double_t norm = 0.0;
194 
195  Double_t componentNorm = fG.NormInf();
196  if (componentNorm > norm) norm = componentNorm;
197 
198  TMatrixDSparse fQ_abs(fQ);
199  componentNorm = (fQ_abs.Abs()).Max();
200  if (componentNorm > norm) norm = componentNorm;
201 
202  componentNorm = fBa.NormInf();
203  if (componentNorm > norm) norm = componentNorm;
204 
205  TMatrixDSparse fA_abs(fA);
206  componentNorm = (fA_abs.Abs()).Max();
207  if (componentNorm > norm) norm = componentNorm;
208 
209  TMatrixDSparse fC_abs(fC);
210  componentNorm = (fC_abs.Abs()).Max();
211  if (componentNorm > norm) norm = componentNorm;
212 
214  componentNorm = fXloBound.NormInf();
215  if (componentNorm > norm) norm = componentNorm;
216 
218  componentNorm = fXupBound.NormInf();
219  if (componentNorm > norm) norm = componentNorm;
220 
222  componentNorm = fCloBound.NormInf();
223  if (componentNorm > norm) norm = componentNorm;
224 
226  componentNorm = fCupBound.NormInf();
227  if (componentNorm > norm) norm = componentNorm;
228 
229  return norm;
230 }
231 
232 
233 ////////////////////////////////////////////////////////////////////////////////
234 /// Print class members
235 
236 void TQpDataSparse::Print(Option_t * /*opt*/) const
237 {
238  fQ.Print("Q");
239  fG.Print("c");
240 
241  fXloBound.Print("xlow");
242  fXloIndex.Print("ixlow");
243 
244  fXupBound.Print("xupp");
245  fXupIndex.Print("ixupp");
246 
247  fA.Print("A");
248  fBa.Print("b");
249  fC.Print("C");
250 
251  fCloBound.Print("clow");
252  fCloIndex.Print("iclow");
253 
254  fCupBound.Print("cupp");
255  fCupIndex.Print("icupp");
256 }
257 
258 
259 ////////////////////////////////////////////////////////////////////////////////
260 /// Insert the Hessian Q into the matrix M at index (row,col) for the fundamental
261 /// linear system
262 
264 {
265  m.SetSub(row,col,fQ);
266 }
267 
268 
269 ////////////////////////////////////////////////////////////////////////////////
270 /// Insert the constraint matrix A into the matrix M at index (row,col) for the fundamental
271 /// linear system
272 
274 {
275  m.SetSub(row,col,fA);
276 }
277 
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// Insert the constraint matrix C into the matrix M at index (row,col) for the fundamental
281 /// linear system
282 
284 {
285  m.SetSub(row,col,fC);
286 }
287 
288 
289 ////////////////////////////////////////////////////////////////////////////////
290 /// Return in vector dq the diagonal of matrix fQ
291 
293 {
294  const Int_t n = TMath::Min(fQ.GetNrows(),fQ.GetNcols());
295  dq.ResizeTo(n);
296  dq = TMatrixDSparseDiag(fQ);
297 }
298 
299 
300 ////////////////////////////////////////////////////////////////////////////////
301 /// Return value of the objective function
302 
304 {
305  TVectorD tmp(fG);
306  this->Qmult(1.0,tmp,0.5,vars->fX);
307 
308  return tmp*vars->fX;
309 }
310 
311 
312 ////////////////////////////////////////////////////////////////////////////////
313 /// Choose randomly a QP problem
314 
316 {
317  Double_t ix = 3074.20374;
318 
319  TVectorD xdual(fNx);
321 
322  TVectorD sprime(fMz);
324 
325  fQ.RandomizePD(0.0,1.0,ix);
326  fA.Randomize(-10.0,10.0,ix);
327  fC.Randomize(-10.0,10.0,ix);
328  y .Randomize(-10.0,10.0,ix);
329 
330  // fG = - fQ x + A^T y + C^T z + xdual
331 
332  fG = xdual;
333  fG -= fQ*x;
334 
337 
338  fBa = fA*x;
339  s = fC*x;
340 
341  // Now compute the real q = s-sprime
342  const TVectorD q = s-sprime;
343 
344  // Adjust fCloBound and fCupBound appropriately
345  Add(fCloBound,1.0,q);
346  Add(fCupBound,1.0,q);
347 
350 }
351 
352 
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Assignment operator
355 
357 {
358  if (this != &source) {
359  TQpDataBase::operator=(source);
360  fQ.ResizeTo(source.fQ); fQ = source.fQ;
361  fA.ResizeTo(source.fA); fA = source.fA;
362  fC.ResizeTo(source.fC); fC = source.fC;
363  }
364  return *this;
365 }
const int nx
Definition: kalman.C:16
TVectorD fXupBound
Definition: TQpDataBase.h:79
virtual void Print(Option_t *opt="") const
Print class members.
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 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...
TVectorD fG
Definition: TQpDataBase.h:77
const char Option_t
Definition: RtypesCore.h:62
Int_t GetNcols() const
Definition: TMatrixTBase.h:125
TMatrixDSparse fQ
Definition: TQpDataSparse.h:65
#define R__ASSERT(e)
Definition: TError.h:96
Int_t GetNrows() const
Definition: TVectorT.h:75
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1361
TVectorD fX
Definition: TQpVar.h:91
TMatrixDSparse fC
Definition: TQpDataSparse.h:67
TVectorD fCupIndex
Definition: TQpDataBase.h:84
virtual void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA^T*x)
double beta(double x, double y)
Calculates the beta function.
void Randomize(Element alpha, Element beta, Double_t &seed)
randomize vector elements value
Definition: TVectorT.cxx:1303
Int_t GetNoElements() const
Definition: TMatrixTBase.h:126
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...
Double_t x[n]
Definition: legend1.C:17
TVectorD fCloBound
Definition: TQpDataBase.h:85
virtual void DataRandom(TVectorD &x, TVectorD &y, TVectorD &z, TVectorD &s)
Choose randomly a QP problem.
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
randomize matrix element values
virtual TMatrixTSparse< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
TMatrixTSparse< Double_t > TMatrixDSparse
void SetNonZeros(Int_t nnzQ, Int_t nnzA, Int_t nnzC)
Allocate space for the appropriate number of non-zeros in the matrices.
virtual TMatrixTBase< Element > & SetSub(Int_t row_lwb, Int_t col_lwb, const TMatrixTBase< Element > &source)=0
virtual Double_t ObjectiveValue(TQpVar *vars)
Return value of the objective function.
TMarker * m
Definition: textangle.C:8
virtual Double_t DataNorm()
Return the largest component of several vectors in the data class.
TMatrixTSparse< Element > & SetSparseIndex(Int_t nelem_new)
Increase/decrease the number of non-zero elements to nelems_new.
virtual void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA*x)
TVectorD fXloBound
Definition: TQpDataBase.h:81
TVectorD fXupIndex
Definition: TQpDataBase.h:80
Linear Algebra Package.
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
void Add(THist< DIMENSIONS, PRECISION_TO, STAT_TO... > &to, const THist< DIMENSIONS, PRECISION_FROM, STAT_FROM... > &from)
Add two histograms.
Definition: THist.hxx:336
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
virtual void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC^T*x)
TMatrixTSparse< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t nr_nonzeros, Int_t *pRowIndex, Int_t *pColIndex, Element *pData)
TVectorD fBa
Definition: TQpDataBase.h:78
TMatrixDSparse fA
Definition: TQpDataSparse.h:66
Double_t y[n]
Definition: legend1.C:17
virtual void GetDiagonalOfQ(TVectorD &dQ)
Return in vector dq the diagonal of matrix fQ.
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.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
TQpDataSparse & operator=(const TQpDataSparse &source)
Assignment operator.
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...
Definition: TQpVar.h:59
TVectorD fCupBound
Definition: TQpDataBase.h:83
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1)
Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries if nr_nonzeros > 0 ...
virtual void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fQ*x)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
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
float * q
Definition: THbookFile.cxx:87
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
const Int_t n
Definition: legend1.C:16
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
TMatrixTSparseDiag< Double_t > TMatrixDSparseDiag
virtual void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC*x)