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