MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
BiCGStab Class Reference

#include <bicg.h>

Inheritance diagram for BiCGStab:
Collaboration diagram for BiCGStab:

Public Member Functions

 BiCGStab (const Preconditioner::PrecondType PType, const integer iPStep, doublereal ITol, integer MaxIt, doublereal etaMx, doublereal T, const NonlinearSolverOptions &options)
 
 ~BiCGStab (void)
 
virtual void Solve (const NonlinearProblem *NLP, Solver *pS, const integer iMaxIter, const doublereal &Tol, integer &iIterCnt, doublereal &dErr, const doublereal &SolTol, doublereal &dSolErr)
 
- Public Member Functions inherited from MatrixFreeSolver
 MatrixFreeSolver (const Preconditioner::PrecondType PType, const integer iPStep, doublereal ITol, integer MaxIt, doublereal etaMx, doublereal T, const NonlinearSolverOptions &options)
 
 ~MatrixFreeSolver (void)
 
- Public Member Functions inherited from NonlinearSolver
 NonlinearSolver (const NonlinearSolverOptions &options)
 
virtual void SetTest (NonlinearSolverTest *pr, NonlinearSolverTest *ps)
 
virtual ~NonlinearSolver (void)
 
virtual bool MakeResTest (Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
 
virtual integer TotalAssembledJacobian (void)
 
virtual NonlinearSolverTestpGetResTest (void)
 
virtual NonlinearSolverTestpGetSolTest (void)
 
- Public Member Functions inherited from SolverDiagnostics
 SolverDiagnostics (unsigned OF=OUTPUT_DEFAULT, DriveCaller *pOM=0)
 
virtual ~SolverDiagnostics (void)
 
void SetNoOutput (void)
 
void SetOutputMeter (DriveCaller *pOM)
 
void SetOutputDriveHandler (const DriveHandler *pDH)
 
void SetOutputFlags (unsigned OF)
 
void AddOutputFlags (unsigned OF)
 
void DelOutputFlags (unsigned OF)
 
MatrixHandler::Norm_t GetCondMatNorm (void) const
 
bool outputMeter (void) const
 
bool outputIters (void) const
 
bool outputRes (void) const
 
bool outputSol (void) const
 
bool outputJac (void) const
 
bool outputStep (void) const
 
bool outputBailout (void) const
 
bool outputCounter (void) const
 
bool outputMatrixConditionNumber (void) const
 
bool outputSolverConditionNumber (void) const
 
bool outputSolverConditionStat (void) const
 
bool outputCPUTime (void) const
 
bool outputMsg (void) const
 

Private Attributes

MyVectorHandler rHat
 
MyVectorHandler p
 
MyVectorHandler pHat
 
MyVectorHandler s
 
MyVectorHandler sHat
 
MyVectorHandler t
 
MyVectorHandler v
 
MyVectorHandler dx
 

Additional Inherited Members

- Public Types inherited from MatrixFreeSolver
enum  SolverType {
  UNKNOWN = -1, BICGSTAB, GMRES, DEFAULT = GMRES,
  LAST
}
 
- Public Types inherited from NonlinearSolver
enum  Type {
  UNKNOWN = -1, NEWTONRAPHSON, MATRIXFREE, LINESEARCH,
  DEFAULT = NEWTONRAPHSON, LASTSOLVERTYPE
}
 
- Protected Types inherited from NonlinearSolver
enum  CPUTimeType { CPU_RESIDUAL, CPU_JACOBIAN, CPU_LINEAR_SOLVER, CPU_LAST_TYPE }
 
- Protected Types inherited from SolverDiagnostics
enum  {
  OUTPUT_NONE = 0x0000, OUTPUT_ITERS = 0x0001, OUTPUT_RES = 0x0002, OUTPUT_SOL = 0x0004,
  OUTPUT_JAC = 0x0008, OUTPUT_BAILOUT = 0x0010, OUTPUT_MSG = 0x0020, OUTPUT_COUNTER = 0x0040,
  OUTPUT_MAT_COND_NUM_1 = 0x0080, OUTPUT_MAT_COND_NUM_INF = 0x0100, OUTPUT_SOLVER_COND_NUM = 0x0200, OUTPUT_SOLVER_COND_STAT = 0x400,
  OUTPUT_CPU_TIME = 0x800, OUTPUT_MAT_COND_NUM = OUTPUT_MAT_COND_NUM_1 | OUTPUT_MAT_COND_NUM_INF, OUTPUT_DEFAULT = OUTPUT_MSG, OUTPUT_STEP = (OUTPUT_ITERS | OUTPUT_RES | OUTPUT_SOL | OUTPUT_JAC | OUTPUT_MAT_COND_NUM | OUTPUT_SOLVER_COND_NUM),
  OUTPUT_MASK = 0x07FF
}
 
- Protected Types inherited from NonlinearSolverOptions
enum  ScaleFlags { SCALE_ALGEBRAIC_EQUATIONS_NO = 0, SCALE_ALGEBRAIC_EQUATIONS_YES = 1 }
 
- Protected Member Functions inherited from NonlinearSolver
virtual bool MakeSolTest (Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
 
doublereal dGetCondMax () const
 
doublereal dGetCondMin () const
 
doublereal dGetCondAvg () const
 
doublereal dGetTimeCPU (CPUTimeType eType) const
 
void AddCond (doublereal dCond)
 
void AddTimeCPU (doublereal dTime, CPUTimeType eType)
 
- Protected Member Functions inherited from NonlinearSolverOptions
 NonlinearSolverOptions (bool bHonorJacRequest=false, enum ScaleFlags eScaleFlags=SCALE_ALGEBRAIC_EQUATIONS_NO, doublereal dScaleAlgebraic=1.)
 
- Protected Attributes inherited from MatrixFreeSolver
PreconditionerpPM
 
VectorHandlerpRes
 
doublereal IterTol
 
integer MaxLinIt
 
doublereal Tau
 
doublereal gamma
 
doublereal etaMax
 
integer PrecondIter
 
bool bBuildMat
 
const NonlinearProblempPrevNLP
 
- Protected Attributes inherited from NonlinearSolver
integer Size
 
integer TotJac
 
NonlinearSolverTestpResTest
 
NonlinearSolverTestpSolTest
 
- Protected Attributes inherited from SolverDiagnostics
unsigned OutputFlags
 
DriveCallerpOutputMeter
 
- Protected Attributes inherited from NonlinearSolverOptions
bool bHonorJacRequest
 
enum
NonlinearSolverOptions::ScaleFlags 
eScaleFlags
 
doublereal dScaleAlgebraic
 

Detailed Description

Definition at line 45 of file bicg.h.

Constructor & Destructor Documentation

BiCGStab::BiCGStab ( const Preconditioner::PrecondType  PType,
const integer  iPStep,
doublereal  ITol,
integer  MaxIt,
doublereal  etaMx,
doublereal  T,
const NonlinearSolverOptions options 
)

Definition at line 57 of file bicg.cc.

References NO_OP.

64 : MatrixFreeSolver(PType, iPStep, ITol, MaxIt, etaMx, T, options)
65 {
66  NO_OP;
67 }
#define NO_OP
Definition: myassert.h:74
MatrixFreeSolver(const Preconditioner::PrecondType PType, const integer iPStep, doublereal ITol, integer MaxIt, doublereal etaMx, doublereal T, const NonlinearSolverOptions &options)
Definition: mfree.cc:47
BiCGStab::~BiCGStab ( void  )

Definition at line 69 of file bicg.cc.

References NO_OP.

70 {
71  NO_OP;
72 }
#define NO_OP
Definition: myassert.h:74

Member Function Documentation

void BiCGStab::Solve ( const NonlinearProblem NLP,
Solver pS,
const integer  iMaxIter,
const doublereal Tol,
integer iIterCnt,
doublereal dErr,
const doublereal SolTol,
doublereal dSolErr 
)
virtual

Implements NonlinearSolver.

Definition at line 75 of file bicg.cc.

References ASSERT, MatrixFreeSolver::bBuildMat, NonlinearSolverOptions::bHonorJacRequest, Solver::CheckTimeStepLimit(), DEBUGCOUT, dx, MatrixFreeSolver::etaMax, NonlinearProblem::EvalProd(), grad::fabs(), MatrixFreeSolver::gamma, VectorHandler::iGetSize(), VectorHandler::InnerProd(), NonlinearProblem::Jacobian(), NonlinearSolver::MakeResTest(), NonlinearSolver::MakeSolTest(), SolutionManager::MatrInitialize(), SolutionManager::MatrReset(), MatrixFreeSolver::MaxLinIt, MBDYN_EXCEPT_ARGS, mbdyn_stop_at_end_of_iteration(), VectorHandler::Norm(), SolverDiagnostics::outputBailout(), SolverDiagnostics::outputIters(), SolverDiagnostics::outputRes(), SolverDiagnostics::outputSol(), p, Solver::pGetSolutionManager(), pHat, SolutionManager::pMatHdl(), MatrixFreeSolver::pPM, MatrixFreeSolver::pPrevNLP, Preconditioner::Precond(), MatrixFreeSolver::PrecondIter, MatrixFreeSolver::pRes, SolutionManager::pResHdl(), Solver::PrintResidual(), Solver::PrintSolution(), VectorHandler::Reset(), MyVectorHandler::Reset(), NonlinearProblem::Residual(), MyVectorHandler::Resize(), rHat, s, VectorHandler::ScalarAddMul(), MyVectorHandler::ScalarAddMul(), sHat, NonlinearSolver::Size, t, MatrixFreeSolver::Tau, NonlinearSolver::TotJac, NonlinearProblem::Update(), and v.

83 {
84  ASSERT(pNLP != NULL);
85  ASSERT(pS != NULL);
86 
88 
89  iIterCnt = 0;
90  dSolErr = 0.;
91 
92  /* external nonlinear iteration */
93 
94  /* riassembla sempre lo jacobiano se l'integratore e' nuovo */
95  if (pNLP != pPrevNLP) {
96  bBuildMat = true;
97  }
98 
99  if (!PrecondIter) {
100  bBuildMat = true;
101  }
102 
103  pPrevNLP = pNLP;
104  pRes = pSM->pResHdl();
105  Size = pRes->iGetSize();
106 
107  doublereal eta = etaMax;
108  doublereal rateo = 0.;
109  doublereal Fnorm = 1.;
110  doublereal resid;
111  doublereal rho_1;
112  doublereal rho_2 = 0.;
113  doublereal alpha;
114  doublereal beta;
115  doublereal omega;
116  VectorHandler* pr;
117 
118  /*
119  * these will be resized (actually allocated)
120  * only the first time they're used, unless
121  * the size of the problem changes
122  *
123  * FIXME: need to review this code.
124  */
125  rHat.Resize(Size);
126  p.Resize(Size);
127  pHat.Resize(Size);
128  s.Resize(Size);
129  sHat.Resize(Size);
130  t.Resize(Size);
131  v.Resize(Size);
132  dx.Resize(Size);
133 
134  integer TotalIter = 0;
135 
136 #ifdef DEBUG_ITERATIVE
137  std::cout << " BiCGStab New Step " <<std::endl;
138 #endif /* DEBUG_ITERATIVE */
139 
140  doublereal dOldErr = 1.; //initialize to silence g++ warning
141  doublereal dErrFactor = 1.;
142  doublereal dErrDiff = 0.;
143 
144  while (true) {
145 
146 #ifdef USE_EXTERNAL
147  SendExternal();
148 #endif /* USE_EXTERNAL */
149 
150  pRes->Reset();
151  try {
152  pNLP->Residual(pRes);
153  }
155  if (bHonorJacRequest) {
156  bBuildMat = true;
157  }
158  }
159 
160  if (outputRes()) {
161  pS->PrintResidual(*pRes, iIterCnt);
162  }
163 
164  bool bTest = MakeResTest(pS, pNLP, *pRes, Tol, dErr, dErrDiff);
165  if (iIterCnt > 0) {
166  dErrFactor *= dErr/dOldErr;
167  }
168  dOldErr = dErr;
169 
170 #ifdef DEBUG_ITERATIVE
171  std::cerr << "dErr " << dErr << std::endl;
172 #endif /* DEBUG_ITERATIVE */
173 
174  if (outputIters()) {
175 #ifdef USE_MPI
176  if (!bParallel || MBDynComm.Get_rank() == 0)
177 #endif /* USE_MPI */
178  {
179  silent_cout("\tIteration(" << iIterCnt << ") " << dErr);
180  if (bBuildMat && !bTest) {
181  silent_cout(" J");
182  }
183  silent_cout(std::endl);
184  }
185  }
186 
187  pS->CheckTimeStepLimit(dErr, dErrDiff);
188 
189  if (bTest) {
190  return;
191  }
192  if (!std::isfinite(dErr)) {
193  throw ErrSimulationDiverged(MBDYN_EXCEPT_ARGS);
194  }
195  if (iIterCnt >= std::abs(iMaxIter)) {
196  if (iMaxIter < 0 && dErrFactor < 1.) {
197  return;
198  }
199  if (outputBailout()) {
200  pS->PrintResidual(*pRes, iIterCnt);
201  }
202  throw NoConvergence(MBDYN_EXCEPT_ARGS);
203  }
204  rateo = dErr*dErr/Fnorm;
205  Fnorm = dErr*dErr;
206 
207  iIterCnt++;
208 
209  /* inner iteration to solve the linear system */
210 
211  /* BiCGSTAB Iterative solver */
212  DEBUGCOUT("Using BiCGStab iterative solver" << std::endl);
213 
214  /* r0 = b- A*x0 but we choose (x0 = 0) => r0 = b */
215  /* N.B. *pRes = -F(0) */
216 
217  pr = pRes;
218  doublereal LocTol = eta * dErr;
219 
220  rHat = *pr;
221 
222 #ifdef DEBUG_ITERATIVE
223  std::cerr << "LocTol " << LocTol << std::endl;
224 #endif /* DEBUG_ITERATIVE */
225 
226  rho_1 = dErr*dErr; /*rhat.InnerProd(r); */
227  resid = dErr;
228  v.Reset();
229  t.Reset();
230  p.Reset();
231  dx.Reset();
232 
233  if (bBuildMat) {
234  pSM->MatrReset();
235 
236 rebuild_matrix:;
237  try {
238  pNLP->Jacobian(pSM->pMatHdl());
239 
241  silent_cout("NewtonRaphsonSolver: "
242  "rebuilding matrix..."
243  << std::endl);
244 
245  /* need to rebuild the matrix... */
246  pSM->MatrInitialize();
247  goto rebuild_matrix;
248 
249  } catch (...) {
250  throw;
251  }
252 
253  bBuildMat = false;
254  TotalIter = 0;
255  TotJac++;
256 
257 #ifdef DEBUG_ITERATIVE
258  std::cerr << "Jacobian " << std::endl;
259 #endif /* DEBUG_ITERATIVE */
260 
261  }
262 
263 #ifdef DEBUG_ITERATIVE
264  std::cerr << "rho_1 " << rho_1 << std::endl;
265 #endif /* DEBUG_ITERATIVE */
266 
267  int It = 0;
268  while ((resid > LocTol) && (It++ < MaxLinIt)) {
269  if (It == 1) {
270  p = *pr;
271  } else {
272  rho_1 = rHat.InnerProd(*pr);
273 
274 #ifdef DEBUG_ITERATIVE
275  std::cerr << "rho_1 " << rho_1 << std::endl;
276 #endif /* DEBUG_ITERATIVE */
277 
278  if (fabs(rho_1) < std::numeric_limits<doublereal>::epsilon()) {
279  silent_cout("Bi-CGStab Iterative "
280  "Solver breakdown"
281  << " rho_1 = 0 "
282  << std::endl);
283  break;
284  }
285  beta = (rho_1/rho_2) * (alpha/omega);
286 
287 #ifdef DEBUG_ITERATIVE
288  std::cerr << "beta " << beta << std::endl;
289 #endif /* DEBUG_ITERATIVE */
290 
291  p.ScalarAddMul(*pr, p.ScalarAddMul(v, -omega), beta);
292  }
293  /* right preconditioning */
294  pPM->Precond(p, pHat, pSM);
295  pNLP->EvalProd(Tau, rHat, pHat, v);
296 #if 0
297  (pSM->pMatHdl())->MatVecMul(v,pHat);
298 #endif
299 #ifdef DEBUG_ITERATIVE
300  std::cout << "v:" << std::endl;
301  for (int iTmpCnt = 1; iTmpCnt <= Size; iTmpCnt++) {
302  std::cout << "Dof" << std::setw(4) << iTmpCnt << ": "
303  << v(iTmpCnt) << std::endl;
304  }
305 #endif /* DEBUG_ITERATIVE */
306 
307  alpha = rHat.InnerProd(v);
308  alpha = rho_1 / alpha;
309 
310 #ifdef DEBUG_ITERATIVE
311  std::cerr << "alpha " << alpha << std::endl;
312 #endif /* DEBUG_ITERATIVE */
313 
314  s.ScalarAddMul(*pr, v, -alpha);
315 
316 #ifdef DEBUG_ITERATIVE
317  std::cerr << "s.Norm() " << s.Norm() << std::endl;
318 #endif /* DEBUG_ITERATIVE */
319 
320  if ((resid = s.Norm()) < LocTol) {
321  dx.ScalarAddMul(pHat, alpha);
322  TotalIter++;
323  break;
324  }
325  pPM->Precond(s, sHat, pSM);
326  pNLP->EvalProd(Tau, rHat, sHat, t);
327 #if 0
328  (pSM->pMatHdl())->MatVecMul(t,sHat);
329 #endif
330  omega = t.Norm();
331  omega = t.InnerProd(s) / omega;
332 
333 #ifdef DEBUG_ITERATIVE
334  std::cerr << "omega " << omega << std::endl;
335 #endif /* DEBUG_ITERATIVE */
336 
337  dx.ScalarAddMul(pHat, alpha);
338  dx.ScalarAddMul(sHat, omega);
339  pr->ScalarAddMul(s, t, -omega);
340  rho_2 = rho_1;
341  resid = pr->Norm();
342 
343 #ifdef DEBUG_ITERATIVE
344  std::cerr << "resid " << resid << std::endl;
345 #endif /* DEBUG_ITERATIVE */
346 
347  TotalIter++;
348  if (fabs(omega) < std::numeric_limits<doublereal>::epsilon()) {
349  silent_cout("Bi-CGStab Iterative Solver "
350  "breakdown omega = 0 " << std::endl);
351  break;
352  }
353  if (It == MaxLinIt) {
354  silent_cerr("Iterative inner solver didn't "
355  "converge. Continuing..."
356  << std::endl);
357  }
358  }
359  /* se ha impiegato troppi passi riassembla lo jacobiano */
360 
361  if (TotalIter >= PrecondIter && PrecondIter) {
362  bBuildMat = true;
363  }
364  /* calcola il nuovo eta */
365 
366  doublereal etaNew = gamma * rateo;
367  doublereal etaBis;
368  if ((etaBis = gamma*eta*eta) > .1) {
369  etaNew = (etaNew > etaBis) ? etaNew : etaBis;
370  }
371  eta = (etaNew < etaMax) ? etaNew : etaMax;
372  /* prevent oversolving */
373  etaBis = .5*Tol/dErr;
374  eta = (eta > etaBis) ? eta : etaBis;
375 
376 #ifdef DEBUG_ITERATIVE
377  std::cerr << "eta " << eta << std::endl;
378 #endif /* DEBUG_ITERATIVE */
379 
380  if (outputSol()) {
381  pS->PrintSolution(dx, iIterCnt);
382  }
383 
384  pNLP->Update(&dx);
385 
386  bTest = MakeSolTest(pS, dx, SolTol, dSolErr);
387  if (outputIters()) {
388 #ifdef USE_MPI
389  if (!bParallel || MBDynComm.Get_rank() == 0)
390 #endif /* USE_MPI */
391  {
392  silent_cout("\t\tSolErr " << dSolErr
393  << std::endl);
394  }
395  }
396 
397  if (bTest) {
398  throw ConvergenceOnSolution(MBDYN_EXCEPT_ARGS);
399  }
400 
401  // allow to bail out in case of multiple CTRL^C
404  }
405  }
406 }
Preconditioner * pPM
Definition: mfree.h:61
virtual void Reset(void)=0
virtual VectorHandler * pResHdl(void) const =0
virtual SolutionManager * pGetSolutionManager(void) const
Definition: solver.h:398
bool outputIters(void) const
integer TotJac
Definition: nonlin.h:252
#define MBDYN_EXCEPT_ARGS
Definition: except.h:63
bool outputBailout(void) const
MyVectorHandler pHat
Definition: bicg.h:48
virtual bool MakeSolTest(Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
Definition: nonlin.cc:497
virtual void PrintResidual(const VectorHandler &Res, integer iIterCnt) const
Definition: solver.cc:5194
virtual bool MakeResTest(Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
Definition: nonlin.cc:474
MyVectorHandler s
Definition: bicg.h:48
bool bBuildMat
Definition: mfree.h:69
virtual void CheckTimeStepLimit(doublereal dErr, doublereal dErrDiff) const
Definition: solver.cc:5205
virtual doublereal InnerProd(const VectorHandler &VH) const
Definition: vh.cc:269
virtual VectorHandler & ScalarAddMul(const VectorHandler &VH, const VectorHandler &VH1, const doublereal &d)
Definition: vh.cc:499
bool outputRes(void) const
integer Size
Definition: nonlin.h:251
MyVectorHandler rHat
Definition: bicg.h:48
MyVectorHandler dx
Definition: bicg.h:48
virtual integer iGetSize(void) const =0
virtual doublereal Norm(void) const
Definition: vh.cc:262
MyVectorHandler p
Definition: bicg.h:48
GradientExpression< UnaryExpr< FuncFabs, Expr > > fabs(const GradientExpression< Expr > &u)
Definition: gradient.h:2973
virtual void Resize(integer iSize)
Definition: vh.cc:347
MyVectorHandler t
Definition: bicg.h:48
int mbdyn_stop_at_end_of_iteration(void)
Definition: solver.cc:172
integer PrecondIter
Definition: mfree.h:68
bool outputSol(void) const
const NonlinearProblem * pPrevNLP
Definition: mfree.h:70
virtual MatrixHandler * pMatHdl(void) const =0
MyVectorHandler sHat
Definition: bicg.h:48
#define DEBUGCOUT(msg)
Definition: myassert.h:232
virtual void MatrReset(void)=0
MyVectorHandler v
Definition: bicg.h:48
#define ASSERT(expression)
Definition: colamd.c:977
virtual void Precond(VectorHandler &b, VectorHandler &x, SolutionManager *pSM) const =0
virtual void MatrInitialize(void)
Definition: solman.cc:72
virtual void Reset(void)
Definition: vh.cc:459
integer MaxLinIt
Definition: mfree.h:64
doublereal Tau
Definition: mfree.h:65
virtual VectorHandler & ScalarAddMul(const VectorHandler &VH, const doublereal &d)
Definition: vh.cc:108
VectorHandler * pRes
Definition: mfree.h:62
virtual void PrintSolution(const VectorHandler &Sol, integer iIterCnt) const
Definition: solver.cc:5200
doublereal etaMax
Definition: mfree.h:67
double doublereal
Definition: colamd.c:52
long int integer
Definition: colamd.c:51
doublereal gamma
Definition: mfree.h:66

Here is the call graph for this function:

Member Data Documentation

MyVectorHandler BiCGStab::dx
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::p
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::pHat
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::rHat
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::s
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::sHat
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::t
private

Definition at line 48 of file bicg.h.

Referenced by Solve().

MyVectorHandler BiCGStab::v
private

Definition at line 48 of file bicg.h.

Referenced by Solve().


The documentation for this class was generated from the following files: