Logo ROOT  
Reference Guide
TQpLinSolverBase.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////////////////////////////////////////////////////////////////////////////////
44///
45/// \class TQpLinSolverBase
46///
47/// Implementation of main solver for linear systems
48///
49////////////////////////////////////////////////////////////////////////////////
50
51#include "TQpLinSolverBase.h"
52
54
55////////////////////////////////////////////////////////////////////////////////
56/// Default constructor
57
59{
60 fNx = 0;
61 fMy = 0;
62 fMz = 0;
63 fNxup = 0;
64 fNxlo = 0;
65 fMcup = 0;
66 fMclo = 0;
67 fFactory = 0;
68}
69
70
71////////////////////////////////////////////////////////////////////////////////
72/// Constructor
73
75{
76 fFactory = factory;
77
78 fNx = data->fNx;
79 fMy = data->fMy;
80 fMz = data->fMz;
81
82 fXloIndex.ResizeTo(data->fXloIndex); fXloIndex = data->fXloIndex;
83 fXupIndex.ResizeTo(data->fXupIndex); fXupIndex = data->fXupIndex;
84 fCloIndex.ResizeTo(data->fCloIndex); fCloIndex = data->fCloIndex;
85 fCupIndex.ResizeTo(data->fCupIndex); fCupIndex = data->fCupIndex;
86
91
92 if (fNxup+fNxlo > 0) {
95 data->GetDiagonalOfQ(fDq);
96 }
99}
100
101
102////////////////////////////////////////////////////////////////////////////////
103/// Copy constructor
104
106 fFactory(another.fFactory)
107{
108 *this = another;
109}
110
111
112////////////////////////////////////////////////////////////////////////////////
113/// Sets up the matrix for the main linear system in "augmented system" form. The
114/// actual factorization is performed by a routine specific to either the sparse
115/// or dense case
116
118{
120
121 if (fNxlo+fNxup > 0) {
123 fDd = fDq;
124 }
126 vars->fT,vars->fLambda,vars->fU,vars->fPi,
127 vars->fV,vars->fGamma,vars->fW,vars->fPhi);
128 if (fNxlo+fNxup > 0) this->PutXDiagonal(fDd);
129
131 fNomegaInv *= -1.;
132
133 if (fMclo+fMcup > 0) this->PutZDiagonal(fNomegaInv);
134}
135
136
137////////////////////////////////////////////////////////////////////////////////
138/// Computes the diagonal matrices in the augmented system from the current set of variables
139
141 TVectorD &t, TVectorD &lambda,
142 TVectorD &u, TVectorD &pi,
144 TVectorD &w, TVectorD &phi)
145{
146 if (fNxup+fNxlo > 0) {
147 if (fNxlo > 0) AddElemDiv(dd,1.0,gamma,v,fXloIndex);
148 if (fNxup > 0) AddElemDiv(dd,1.0,phi ,w,fXupIndex);
149 }
150 omega.Zero();
151 if (fMclo > 0) AddElemDiv(omega,1.0,lambda,t,fCloIndex);
152 if (fMcup > 0) AddElemDiv(omega,1.0,pi, u,fCupIndex);
153}
154
155
156////////////////////////////////////////////////////////////////////////////////
157/// Solves the system for a given set of residuals. Assembles the right-hand side appropriate
158/// to the matrix factored in factor, solves the system using the factorization produced there,
159/// partitions the solution vector into step components, then recovers the step components
160/// eliminated during the block elimination that produced the augmented system form .
161
163{
165 R__ASSERT(res ->ValidNonZeroPattern());
166
167 (step->fX).ResizeTo(res->fRQ); step->fX = res->fRQ;
168 if (fNxlo > 0) {
169 TVectorD &vInvGamma = step->fV;
170 vInvGamma.ResizeTo(vars->fGamma); vInvGamma = vars->fGamma;
171 ElementDiv(vInvGamma,vars->fV,fXloIndex);
172
173 AddElemMult(step->fX,1.0,vInvGamma,res->fRv);
174 AddElemDiv (step->fX,1.0,res->fRgamma,vars->fV,fXloIndex);
175 }
176
177 if (fNxup > 0) {
178 TVectorD &wInvPhi = step->fW;
179 wInvPhi.ResizeTo(vars->fPhi); wInvPhi = vars->fPhi;
180 ElementDiv(wInvPhi,vars->fW,fXupIndex);
181
182 AddElemMult(step->fX,1.0,wInvPhi,res->fRw);
183 AddElemDiv (step->fX,-1.0,res->fRphi,vars->fW,fXupIndex);
184 }
185
186 // start by partially computing step->fS
187 (step->fS).ResizeTo(res->fRz); step->fS = res->fRz;
188 if (fMclo > 0) {
189 TVectorD &tInvLambda = step->fT;
190 tInvLambda.ResizeTo(vars->fLambda); tInvLambda = vars->fLambda;
191 ElementDiv(tInvLambda,vars->fT,fCloIndex);
192
193 AddElemMult(step->fS,1.0,tInvLambda,res->fRt);
194 AddElemDiv (step->fS,1.0,res->fRlambda,vars->fT,fCloIndex);
195 }
196
197 if (fMcup > 0) {
198 TVectorD &uInvPi = step->fU;
199
200 uInvPi.ResizeTo(vars->fPi); uInvPi = vars->fPi;
201 ElementDiv(uInvPi,vars->fU,fCupIndex);
202
203 AddElemMult(step->fS,1.0,uInvPi,res->fRu);
204 AddElemDiv (step->fS,-1.0,res->fRpi,vars->fU,fCupIndex);
205 }
206
207 (step->fY).ResizeTo(res->fRA); step->fY = res->fRA;
208 (step->fZ).ResizeTo(res->fRC); step->fZ = res->fRC;
209
210 if (fMclo > 0)
211 this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fLambda,prob);
212 else
213 this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fPi,prob);
214
215 if (fMclo > 0) {
216 (step->fT).ResizeTo(step->fS); step->fT = step->fS;
217 Add(step->fT,-1.0,res->fRt);
218 (step->fT).SelectNonZeros(fCloIndex);
219
220 (step->fLambda).ResizeTo(res->fRlambda); step->fLambda = res->fRlambda;
221 AddElemMult(step->fLambda,-1.0,vars->fLambda,step->fT);
222 ElementDiv(step->fLambda,vars->fT,fCloIndex);
223 }
224
225 if (fMcup > 0) {
226 (step->fU).ResizeTo(res->fRu); step->fU = res->fRu;
227 Add(step->fU,-1.0,step->fS);
228 (step->fU).SelectNonZeros(fCupIndex);
229
230 (step->fPi).ResizeTo(res->fRpi); step->fPi = res->fRpi;
231 AddElemMult(step->fPi,-1.0,vars->fPi,step->fU);
232 ElementDiv(step->fPi,vars->fU,fCupIndex);
233 }
234
235 if (fNxlo > 0) {
236 (step->fV).ResizeTo(step->fX); step->fV = step->fX;
237 Add(step->fV,-1.0,res->fRv);
238 (step->fV).SelectNonZeros(fXloIndex);
239
240 (step->fGamma).ResizeTo(res->fRgamma); step->fGamma = res->fRgamma;
241 AddElemMult(step->fGamma,-1.0,vars->fGamma,step->fV);
242 ElementDiv(step->fGamma,vars->fV,fXloIndex);
243 }
244
245 if (fNxup > 0) {
246 (step->fW).ResizeTo(res->fRw); step->fW = res->fRw;
247 Add(step->fW,-1.0,step->fX);
248 (step->fW).SelectNonZeros(fXupIndex);
249
250 (step->fPhi).ResizeTo(res->fRphi); step->fPhi = res->fRphi;
251 AddElemMult(step->fPhi,-1.0,vars->fPhi,step->fW);
252 ElementDiv(step->fPhi,vars->fW,fXupIndex);
253 }
255}
256
257
258////////////////////////////////////////////////////////////////////////////////
259/// Assemble right-hand side of augmented system and call SolveCompressed to solve it
260
262 TVectorD &stepz,TVectorD &steps,
263 TVectorD & /* ztemp */, TQpDataBase * /* prob */ )
264{
265 AddElemMult(stepz,-1.0,fNomegaInv,steps);
266 this->JoinRHS(fRhs,stepx,stepy,stepz);
267
268 this->SolveCompressed(fRhs);
269
270 this->SeparateVars(stepx,stepy,stepz,fRhs);
271
272 stepy *= -1.;
273 stepz *= -1.;
274
275 Add(steps,-1.0,stepz);
276 ElementMult(steps,fNomegaInv);
277 steps *= -1.;
278}
279
280
281////////////////////////////////////////////////////////////////////////////////
282/// Assembles a single vector object from three given vectors .
283/// rhs_out (output) final joined vector
284/// rhs1_in (input) first part of rhs
285/// rhs2_in (input) middle part of rhs
286/// rhs3_in (input) last part of rhs .
287
289 TVectorD &rhs2_in,TVectorD &rhs3_in)
290{
291 fFactory->JoinRHS(rhs_out,rhs1_in,rhs2_in,rhs3_in);
292}
293
294
295////////////////////////////////////////////////////////////////////////////////
296/// Extracts three component vectors from a given aggregated vector.
297/// vars_in (input) aggregated vector
298/// x_in (output) first part of vars
299/// y_in (output) middle part of vars
300/// z_in (output) last part of vars
301
303 TVectorD &z_in,TVectorD &vars_in)
304{
305 fFactory->SeparateVars(x_in,y_in,z_in,vars_in);
306}
307
308
309////////////////////////////////////////////////////////////////////////////////
310/// Assignment opeartor
311
313{
314 if (this != &source) {
315 TObject::operator=(source);
316
317 fNx = source.fNx;
318 fMy = source.fMy;
319 fMz = source.fMz;
320 fNxup = source.fNxup;
321 fNxlo = source.fNxlo;
322 fMcup = source.fMcup;
323 fMclo = source.fMclo;
324
326 fRhs .ResizeTo(source.fRhs); fRhs = source.fRhs;
327
328 fDd .ResizeTo(source.fDd); fDd = source.fDd;
329 fDq .ResizeTo(source.fDq); fDq = source.fDq;
330
335
336 // LM : copy also pointer data member
337 fFactory = source.fFactory;
338 }
339 return *this;
340}
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
ClassImp(TQpLinSolverBase)
Mother of all ROOT objects.
Definition: TObject.h:41
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:298
Data for the general QP formulation.
Definition: TQpDataBase.h:61
Implementation of main solver for linear systems.
virtual void JoinRHS(TVectorD &rhs, TVectorD &rhs1, TVectorD &rhs2, TVectorD &rhs3)
Assembles a single vector object from three given vectors .
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
TQpLinSolverBase & operator=(const TQpLinSolverBase &source)
Assignment opeartor.
virtual void SolveCompressed(TVectorD &rhs)=0
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
virtual void ComputeDiagonals(TVectorD &dd, TVectorD &omega, TVectorD &t, TVectorD &lambda, TVectorD &u, TVectorD &pi, TVectorD &v, TVectorD &gamma, TVectorD &w, TVectorD &phi)
Computes the diagonal matrices in the augmented system from the current set of variables.
virtual void PutXDiagonal(TVectorD &xdiag)=0
TQpLinSolverBase()
Default constructor.
TQpProbBase * fFactory
virtual void SeparateVars(TVectorD &vars1, TVectorD &vars2, TVectorD &vars3, TVectorD &vars)
Extracts three component vectors from a given aggregated vector.
virtual void SolveXYZS(TVectorD &stepx, TVectorD &stepy, TVectorD &stepz, TVectorD &steps, TVectorD &ztemp, TQpDataBase *data)
Assemble right-hand side of augmented system and call SolveCompressed to solve it.
virtual void PutZDiagonal(TVectorD &zdiag)=0
default general problem formulation:
Definition: TQpProbBase.h:89
virtual void JoinRHS(TVectorD &rhs_in, TVectorD &rhs1_in, TVectorD &rhs2_in, TVectorD &rhs3_in)=0
virtual void SeparateVars(TVectorD &x_in, TVectorD &y_in, TVectorD &z_in, TVectorD &vars_in)=0
The Residuals class calculates and stores the quantities that appear on the right-hand side of the li...
Definition: TQpResidual.h:62
TVectorD fRgamma
Definition: TQpResidual.h:96
TVectorD fRu
Definition: TQpResidual.h:95
TVectorD fRpi
Definition: TQpResidual.h:99
TVectorD fRC
Definition: TQpResidual.h:90
TVectorD fRQ
Definition: TQpResidual.h:88
TVectorD fRt
Definition: TQpResidual.h:94
TVectorD fRw
Definition: TQpResidual.h:93
TVectorD fRlambda
Definition: TQpResidual.h:98
TVectorD fRA
Definition: TQpResidual.h:89
TVectorD fRv
Definition: TQpResidual.h:92
TVectorD fRphi
Definition: TQpResidual.h:97
TVectorD fRz
Definition: TQpResidual.h:91
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
virtual Bool_t ValidNonZeroPattern()
Check that the variables conform to the non-zero indices.
Definition: TQpVar.cxx:740
TVectorD fLambda
Definition: TQpVar.h:103
TVectorD fGamma
Definition: TQpVar.h:100
TVectorD fT
Definition: TQpVar.h:102
TVectorD fPi
Definition: TQpVar.h:106
TVectorD fS
Definition: TQpVar.h:92
TVectorD fPhi
Definition: TQpVar.h:97
TVectorD fV
Definition: TQpVar.h:96
TVectorD fW
Definition: TQpVar.h:99
TVectorD fY
Definition: TQpVar.h:93
TVectorD fZ
Definition: TQpVar.h:94
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:294
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:620
TVectorT< Element > & Invert()
v[i] = 1/v[i]
Definition: TVectorT.cxx:522
double gamma(double x)
static constexpr double pi
TVectorT< Element > & AddElemDiv(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementDiv(source1,source2) .
Definition: TVectorT.cxx:1917
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
Definition: TMatrixT.cxx:2982
TVectorT< Element > & AddElemMult(TVectorT< Element > &target, Element scalar, const TVectorT< Element > &source1, const TVectorT< Element > &source2)
Modify addition: target += scalar * ElementMult(source1,source2) .
Definition: TVectorT.cxx:1844
TMatrixT< Element > & ElementDiv(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Divide target by the source, element-by-element.
Definition: TMatrixT.cxx:3022