Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 "TQpDataSparse.h"
44
45////////////////////////////////////////////////////////////////////////////////
46///
47/// \class TQpDataSparse
48///
49/// Data for the sparse QP formulation
50///
51////////////////////////////////////////////////////////////////////////////////
52
54
55////////////////////////////////////////////////////////////////////////////////
56/// Constructor
57
59: TQpDataBase(nx,my,mz)
60{
64}
65
66
67////////////////////////////////////////////////////////////////////////////////
68/// Constructor
69
71 TVectorD &xlow_in,TVectorD &ixlow_in,
72 TVectorD &xupp_in,TVectorD &ixupp_in,
73 TMatrixDSparse &A_in, TVectorD &bA_in,
74 TMatrixDSparse &C_in,
75 TVectorD &clow_in,TVectorD &iclow_in,
76 TVectorD &cupp_in,TVectorD &icupp_in)
77{
78 fG .ResizeTo(c_in) ; fG = c_in;
79 fBa .ResizeTo(bA_in) ; fBa = bA_in;
80 fXloBound.ResizeTo(xlow_in) ; fXloBound = xlow_in;
81 fXloIndex.ResizeTo(ixlow_in); fXloIndex = ixlow_in;
82 fXupBound.ResizeTo(xupp_in) ; fXupBound = xupp_in;
83 fXupIndex.ResizeTo(ixupp_in); fXupIndex = ixupp_in;
84 fCloBound.ResizeTo(clow_in) ; fCloBound = clow_in;
85 fCloIndex.ResizeTo(iclow_in); fCloIndex = iclow_in;
86 fCupBound.ResizeTo(cupp_in) ; fCupBound = cupp_in;
87 fCupIndex.ResizeTo(icupp_in); fCupIndex = icupp_in;
88
89 fNx = fG.GetNrows();
90 fQ.Use(Q_in);
91
92 if (A_in.GetNrows() > 0) {
93 fA.Use(A_in);
94 fMy = fA.GetNrows();
95 } else
96 fMy = 0;
97
98 if (C_in.GetNrows()) {
99 fC.Use(C_in);
100 fMz = fC.GetNrows();
101 } else
102 fMz = 0;
103 // fQ.Print();
104 // fA.Print();
105 // fC.Print();
106 // printf("fNx: %d\n",fNx);
107 // printf("fMy: %d\n",fMy);
108 // printf("fMz: %d\n",fMz);
109}
110
111
112////////////////////////////////////////////////////////////////////////////////
113/// Copy constructor
114
116{
117 *this = another;
118}
119
120
121////////////////////////////////////////////////////////////////////////////////
122/// Allocate space for the appropriate number of non-zeros in the matrices
123
125{
126 fQ.SetSparseIndex(nnzQ);
127 fA.SetSparseIndex(nnzA);
128 fC.SetSparseIndex(nnzC);
129}
130
131
132////////////////////////////////////////////////////////////////////////////////
133/// calculate y = beta*y + alpha*(fQ*x)
134
136{
137 y *= beta;
138 if (fQ.GetNoElements() > 0)
139 y += alpha*(fQ*x);
140}
141
142
143////////////////////////////////////////////////////////////////////////////////
144/// calculate y = beta*y + alpha*(fA*x)
145
147{
148 y *= beta;
149 if (fA.GetNoElements() > 0)
150 y += alpha*(fA*x);
151}
152
153
154////////////////////////////////////////////////////////////////////////////////
155/// calculate y = beta*y + alpha*(fC*x)
156
158{
159 y *= beta;
160 if (fC.GetNoElements() > 0)
161 y += alpha*(fC*x);
162}
163
164
165////////////////////////////////////////////////////////////////////////////////
166/// calculate y = beta*y + alpha*(fA^T*x)
167
169{
170 y *= beta;
171 if (fA.GetNoElements() > 0)
173}
174
175
176////////////////////////////////////////////////////////////////////////////////
177/// calculate y = beta*y + alpha*(fC^T*x)
178
180{
181 y *= beta;
182 if (fC.GetNoElements() > 0)
184}
185
186
187////////////////////////////////////////////////////////////////////////////////
188/// Return the largest component of several vectors in the data class
189
191{
192 Double_t norm = 0.0;
193
194 Double_t componentNorm = fG.NormInf();
195 if (componentNorm > norm) norm = componentNorm;
196
197 TMatrixDSparse fQ_abs(fQ);
198 componentNorm = (fQ_abs.Abs()).Max();
199 if (componentNorm > norm) norm = componentNorm;
200
201 componentNorm = fBa.NormInf();
202 if (componentNorm > norm) norm = componentNorm;
203
204 TMatrixDSparse fA_abs(fA);
205 componentNorm = (fA_abs.Abs()).Max();
206 if (componentNorm > norm) norm = componentNorm;
207
208 TMatrixDSparse fC_abs(fC);
209 componentNorm = (fC_abs.Abs()).Max();
210 if (componentNorm > norm) norm = componentNorm;
211
213 componentNorm = fXloBound.NormInf();
214 if (componentNorm > norm) norm = componentNorm;
215
217 componentNorm = fXupBound.NormInf();
218 if (componentNorm > norm) norm = componentNorm;
219
221 componentNorm = fCloBound.NormInf();
222 if (componentNorm > norm) norm = componentNorm;
223
225 componentNorm = fCupBound.NormInf();
226 if (componentNorm > norm) norm = componentNorm;
227
228 return norm;
229}
230
231
232////////////////////////////////////////////////////////////////////////////////
233/// Print class members
234
235void TQpDataSparse::Print(Option_t * /*opt*/) const
236{
237 fQ.Print("Q");
238 fG.Print("c");
239
240 fXloBound.Print("xlow");
241 fXloIndex.Print("ixlow");
242
243 fXupBound.Print("xupp");
244 fXupIndex.Print("ixupp");
245
246 fA.Print("A");
247 fBa.Print("b");
248 fC.Print("C");
249
250 fCloBound.Print("clow");
251 fCloIndex.Print("iclow");
252
253 fCupBound.Print("cupp");
254 fCupIndex.Print("icupp");
255}
256
257
258////////////////////////////////////////////////////////////////////////////////
259/// Insert the Hessian Q into the matrix M at index (row,col) for the fundamental
260/// linear system
261
263{
264 m.SetSub(row,col,fQ);
265}
266
267
268////////////////////////////////////////////////////////////////////////////////
269/// Insert the constraint matrix A into the matrix M at index (row,col) for the fundamental
270/// linear system
271
273{
274 m.SetSub(row,col,fA);
275}
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// Insert the constraint matrix C into the matrix M at index (row,col) for the fundamental
280/// linear system
281
283{
284 m.SetSub(row,col,fC);
285}
286
287
288////////////////////////////////////////////////////////////////////////////////
289/// Return in vector dq the diagonal of matrix fQ
290
292{
293 const Int_t n = TMath::Min(fQ.GetNrows(),fQ.GetNcols());
294 dq.ResizeTo(n);
296}
297
298
299////////////////////////////////////////////////////////////////////////////////
300/// Return value of the objective function
301
303{
304 TVectorD tmp(fG);
305 this->Qmult(1.0,tmp,0.5,vars->fX);
306
307 return tmp*vars->fX;
308}
309
310
311////////////////////////////////////////////////////////////////////////////////
312/// Choose randomly a QP problem
313
315{
316 Double_t ix = 3074.20374;
317
318 TVectorD xdual(fNx);
320
321 TVectorD sprime(fMz);
323
324 fQ.RandomizePD(0.0,1.0,ix);
325 fA.Randomize(-10.0,10.0,ix);
326 fC.Randomize(-10.0,10.0,ix);
327 y .Randomize(-10.0,10.0,ix);
328
329 // fG = - fQ x + A^T y + C^T z + xdual
330
331 fG = xdual;
332 fG -= fQ*x;
333
336
337 fBa = fA*x;
338 s = fC*x;
339
340 // Now compute the real q = s-sprime
341 const TVectorD q = s-sprime;
342
343 // Adjust fCloBound and fCupBound appropriately
344 Add(fCloBound,1.0,q);
345 Add(fCupBound,1.0,q);
346
349}
350
351
352////////////////////////////////////////////////////////////////////////////////
353/// Assignment operator
354
356{
357 if (this != &source) {
359 fQ.ResizeTo(source.fQ); fQ = source.fQ;
360 fA.ResizeTo(source.fA); fA = source.fA;
361 fC.ResizeTo(source.fC); fC = source.fC;
362 }
363 return *this;
364}
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
float * q
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixTSparseDiag< Double_t > TMatrixDSparseDiag
virtual TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed)
Randomize matrix element values.
Int_t GetNrows() const
void Print(Option_t *name="") const override
Print the matrix as a table of elements.
Int_t GetNoElements() const
Int_t GetNcols() const
virtual TMatrixTBase< Element > & Abs()
Take an absolute value of a matrix, i.e. apply Abs() to each element.
TMatrixTBase< Element > & Randomize(Element alpha, Element beta, Double_t &seed) override
randomize matrix element values
TMatrixTSparse< Element > & SetSparseIndex(Int_t nelem_new)
Increase/decrease the number of non-zero elements to nelems_new.
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)
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t nr_nonzeros=-1) override
Set size of the matrix to nrows x ncols with nr_nonzeros non-zero entries if nr_nonzeros > 0 .
virtual TMatrixTSparse< Element > & RandomizePD(Element alpha, Element beta, Double_t &seed)
randomize matrix element values but keep matrix symmetric positive definite
Data for the general QP formulation.
Definition TQpDataBase.h:61
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.
Data for the sparse QP formulation.
Double_t ObjectiveValue(TQpVar *vars) override
Return value of the objective function.
Double_t DataNorm() override
Return the largest component of several vectors in the data class.
void PutAIntoAt(TMatrixDBase &M, Int_t row, Int_t col) override
Insert the constraint matrix A into the matrix M at index (row,col) for the fundamental linear system...
void Amult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fA*x)
void PutQIntoAt(TMatrixDBase &M, Int_t row, Int_t col) override
Insert the Hessian Q into the matrix M at index (row,col) for the fundamental linear system.
void DataRandom(TVectorD &x, TVectorD &y, TVectorD &z, TVectorD &s) override
Choose randomly a QP problem.
void Print(Option_t *opt="") const override
Print class members.
void PutCIntoAt(TMatrixDBase &M, Int_t row, Int_t col) override
Insert the constraint matrix C into the matrix M at index (row,col) for the fundamental linear system...
TQpDataSparse & operator=(const TQpDataSparse &source)
Assignment operator.
TMatrixDSparse fC
TMatrixDSparse fQ
void SetNonZeros(Int_t nnzQ, Int_t nnzA, Int_t nnzC)
Allocate space for the appropriate number of non-zeros in the matrices.
void Cmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fC*x)
void Qmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fQ*x)
void GetDiagonalOfQ(TVectorD &dQ) override
Return in vector dq the diagonal of matrix fQ.
TMatrixDSparse fA
void CTransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fC^T*x)
void ATransmult(Double_t beta, TVectorD &y, Double_t alpha, const TVectorD &x) override
calculate y = beta*y + alpha*(fA^T*x)
Class containing the variables for the general QP formulation.
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:603
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 > & SelectNonZeros(const TVectorT< Element > &select)
Keep only element as selected through array select non-zero.
Definition TVectorT.cxx:544
void Print(Option_t *option="") const override
Print the vector as a list of elements.
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
TMarker m
Definition textangle.C:8