Logo ROOT   6.16/01
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{
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
236void 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);
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) {
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}
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:363
#define R__ASSERT(e)
Definition: TError.h:96
float * q
Definition: THbookFile.cxx:87
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixTSparseDiag< Double_t > TMatrixDSparseDiag
Linear Algebra Package.
Definition: TMatrixTBase.h:85
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Int_t GetNoElements() const
Definition: TMatrixTBase.h:128
Int_t GetNcols() const
Definition: TMatrixTBase.h:127
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
TMatrixTSparse< Element > & SetSparseIndex(Int_t nelem_new)
Increase/decrease the number of non-zero elements to nelems_new.
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 .
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)
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
TQpDataBase & operator=(const TQpDataBase &source)
Assignment operator.
TVectorD fXupIndex
Definition: TQpDataBase.h:80
TVectorD fXloBound
Definition: TQpDataBase.h:81
TVectorD fCloIndex
Definition: TQpDataBase.h:86
TVectorD fCupIndex
Definition: TQpDataBase.h:84
TVectorD fG
Definition: TQpDataBase.h:77
TVectorD fXupBound
Definition: TQpDataBase.h:79
TVectorD fXloIndex
Definition: TQpDataBase.h:82
TVectorD fCupBound
Definition: TQpDataBase.h:83
TVectorD fCloBound
Definition: TQpDataBase.h:85
TVectorD fBa
Definition: TQpDataBase.h:78
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.
virtual void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC^T*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...
virtual void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fQ*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 Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA*x)
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.
virtual void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fC*x)
virtual void GetDiagonalOfQ(TVectorD &dQ)
Return in vector dq the diagonal of matrix fQ.
TQpDataSparse & operator=(const TQpDataSparse &source)
Assignment operator.
TMatrixDSparse fC
Definition: TQpDataSparse.h:67
TMatrixDSparse fQ
Definition: TQpDataSparse.h:65
void SetNonZeros(Int_t nnzQ, Int_t nnzA, Int_t nnzC)
Allocate space for the appropriate number of non-zeros in the matrices.
virtual Double_t DataNorm()
Return the largest component of several vectors in the data class.
TMatrixDSparse fA
Definition: TQpDataSparse.h:66
virtual void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x)
calculate y = beta*y + alpha*(fA^T*x)
virtual Double_t ObjectiveValue(TQpVar *vars)
Return value of the objective function.
virtual void Print(Option_t *opt="") const
Print class members.
virtual void DataRandom(TVectorD &x, TVectorD &y, TVectorD &z, TVectorD &s)
Choose randomly a QP problem.
Definition: TQpVar.h:60
TVectorD fX
Definition: TQpVar.h:91
Element NormInf() const
Compute the infinity-norm of the vector MAX{ |v[i]| }.
Definition: TVectorT.cxx:601
Bool_t MatchesNonZeroPattern(const TVectorT< Element > &select)
Check if vector elements as selected through array select are non-zero.
Definition: TVectorT.cxx:1236
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
Int_t GetNrows() const
Definition: TVectorT.h:75
TVectorT< Element > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition: TVectorT.cxx:542
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1361
double beta(double x, double y)
Calculates the beta function.
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
void Add(RHist< DIMENSIONS, PRECISION_TO, STAT_TO... > &to, const RHist< DIMENSIONS, PRECISION_FROM, STAT_FROM... > &from)
Add two histograms.
Definition: RHist.hxx:310
static constexpr double s
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
auto * m
Definition: textangle.C:8