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 
86 
91 
92  if (fNxup+fNxlo > 0) {
93  fDd.ResizeTo(fNx);
94  fDq.ResizeTo(fNx);
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) {
122  fDd.ResizeTo(fDq);
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 
130  fNomegaInv.Invert();
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 
331  fXupIndex .ResizeTo(source.fXupIndex); fXupIndex = source.fXupIndex;
332  fCupIndex .ResizeTo(source.fCupIndex); fCupIndex = source.fCupIndex;
333  fXloIndex .ResizeTo(source.fXloIndex); fXloIndex = source.fXloIndex;
334  fCloIndex .ResizeTo(source.fCloIndex); fCloIndex = source.fCloIndex;
335 
336  // LM : copy also pointer data member
337  fFactory = source.fFactory;
338  }
339  return *this;
340 }
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
TQpLinSolverBase::fNx
Int_t fNx
Definition: TQpLinSolverBase.h:74
AddElemDiv
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
TQpResidual::fRC
TVectorD fRC
Definition: TQpResidual.h:125
TQpResidual::fRpi
TVectorD fRpi
Definition: TQpResidual.h:134
TQpVar::fZ
TVectorD fZ
Definition: TQpVar.h:129
TQpLinSolverBase::fMy
Int_t fMy
Definition: TQpLinSolverBase.h:75
TQpResidual::fRt
TVectorD fRt
Definition: TQpResidual.h:129
TQpVar::fS
TVectorD fS
Definition: TQpVar.h:127
TQpLinSolverBase.h
TQpVar::fT
TVectorD fT
Definition: TQpVar.h:137
TQpLinSolverBase::fXupIndex
TVectorD fXupIndex
Definition: TQpLinSolverBase.h:81
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
TQpLinSolverBase::PutZDiagonal
virtual void PutZDiagonal(TVectorD &zdiag)=0
TQpLinSolverBase::TQpLinSolverBase
TQpLinSolverBase()
Default constructor.
Definition: TQpLinSolverBase.cxx:58
TQpDataBase::fXloIndex
TVectorD fXloIndex
Definition: TQpDataBase.h:117
TQpVar::fV
TVectorD fV
Definition: TQpVar.h:131
TQpDataBase::GetDiagonalOfQ
virtual void GetDiagonalOfQ(TVectorD &dQ)=0
TQpResidual::fRgamma
TVectorD fRgamma
Definition: TQpResidual.h:131
TQpLinSolverBase::SeparateVars
virtual void SeparateVars(TVectorD &vars1, TVectorD &vars2, TVectorD &vars3, TVectorD &vars)
Extracts three component vectors from a given aggregated vector.
Definition: TQpLinSolverBase.cxx:302
TQpLinSolverBase::PutXDiagonal
virtual void PutXDiagonal(TVectorD &xdiag)=0
TVectorT::Zero
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
TQpDataBase
Definition: TQpDataBase.h:60
TVectorT::ResizeTo
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:294
TQpDataBase::fMy
Int_t fMy
Definition: TQpDataBase.h:109
TQpLinSolverBase::fXloIndex
TVectorD fXloIndex
Definition: TQpLinSolverBase.h:83
AddElemMult
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
TQpResidual::fRlambda
TVectorD fRlambda
Definition: TQpResidual.h:133
v
@ v
Definition: rootcling_impl.cxx:3635
TQpLinSolverBase::fDq
TVectorD fDq
Definition: TQpLinSolverBase.h:79
TQpResidual::fRv
TVectorD fRv
Definition: TQpResidual.h:127
TQpLinSolverBase::fDd
TVectorD fDd
Definition: TQpLinSolverBase.h:78
ElementMult
TMatrixT< Element > & ElementMult(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Multiply target by the source, element-by-element.
Definition: TMatrixT.cxx:2982
TQpVar::ValidNonZeroPattern
virtual Bool_t ValidNonZeroPattern()
Check that the variables conform to the non-zero indices.
Definition: TQpVar.cxx:740
TQpLinSolverBase::Solve
virtual void Solve(TQpDataBase *prob, TQpVar *vars, TQpResidual *resids, TQpVar *step)
Solves the system for a given set of residuals.
Definition: TQpLinSolverBase.cxx:162
TQpDataBase::fCupIndex
TVectorD fCupIndex
Definition: TQpDataBase.h:119
TQpResidual::fRz
TVectorD fRz
Definition: TQpResidual.h:126
TQpResidual::fRA
TVectorD fRA
Definition: TQpResidual.h:124
TQpLinSolverBase
Definition: TQpLinSolverBase.h:66
TQpVar::fW
TVectorD fW
Definition: TQpVar.h:134
TQpVar::fLambda
TVectorD fLambda
Definition: TQpVar.h:138
TQpLinSolverBase::JoinRHS
virtual void JoinRHS(TVectorD &rhs, TVectorD &rhs1, TVectorD &rhs2, TVectorD &rhs3)
Assembles a single vector object from three given vectors .
Definition: TQpLinSolverBase.cxx:288
TQpLinSolverBase::fMclo
Int_t fMclo
Definition: TQpLinSolverBase.h:89
TQpVar::fPhi
TVectorD fPhi
Definition: TQpVar.h:132
TQpLinSolverBase::fMz
Int_t fMz
Definition: TQpLinSolverBase.h:76
TQpLinSolverBase::fCupIndex
TVectorD fCupIndex
Definition: TQpLinSolverBase.h:82
TQpLinSolverBase::fCloIndex
TVectorD fCloIndex
Definition: TQpLinSolverBase.h:84
TQpResidual
Definition: TQpResidual.h:61
TQpVar::fPi
TVectorD fPi
Definition: TQpVar.h:141
TQpLinSolverBase::fFactory
TQpProbBase * fFactory
Definition: TQpLinSolverBase.h:91
TQpProbBase::SeparateVars
virtual void SeparateVars(TVectorD &x_in, TVectorD &y_in, TVectorD &z_in, TVectorD &vars_in)=0
TVectorT::Invert
TVectorT< Element > & Invert()
v[i] = 1/v[i]
Definition: TVectorT.cxx:522
TQpVar
Definition: TQpVar.h:59
TQpLinSolverBase::fMcup
Int_t fMcup
Definition: TQpLinSolverBase.h:88
TQpLinSolverBase::fNomegaInv
TVectorD fNomegaInv
Definition: TQpLinSolverBase.h:71
TVectorT::NonZeros
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:620
TQpDataBase::fNx
Int_t fNx
Definition: TQpDataBase.h:108
TQpProbBase::JoinRHS
virtual void JoinRHS(TVectorD &rhs_in, TVectorD &rhs1_in, TVectorD &rhs2_in, TVectorD &rhs3_in)=0
TVectorT
Definition: TMatrixTBase.h:78
TQpResidual::fRu
TVectorD fRu
Definition: TQpResidual.h:130
TQpLinSolverBase::SolveXYZS
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.
Definition: TQpLinSolverBase.cxx:261
TQpProbBase
Definition: TQpProbBase.h:88
TQpLinSolverBase::fNxup
Int_t fNxup
Definition: TQpLinSolverBase.h:86
TObject::operator=
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition: TObject.h:268
R__ASSERT
#define R__ASSERT(e)
Definition: TError.h:120
TQpDataBase::fMz
Int_t fMz
Definition: TQpDataBase.h:110
TGeant4Unit::pi
static constexpr double pi
Definition: TGeant4SystemOfUnits.h:73
TQpVar::fU
TVectorD fU
Definition: TQpVar.h:140
TQpDataBase::fXupIndex
TVectorD fXupIndex
Definition: TQpDataBase.h:115
TQpVar::fGamma
TVectorD fGamma
Definition: TQpVar.h:135
TObject
Definition: TObject.h:37
TQpDataBase::fCloIndex
TVectorD fCloIndex
Definition: TQpDataBase.h:121
TQpLinSolverBase::SolveCompressed
virtual void SolveCompressed(TVectorD &rhs)=0
TQpResidual::fRw
TVectorD fRw
Definition: TQpResidual.h:128
TQpLinSolverBase::operator=
TQpLinSolverBase & operator=(const TQpLinSolverBase &source)
Assignment opeartor.
Definition: TQpLinSolverBase.cxx:312
TQpLinSolverBase::ComputeDiagonals
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.
Definition: TQpLinSolverBase.cxx:140
TQpVar::fY
TVectorD fY
Definition: TQpVar.h:128
TQpLinSolverBase::fRhs
TVectorD fRhs
Definition: TQpLinSolverBase.h:72
TQpLinSolverBase::Factor
virtual void Factor(TQpDataBase *prob, TQpVar *vars)
Sets up the matrix for the main linear system in "augmented system" form.
Definition: TQpLinSolverBase.cxx:117
TQpResidual::fRQ
TVectorD fRQ
Definition: TQpResidual.h:123
ROOT::Experimental::Add
void Add(RHist< DIMENSIONS, PRECISION, STAT_TO... > &to, const RHist< DIMENSIONS, PRECISION, STAT_FROM... > &from)
Add two histograms.
Definition: RHist.hxx:335
ElementDiv
TMatrixT< Element > & ElementDiv(TMatrixT< Element > &target, const TMatrixT< Element > &source)
Divide target by the source, element-by-element.
Definition: TMatrixT.cxx:3022
TQpVar::fX
TVectorD fX
Definition: TQpVar.h:126
TQpResidual::fRphi
TVectorD fRphi
Definition: TQpResidual.h:132
TQpLinSolverBase::fNxlo
Int_t fNxlo
Definition: TQpLinSolverBase.h:87