MBDyn-1.7.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups
solver.cc
Go to the documentation of this file.
1 /* $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/base/solver.cc,v 1.301 2017/10/14 23:58:06 masarati Exp $ */
2 /*
3  * MBDyn (C) is a multibody analysis code.
4  * http://www.mbdyn.org
5  *
6  * Copyright (C) 1996-2017
7  *
8  * Pierangelo Masarati <masarati@aero.polimi.it>
9  * Paolo Mantegazza <mantegazza@aero.polimi.it>
10  *
11  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
12  * via La Masa, 34 - 20156 Milano, Italy
13  * http://www.aero.polimi.it
14  *
15  * Changing this copyright notice is forbidden.
16  *
17  * This program is free software; you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation (version 2 of the License).
20  *
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program; if not, write to the Free Software
29  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
30  */
31 
32 /*
33  * Copyright 1999-2017 Giuseppe Quaranta <quaranta@aero.polimi.it>
34  * Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano
35  *
36  * This copyright statement applies to the MPI related code, which was
37  * merged from files schur.h/schur.cc
38  */
39 
40 /*
41  *
42  * Copyright (C) 2003-2017
43  * Giuseppe Quaranta <quaranta@aero.polimi.it>
44  *
45  */
46 
47 /* metodo per la soluzione del modello */
48 
49 #include "mbconfig.h" /* This goes first in every *.c,*.cc file */
50 
51 /* required for configure time macros with paths */
52 #include "mbdefs.h"
53 
54 #include <cstdlib>
55 #include <cstring>
56 #include <limits>
57 #include <ac/unistd.h>
58 #include <cerrno>
59 #include <csignal>
60 #include <cfloat>
61 #include <cmath>
62 #include <vector>
63 #include <algorithm>
64 #include "ac/sys_sysinfo.h"
65 
66 #include "solver.h"
67 #include "dataman.h"
68 #include "mtdataman.h"
69 #include "thirdorderstepsol.h"
70 #include "nr.h"
71 #include "linesearch.h"
72 #include "bicg.h"
73 #include "gmres.h"
74 #include "solman.h"
75 #include "readlinsol.h"
76 #include "ls.h"
77 #include "naivewrap.h"
78 #include "Rot.hh"
79 #include "cleanup.h"
80 #include "drive_.h"
81 #include "TimeStepControl.h"
82 #include "solver_impl.h"
83 
84 #include "ac/lapack.h"
85 #include "ac/arpack.h"
86 #include "eigjdqz.h"
87 
88 
89 #ifdef HAVE_SIGNAL
90 /*
91  * MBDyn starts with mbdyn_keep_going == MBDYN_KEEP_GOING.
92  *
93  * A single CTRL^C sets it to MBDYN_STOP_AT_END_OF_TIME_STEP,
94  * which results in exiting at the end of the time step,
95  * after the output in case of success.
96  *
97  * A second CTRL^C sets it to MBDYN_STOP_AT_END_OF_ITERATION,
98  * which results in exiting at the end of the current iteration,
99  * after printing debug output if required.
100  *
101  * A further CTRL^C sets it to MBDYN_STOP_NOW and in throwing
102  * an exception.
103  */
104 enum {
105  MBDYN_KEEP_GOING = 0,
106  MBDYN_STOP_AT_END_OF_TIME_STEP = 1,
107  MBDYN_STOP_AT_END_OF_ITERATION = 2,
108  MBDYN_STOP_NOW = 3
109 
110 };
111 
112 volatile sig_atomic_t mbdyn_keep_going = MBDYN_KEEP_GOING;
113 #if defined(__FreeBSD__)
114 __sighandler_t *mbdyn_sh_term = SIG_DFL;
115 __sighandler_t *mbdyn_sh_int = SIG_DFL;
116 __sighandler_t *mbdyn_sh_hup = SIG_DFL;
117 __sighandler_t *mbdyn_sh_pipe = SIG_DFL;
118 #else // !defined(__FreeBSD__)
119 __sighandler_t mbdyn_sh_term = SIG_DFL;
120 __sighandler_t mbdyn_sh_int = SIG_DFL;
121 __sighandler_t mbdyn_sh_hup = SIG_DFL;
122 __sighandler_t mbdyn_sh_pipe = SIG_DFL;
123 #endif // !defined(__FreeBSD__)
124 
125 extern "C" void
126 mbdyn_really_exit_handler(int signum)
127 {
128  ::mbdyn_keep_going = MBDYN_STOP_NOW;
129  switch (signum) {
130  case SIGTERM:
131  signal(signum, ::mbdyn_sh_term);
132  break;
133 
134  case SIGINT:
135  signal(signum, ::mbdyn_sh_int);
136  break;
137 
138 #ifdef SIGHUP
139  case SIGHUP:
140  signal(signum, ::mbdyn_sh_hup);
141  break;
142 #endif // SIGHUP
143 
144 #ifdef SIGPIPE
145  case SIGPIPE:
146  signal(signum, ::mbdyn_sh_pipe);
147  break;
148 #endif // SIGPIPE
149  }
150 
151  mbdyn_cleanup();
152 
154 }
155 
156 extern "C" void
157 mbdyn_modify_last_iteration_handler(int signum)
158 {
159  ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_ITERATION;
160  signal(signum, mbdyn_really_exit_handler);
161 }
162 
163 extern "C" void
164 mbdyn_modify_final_time_handler(int signum)
165 {
166  ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_TIME_STEP;
167  signal(signum, mbdyn_modify_last_iteration_handler);
168 }
169 #endif /* HAVE_SIGNAL */
170 
171 extern "C" int
173 {
174 #ifdef HAVE_SIGNAL
175  return ::mbdyn_keep_going >= MBDYN_STOP_AT_END_OF_ITERATION;
176 #else // ! HAVE_SIGNAL
177  return 0;
178 #endif // ! HAVE_SIGNAL
179 }
180 
181 extern "C" int
183 {
184 #ifdef HAVE_SIGNAL
185  return ::mbdyn_keep_going >= MBDYN_STOP_AT_END_OF_TIME_STEP;
186 #else // ! HAVE_SIGNAL
187  return 0;
188 #endif // ! HAVE_SIGNAL
189 }
190 
191 extern "C" void
193 {
194 #ifdef HAVE_SIGNAL
195  ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_ITERATION;
196 #endif // HAVE_SIGNAL
197 }
198 
199 extern "C" void
201 {
202 #ifdef HAVE_SIGNAL
203  ::mbdyn_keep_going = MBDYN_STOP_AT_END_OF_TIME_STEP;
204 #endif // HAVE_SIGNAL
205 }
206 
207 extern "C" void
209 {
210 #ifdef HAVE_SIGNAL
211 #if defined(__FreeBSD__)
212  __sighandler_t *hdl;
213 #else // ! defined(__FreeBSD__)
214  __sighandler_t hdl;
215 #endif // ! defined(__FreeBSD__)
216 
217  if (pre) {
218  hdl = mbdyn_really_exit_handler;
219 
220  } else {
221  hdl = mbdyn_modify_final_time_handler;
222  }
223  /*
224  * FIXME: don't do this if compiling with USE_RTAI
225  * Re FIXME: use sigaction() ...
226  */
227  ::mbdyn_sh_term = signal(SIGTERM, hdl);
228  if (::mbdyn_sh_term == SIG_IGN) {
229  signal(SIGTERM, SIG_IGN);
230  }
231 
232  ::mbdyn_sh_int = signal(SIGINT, hdl);
233  if (::mbdyn_sh_int == SIG_IGN) {
234  signal(SIGINT, SIG_IGN);
235  }
236 
237 #ifdef SIGHUP
238  ::mbdyn_sh_hup = signal(SIGHUP, hdl);
239  if (::mbdyn_sh_hup == SIG_IGN) {
240  signal(SIGHUP, SIG_IGN);
241  }
242 #endif // SIGHUP
243 
244 #ifdef SIGPIPE
245  ::mbdyn_sh_pipe = signal(SIGPIPE, hdl);
246  if (::mbdyn_sh_pipe == SIG_IGN) {
247  signal(SIGPIPE, SIG_IGN);
248  }
249 #endif // SIGPIPE
250 #endif /* HAVE_SIGNAL */
251 }
252 
253 int
254 mbdyn_reserve_stack(unsigned long size)
255 {
256  int buf[size];
257 
258 #ifdef HAVE_MEMSET
259  memset(buf, 0, size*sizeof(int));
260 #else /* !HAVE_MEMSET */
261  for (unsigned long i = 0; i < size; i++) {
262  buf[i] = 0;
263  }
264 #endif /* !HAVE_MEMSET */
265 
266 #ifdef HAVE_MLOCKALL
267  return mlockall(MCL_CURRENT | MCL_FUTURE);
268 #else /* !HAVE_MLOCKALL */
269  return 0;
270 #endif /* !HAVE_MLOCKALL */
271 }
272 
273 /* Costruttore: esegue la simulazione */
275  const std::string& sInFName,
276  const std::string& sOutFName,
277  unsigned int nThreads,
278  bool bPar)
279 :
280 #ifdef USE_MULTITHREAD
281 nThreads(nThreads),
282 #endif /* USE_MULTITHREAD */
283 pTSC(0),
284 dCurrTimeStep(0.),
285 iStIter(0),
286 dTime(0.),
287 MaxTimeStep(new ConstDriveCaller(::dDefaultMaxTimeStep)),
288 dMinTimeStep(::dDefaultMinTimeStep),
289 CurrStep(StepIntegrator::NEWSTEP),
290 sInputFileName(sInFName),
291 sOutputFileName(sOutFName),
292 HP(HPar),
293 iMaxIterations(::iDefaultMaxIterations),
294 EigAn(),
295 pRTSolver(0),
296 iNumPreviousVectors(2),
297 iUnkStates(1),
298 pdWorkSpace(0),
299 qX(),
300 qXPrime(),
301 pX(0),
302 pXPrime(0),
303 dInitialTime(0.),
304 dFinalTime(0.),
305 dRefTimeStep(0.),
306 dInitialTimeStep(1.),
307 dMaxResidual(std::numeric_limits<doublereal>::max()),
308 dMaxResidualDiff(std::numeric_limits<doublereal>::max()),
309 eTimeStepLimit(TS_SOFT_LIMIT),
310 iDummyStepsNumber(::iDefaultDummyStepsNumber),
311 dDummyStepsRatio(::dDefaultDummyStepsRatio),
312 eAbortAfter(AFTER_UNKNOWN),
313 RegularType(INT_UNKNOWN),
314 DummyType(INT_UNKNOWN),
315 pDerivativeSteps(0),
316 pFirstDummyStep(0),
317 pDummySteps(0),
318 pFirstRegularStep(0),
319 pRegularSteps(0),
320 pCurrStepIntegrator(0),
321 pRhoRegular(0),
322 pRhoAlgebraicRegular(0),
323 pRhoDummy(0),
324 pRhoAlgebraicDummy(0),
325 dDerivativesCoef(::dDefaultDerivativesCoefficient),
326 CurrLinearSolver(),
327 ResTest(NonlinearSolverTest::NORM),
328 SolTest(NonlinearSolverTest::NONE),
329 bScale(false),
330 bTrueNewtonRaphson(true),
331 bKeepJac(false),
332 iIterationsBeforeAssembly(0),
333 NonlinearSolverType(NonlinearSolver::UNKNOWN),
334 /* for matrix-free solvers */
335 MFSolverType(MatrixFreeSolver::UNKNOWN),
336 dIterTol(::dDefaultTol),
337 PcType(Preconditioner::FULLJACOBIANMATRIX),
338 iPrecondSteps(::iDefaultPreconditionerSteps),
339 iIterativeMaxSteps(::iDefaultPreconditionerSteps),
340 dIterertiveEtaMax(defaultIterativeEtaMax),
341 dIterertiveTau(defaultIterativeTau),
342 /* end of matrix-free solvers */
343 /* for line search solver */
344 LineSearch(),
345 /* end of line search solver */
346 /* for parallel solvers */
347 bParallel(bPar),
348 pSDM(0),
349 iNumLocDofs(0),
350 iNumIntDofs(0),
351 pLocDofs(0),
352 pIntDofs(0),
353 pDofs(0),
354 pLocalSM(0),
355 /* end of parallel solvers */
356 pDM(0),
357 iNumDofs(0),
358 pSM(0),
359 pNLS(0),
360 eStatus(SOLVER_STATUS_UNINITIALIZED),
361 bOutputCounter(false),
362 outputCounterPrefix(),
363 outputCounterPostfix(),
364 iTotIter(0),
365 dTotErr(0.),
366 dTest(std::numeric_limits<double>::max()),
367 dSolTest(std::numeric_limits<double>::max()),
368 bSolConv(false),
369 bOut(false),
370 lStep(0)
371 {
372  DEBUGCOUTFNAME("Solver::Solver");
374  ASSERT(!sInFName.empty());
375 }
376 
377 #ifdef USE_MULTITHREAD
378 void
379 Solver::ThreadPrepare(void)
380 {
381  /* check for thread potential */
382  if (nThreads == 0) {
383  int n = get_nprocs();
384 
385  if (n > 1) {
386  silent_cout("no multithread requested "
387  "with a potential of " << n
388  << " CPUs" << std::endl);
389  nThreads = n;
390 
391  } else {
392  nThreads = 1;
393  }
394  }
395 }
396 #endif /* USE_MULTITHREAD */
397 
398 bool
400 {
401  DEBUGCOUTFNAME("Solver::Prepare");
402 
403  // consistency check
405  silent_cerr("Prepare() must be called first" << std::endl);
407  }
408 
410 
411  /* Legge i dati relativi al metodo di integrazione */
412  ReadData(HP);
413 
414 #ifdef USE_MULTITHREAD
415  ThreadPrepare();
416 #endif /* USE_MULTITHREAD */
417 
418  if (pRTSolver) {
419  pRTSolver->Setup();
420  }
421 
422 #ifdef USE_MPI
423  int mpi_finalize = 0;
424 #endif
425 
426 #ifdef USE_SCHUR
427  if (bParallel) {
428  DEBUGLCOUT(MYDEBUG_MEM, "creating parallel SchurDataManager"
429  << std::endl);
430 
434  OutputFlags,
435  this,
436  dInitialTime,
437  sOutputFileName.c_str(),
438  sInputFileName.c_str(),
440 
441  pDM = pSDM;
442 
443  } else
444 #endif // USE_SCHUR
445  {
446  /* chiama il gestore dei dati generali della simulazione */
447 #ifdef USE_MULTITHREAD
448  if (nThreads > 1) {
450  /* conservative: dir may use too much memory */
452  bool b;
453 
454 #if defined(USE_UMFPACK)
457 #elif defined(USE_Y12)
460 #else
461  b = false;
462 #endif
463  if (!b) {
464  silent_cerr("unable to select a CC-capable solver"
465  << std::endl);
467  }
468  }
469  }
470 
471  silent_cout("Creating multithread solver "
472  "with " << nThreads << " threads "
473  "and "
475  << " linear solver"
476  << std::endl);
477 
479  MultiThreadDataManager,
480  MultiThreadDataManager(HP,
481  OutputFlags,
482  this,
483  dInitialTime,
484  sOutputFileName.c_str(),
485  sInputFileName.c_str(),
487  nThreads));
488 
489  } else
490 #endif /* USE_MULTITHREAD */
491  {
492  DEBUGLCOUT(MYDEBUG_MEM, "creating DataManager"
493  << std::endl);
494 
495  silent_cout("Creating scalar solver "
496  "with "
498  << " linear solver"
499  << std::endl);
500 
502  DataManager,
503  DataManager(HP,
504  OutputFlags,
505  this,
506  dInitialTime,
507  sOutputFileName.c_str(),
508  sInputFileName.c_str(),
510  }
511  }
512 
513  // log symbol table
514  std::ostream& log = pDM->GetLogFile();
515  log << "Symbol table:" << std::endl;
516  log << HP.GetMathParser().GetSymbolTable();
517 
518 #ifdef HAVE_ENVIRON
519  // log environment
520  log << "Environment:" << std::endl;
521  for (int i = 0; environ[i] != NULL; i++) {
522  log << " " << environ[i] << std::endl;
523  }
524 #endif // HAVE_ENVIRON
525 
526  // close input stream
527  HP.Close();
528 
529  /* Si fa dare il DriveHandler e linka i drivers di rho ecc. */
530  const DriveHandler* pDH = pDM->pGetDrvHdl();
532 
533  bOutputCounter = outputCounter() && isatty(fileno(stderr));
534  outputCounterPrefix = bOutputCounter ? "\n" : "";
535  outputCounterPostfix = outputStep() ? "\n" : "\r";
536 
537  /* Si fa dare l'std::ostream al file di output per il log */
538  std::ostream& Out = pDM->GetOutFile();
539 
540  if (eAbortAfter == AFTER_INPUT) {
541  /* Esce */
542  pDM->Output(0, dTime, 0., true);
543  Out << "End of Input; no simulation or assembly is required."
544  << std::endl;
545  return false;
546 
547  } else if (eAbortAfter == AFTER_ASSEMBLY) {
548  /* Fa l'output dell'assemblaggio iniziale e poi esce */
549  pDM->Output(0, dTime, 0., true);
550  Out << "End of Initial Assembly; no simulation is required."
551  << std::endl;
552  return false;
553  }
554 
555 #ifdef USE_SCHUR
556  /* Qui crea le partizioni: principale fra i processi, se parallelo */
557  if (bParallel) {
559  }
560 #endif // USE_SCHUR
561 
563  if (iDummyStepsNumber) {
565  }
566 
567  /* Costruisce i vettori della soluzione ai vari passi */
568  DEBUGLCOUT(MYDEBUG_MEM, "creating solution vectors" << std::endl);
569 
570 #ifdef USE_SCHUR
571  if (bParallel) {
573  pDofs = pSDM->pGetDofsList();
574 
577 
580 
581  } else
582 #endif // USE_SCHUR
583  {
584  iNumDofs = pDM->iGetNumDofs();
585  }
586 
587  /* relink those known drive callers that might need
588  * the data manager, but were verated ahead of it */
590 
591  ASSERT(iNumDofs > 0);
592 
594  integer iFSteps = 0;
595  if (iDummyStepsNumber) {
597  }
598  iNumPreviousVectors = (iRSteps < iFSteps) ? iFSteps : iRSteps;
599 
601  integer iFUnkStates = 0;
602  if (iDummyStepsNumber) {
604  }
605  iUnkStates = (iRUnkStates < iFUnkStates) ? iFUnkStates : iRUnkStates;
606 
607  /* allocate workspace for previous time steps */
610  /* allocate MyVectorHandlers for previous time steps: use workspace */
611  for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
614  MyVectorHandler(iNumDofs, pdWorkSpace+ivec*iNumDofs));
615  qX.push_back(pX);
618  MyVectorHandler(iNumDofs,
619  pdWorkSpace+((iNumPreviousVectors)+ivec)*iNumDofs));
620  qXPrime.push_back(pXPrime);
621  pX = 0;
622  pXPrime = 0;
623  }
624  /* allocate MyVectorHandlers for unknown time step(s): own memory */
627  MyVectorHandler(iUnkStates*iNumDofs));
630  MyVectorHandler(iUnkStates*iNumDofs));
631 
632 
633  /* Resetta i vettori */
634  for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
635  qX[ivec]->Reset();
636  qXPrime[ivec]->Reset();
637  }
638  pX->Reset();
639  pXPrime->Reset();
640 
641  /*
642  * Immediately link DataManager to current solution
643  *
644  * this should work as long as the last unknown time step is put
645  * at the beginning of pX, pXPrime
646  */
648 
649  /* a questo punto si costruisce il nonlinear solver */
651 
652  Scale.ResizeReset(iNumDofs);
653  if (bScale) {
654  /* collects scale factors from data manager */
655  pDM->SetScale(Scale);
656  }
657 
658  /*
659  * prepare tests for nonlinear solver;
660  *
661  * test on residual may allow pre-scaling;
662  * test on solution (difference between two iterations) does not
663  */
664  NonlinearSolverTest *pResTest = 0;
665  if (bScale) {
666  NonlinearSolverTestScale *pResTestScale = 0;
667 
668  switch (ResTest) {
670  SAFENEW(pResTestScale, NonlinearSolverTestScaleNorm);
671  break;
672 
674  SAFENEW(pResTestScale, NonlinearSolverTestScaleMinMax);
675  break;
676 
677  default:
678  ASSERT(0);
680  }
681 
682  /* registers scale factors at nonlinear solver */
683  pResTestScale->SetScale(&Scale);
684 
685  pResTest = pResTestScale;
686 
687 
688  } else {
689  switch (ResTest) {
691  SAFENEW(pResTest, NonlinearSolverTestNone);
692  break;
693 
695  SAFENEW(pResTest, NonlinearSolverTestNorm);
696  break;
697 
700  break;
701 
702  default:
703  ASSERT(0);
705  }
706  }
707 
708  NonlinearSolverTest *pSolTest = 0;
709  switch (SolTest) {
711  SAFENEW(pSolTest, NonlinearSolverTestNone);
712  break;
713 
715  SAFENEW(pSolTest, NonlinearSolverTestNorm);
716  break;
717 
720  break;
721 
722  default:
723  ASSERT(0);
725  }
726 
727  /* registers tests in nonlinear solver */
728  pNLS->SetTest(pResTest, pSolTest);
729 
730  /*
731  * Dell'assemblaggio iniziale dei vincoli se ne occupa il DataManager
732  * in quanto e' lui il responsabile dei dati della simulazione,
733  * e quindi anche della loro coerenza. Inoltre e' lui a sapere
734  * quali equazioni sono di vincolo o meno.
735  */
736 
737  pDM->SetValue(*pX, *pXPrime);
738 
739 
740  /*
741  * Prepare output
742  */
743  pDM->OutputPrepare();
744 
745  /*
746  * If eigenanalysis is requested, prepare output for it
747  */
748  if (EigAn.bAnalysis) {
751  EigAn.iFNameWidth = int(std::log10(double(EigAn.Analyses.size()))) + 1;
752  }
753  }
754 
755  /*
756  * Dialoga con il DataManager per dargli il tempo iniziale
757  * e per farsi inizializzare i vettori di soluzione e derivata
758  */
759  /* FIXME: the time is already set by DataManager, but FileDrivers
760  * have not been ServePending'd
761  */
764 
765  EigAn.currAnalysis = std::find_if(EigAn.Analyses.begin(), EigAn.Analyses.end(),
766  bind2nd(std::greater<doublereal>(), dTime));
767  if (EigAn.currAnalysis != EigAn.Analyses.end() && EigAn.currAnalysis != EigAn.Analyses.begin()) {
769  }
770 
771  // if eigenanalysis is requested and currAnalysis points
772  // past the end of the array, the analysis was requested
773  // at Time < initial time; perform *before* derivatives
774  if (EigAn.bAnalysis
775  && ((EigAn.currAnalysis == EigAn.Analyses.end()
776  && EigAn.Analyses.back() < dTime)
777  || (EigAn.currAnalysis != EigAn.Analyses.end()
778  && *EigAn.currAnalysis < dTime)))
779  {
780  Eig();
781  if (EigAn.currAnalysis != EigAn.Analyses.end()) {
783  }
784  }
785 
786  /* calcolo delle derivate */
787  DEBUGLCOUT(MYDEBUG_DERIVATIVES, "derivatives solution step"
788  << std::endl);
789 
791 
792  /* settaggio degli output Types */
793  unsigned OF = OutputFlags;
795  OF |= OUTPUT_RES;
796  }
798  OF |= OUTPUT_JAC;
799  }
801  OF |= OUTPUT_SOL;
802  }
803  pNLS->SetOutputFlags(OF);
804  if (pOutputMeter) {
807  }
808 
811  if (iDummyStepsNumber) {
816  }
821 
822 #ifdef USE_EXTERNAL
823  pNLS->SetExternal(External::EMPTY);
824 #endif /* USE_EXTERNAL */
825  /* Setup SolutionManager(s) */
827 
828  /* Derivative steps */
830  try {
831  if (outputStep()) {
832  if (outputCounter()) {
833  silent_cout(std::endl);
834  }
835  silent_cout("Derivatives t=" << dTime << " coef=" << dDerivativesCoef << std::endl);
836  }
839  qX, qXPrime, pX, pXPrime,
841  }
843  silent_cerr("Initial derivatives calculation " << iStIter
844  << " does not converge; aborting..." << std::endl
845  << "(hint: try playing with the \"derivatives coefficient\" value)" << std::endl);
846  pDM->Output(0, dTime, 0., true);
848  }
850  /*
851  * Mettere qui eventuali azioni speciali
852  * da intraprendere in caso di errore ...
853  */
855  }
856  catch (LinearSolver::ErrFactor& err) {
857  /*
858  * Mettere qui eventuali azioni speciali
859  * da intraprendere in caso di errore ...
860  */
861  silent_cerr("Initial derivatives failed because no pivot element "
862  "could be found for column " << err.iCol
863  << " (" << pDM->GetDofDescription(err.iCol) << "); "
864  "aborting..." << std::endl);
866  }
868  bSolConv = true;
869  }
870  catch (EndOfSimulation& eos) {
871  silent_cerr("Simulation ended during the derivatives steps:\n" << eos.what() << "\n");
872  return false;
873  }
874 
876  pDerivativeSteps = 0;
877 
878 #if 0
879  /* don't sum up the derivatives error */
880  dTotErr += dTest;
881 #endif
882  iTotIter += iStIter;
883 
884  if (outputMsg()) {
885  Out << "# Derivatives solution step at time " << dInitialTime
886  << " performed in " << iStIter
887  << " iterations with " << dTest
888  << " error" << std::endl;
889  }
890 
891  DEBUGCOUT("Derivatives solution step has been performed successfully"
892  " in " << iStIter << " iterations" << std::endl);
893 
894 #ifdef USE_EXTERNAL
895  /* comunica che gli ultimi dati inviati sono la condizione iniziale */
896  if (iDummyStepsNumber == 0) {
897  External::SendInitial();
898  }
899 #endif /* USE_EXTERNAL */
900 
902  /*
903  * Fa l'output della soluzione delle derivate iniziali ed esce
904  */
905 
906 
907  pDM->Output(0, dTime, 0., true);
908  Out << "End of derivatives; no simulation is required."
909  << std::endl;
910  return false;
911 
912  } else if (mbdyn_stop_at_end_of_time_step()) {
913  /*
914  * Fa l'output della soluzione delle derivate iniziali ed esce
915  */
916  pDM->Output(0, dTime, 0., true);
917  Out << "Interrupted during derivatives computation." << std::endl;
919  }
920 
921  // if eigenanalysis is requested and currAnalysis points
922  // past the end of the array, the analysis was requested
923  // at Time == initial time; perform *after* derivatives
924  if (EigAn.bAnalysis) {
925  ASSERT(EigAn.Analyses.size() > 0);
926 
927  if ((EigAn.currAnalysis == EigAn.Analyses.end()
928  && EigAn.Analyses.back() == dTime)
929  || (EigAn.currAnalysis != EigAn.Analyses.end()
930  && *EigAn.currAnalysis == dTime))
931  {
932  Eig();
933  if (EigAn.currAnalysis != EigAn.Analyses.end()) {
935  }
936  }
937  }
938 
939  /* Dati comuni a passi fittizi e normali */
940  lStep = 1;
941 
942  if (iDummyStepsNumber > 0) {
943  /* passi fittizi */
944 
945  /*
946  * inizio integrazione: primo passo a predizione lineare
947  * con sottopassi di correzione delle accelerazioni
948  * e delle reazioni vincolari
949  */
950  pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
951  Flip();
952 
955  /* FIXME: do we need to serve pending drives in dummy steps? */
957 
958  DEBUGLCOUT(MYDEBUG_FSTEPS, "Current time step: "
959  << dCurrTimeStep << std::endl);
960 
961  if (outputStep()) {
962  silent_cout("Dummy Step(" << lStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
963  }
964 
965  ASSERT(pFirstDummyStep != 0);
966 
967  /* Setup SolutionManager(s) */
969  /* pFirstDummyStep */
970  pCurrStepIntegrator = pFirstDummyStep;
971  try {
972  dTest = pFirstDummyStep->Advance(this,
973  dRefTimeStep, 1.,
975  qX, qXPrime, pX, pXPrime,
977  }
979  silent_cerr("First dummy step does not converge; "
980  "TimeStep=" << dCurrTimeStep
981  << " cannot be reduced further; "
982  "aborting..." << std::endl);
983  pDM->Output(0, dTime, dCurrTimeStep, true);
985  }
987  /*
988  * Mettere qui eventuali azioni speciali
989  * da intraprendere in caso di errore ...
990  */
992  }
993  catch (LinearSolver::ErrFactor& err) {
994  /*
995  * Mettere qui eventuali azioni speciali
996  * da intraprendere in caso di errore ...
997  */
998  silent_cerr("First dummy step failed because no pivot element "
999  "could be found for column " << err.iCol
1000  << " (" << pDM->GetDofDescription(err.iCol) << "); "
1001  "aborting..." << std::endl);
1003  }
1005  bSolConv = true;
1006  }
1007  catch (EndOfSimulation& eos) {
1008  silent_cerr("Simulation ended during the first dummy step:\n"
1009  << eos.what() << "\n");
1010  return false;
1011  }
1012 
1014  pFirstDummyStep = 0;
1015 
1017  dTime += dRefTimeStep;
1018 
1019 #if 0
1020  /* don't sum up the derivatives error */
1021  dTotErr += dTest;
1022 #endif
1023  iTotIter += iStIter;
1024 
1026  /*
1027  * Fa l'output della soluzione delle derivate iniziali
1028  * ed esce
1029  */
1030 #ifdef DEBUG_FICTITIOUS
1031  pDM->Output(0, dTime, dCurrTimeStep, true);
1032 #endif /* DEBUG_FICTITIOUS */
1033  Out << "Interrupted during first dummy step." << std::endl;
1035  }
1036 
1037 #ifdef DEBUG_FICTITIOUS
1038  pDM->Output(0, dTime, dCurrTimeStep, true);
1039 #endif /* DEBUG_FICTITIOUS */
1040 
1041  /* Passi fittizi successivi */
1042  if (iDummyStepsNumber > 1) {
1043  /* Setup SolutionManager(s) */
1045  }
1046 
1047  for (int iSubStep = 2;
1048  iSubStep <= iDummyStepsNumber;
1049  iSubStep++)
1050  {
1052  *qX[0], *qXPrime[0]);
1053  Flip();
1054 
1055  DEBUGLCOUT(MYDEBUG_FSTEPS, "Dummy step "
1056  << iSubStep
1057  << "; current time step: " << dCurrTimeStep
1058  << std::endl);
1059 
1060  if (outputStep()) {
1061  silent_cout("Dummy Step(" << iSubStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1062  }
1063 
1064  pCurrStepIntegrator = pDummySteps;
1065  ASSERT(pDummySteps!= 0);
1066  try {
1068  dTest = pDummySteps->Advance(this,
1069  dRefTimeStep,
1072  qX, qXPrime, pX, pXPrime,
1073  iStIter, dTest, dSolTest);
1074  }
1076  silent_cerr("Dummy step " << iSubStep
1077  << " does not converge; "
1078  "TimeStep=" << dCurrTimeStep
1079  << " cannot be reduced further; "
1080  "aborting..." << std::endl);
1081  pDM->Output(0, dTime, dCurrTimeStep, true);
1083  }
1084 
1086  /*
1087  * Mettere qui eventuali azioni speciali
1088  * da intraprendere in caso di errore ...
1089  */
1091  }
1092  catch (LinearSolver::ErrFactor& err) {
1093  /*
1094  * Mettere qui eventuali azioni speciali
1095  * da intraprendere in caso di errore ...
1096  */
1097  silent_cerr("Dummy step " << iSubStep
1098  << " failed because no pivot element "
1099  "could be found for column " << err.iCol
1100  << " (" << pDM->GetDofDescription(err.iCol) << "); "
1101  "aborting..." << std::endl);
1103  }
1105  bSolConv = true;
1106  }
1107  catch (EndOfSimulation& eos) {
1108  silent_cerr("Simulation ended during the dummy steps:\n"
1109  << eos.what() << "\n");
1110  return false;
1111  }
1112 
1113 #if 0
1114  /* don't sum up the derivatives error */
1115  dTotErr += dTest;
1116 #endif
1117  iTotIter += iStIter;
1118 
1119 #ifdef DEBUG
1120  if (DEBUG_LEVEL(MYDEBUG_FSTEPS)) {
1121  Out << "Step " << lStep
1122  << " time " << dTime + dCurrTimeStep
1123  << " step " << dCurrTimeStep
1124  << " iterations " << iStIter
1125  << " error " << dTest << std::endl;
1126  }
1127 #endif /* DEBUG */
1128 
1129  DEBUGLCOUT(MYDEBUG_FSTEPS, "Substep " << iSubStep
1130  << " of step " << lStep
1131  << " has been successfully completed "
1132  "in " << iStIter << " iterations"
1133  << std::endl);
1134 
1136  /* */
1137 #ifdef DEBUG_FICTITIOUS
1138  pDM->Output(0, dTime, dCurrTimeStep);
1139 #endif /* DEBUG_FICTITIOUS */
1140  Out << "Interrupted during dummy steps."
1141  << std::endl;
1143  }
1144 
1145  dTime += dRefTimeStep;
1146  }
1147  if (outputMsg()) {
1148  Out << "# Initial solution after dummy steps "
1149  "at time " << dTime
1150  << " performed in " << iStIter
1151  << " iterations with " << dTest
1152  << " error" << std::endl;
1153  }
1154 
1156  "Dummy steps have been successfully completed "
1157  "in " << iStIter << " iterations" << std::endl);
1158 #ifdef USE_EXTERNAL
1159  /* comunica che gli ultimi dati inviati sono la condizione iniziale */
1160  External::SendInitial();
1161 #endif /* USE_EXTERNAL */
1162  } /* Fine dei passi fittizi */
1163 
1164  /* Output delle "condizioni iniziali" */
1165  bOut = pDM->Output(0, dTime, dCurrTimeStep);
1166 
1167  if (outputMsg()) {
1168  Out
1169  << "# Key for lines starting with \"Step\":"
1170  << std::endl
1171  << "# Step Time TStep NIter ResErr SolErr SolConv Out"
1172  << std::endl
1173  << "Step " << 0
1174  << " " << dTime + dCurrTimeStep
1175  << " " << dCurrTimeStep
1176  << " " << iStIter
1177  << " " << dTest
1178  << " " << dSolTest
1179  << " " << bSolConv
1180  << " " << bOut
1181  << std::endl;
1182  }
1183 
1184 
1185  if (eAbortAfter == AFTER_DUMMY_STEPS) {
1186  Out << "End of dummy steps; no simulation is required."
1187  << std::endl;
1188  return false;
1189 
1190  } else if (mbdyn_stop_at_end_of_time_step()) {
1191  /* Fa l'output della soluzione ed esce */
1192  Out << "Interrupted during dummy steps." << std::endl;
1194  }
1195 
1197 
1198  return true;
1199 }
1200 
1201 bool
1203 {
1204  DEBUGCOUTFNAME("Solver::Start");
1205 
1206  // consistency check
1208  silent_cerr("Start() must be called after Prepare()" << std::endl);
1210  }
1211 
1212  /* primo passo regolare */
1213 
1214 #ifdef USE_EXTERNAL
1215  /* il prossimo passo e' un regular */
1216  pNLS->SetExternal(External::REGULAR);
1217 #endif /* USE_EXTERNAL */
1218 
1219  lStep = 1; /* Resetto di nuovo lStep */
1220 
1221  DEBUGCOUT("Step " << lStep << " has been successfully completed "
1222  "in " << iStIter << " iterations" << std::endl);
1223 
1224 
1225  DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
1226 
1227  pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
1228 
1229  Flip();
1232 
1235 
1236  /* Setup SolutionManager(s) */
1240 
1241 IfFirstStepIsToBeRepeated:
1242  try {
1244  if (outputStep()) {
1245  if (outputCounter()) {
1246  silent_cout(std::endl);
1247  }
1248  silent_cout("Step(" << 1 << ':' << 0 << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1249  }
1252  qX, qXPrime, pX, pXPrime,
1253  iStIter, dTest, dSolTest);
1254  }
1256  if (dCurrTimeStep > dMinTimeStep) {
1257  /* Riduce il passo */
1259  doublereal dOldCurrTimeStep = dCurrTimeStep;
1261  if (dCurrTimeStep < dOldCurrTimeStep) {
1262  DEBUGCOUT("Changing time step"
1263  " from " << dOldCurrTimeStep
1264  << " to " << dCurrTimeStep
1265  << " during first step after "
1266  << iStIter << " iterations"
1267  << std::endl);
1268  goto IfFirstStepIsToBeRepeated;
1269  }
1270  }
1271 
1272  silent_cerr("Max iterations number "
1274  << " has been reached during "
1275  "first step, Time=" << dTime << "; "
1276  << "TimeStep=" << dCurrTimeStep
1277  << " cannot be reduced further; "
1278  "aborting..." << std::endl);
1279  pDM->Output(0, dTime, dCurrTimeStep, true);
1280 
1282  }
1284  /*
1285  * Mettere qui eventuali azioni speciali
1286  * da intraprendere in caso di errore ...
1287  */
1288 
1290  }
1291  catch (LinearSolver::ErrFactor& err) {
1292  /*
1293  * Mettere qui eventuali azioni speciali
1294  * da intraprendere in caso di errore ...
1295  */
1296  silent_cerr("First step failed because no pivot element "
1297  "could be found for column " << err.iCol
1298  << " (" << pDM->GetDofDescription(err.iCol) << "); "
1299  "aborting..." << std::endl);
1301  }
1303  bSolConv = true;
1304  }
1305  catch (EndOfSimulation& eos) {
1306  silent_cerr("Simulation ended during the first regular step:\n"
1307  << eos.what() << "\n");
1308  return false;
1309  }
1310 
1312  pFirstRegularStep = 0;
1313 
1315 
1316  /* Si fa dare l'std::ostream al file di output per il log */
1317  std::ostream& Out = pDM->GetOutFile();
1318 
1320  /* Fa l'output della soluzione al primo passo ed esce */
1321  Out << "Interrupted during first step." << std::endl;
1323  }
1324 
1325  if (outputMsg()) {
1326  Out
1327  << "Step " << lStep
1328  << " " << dTime + dCurrTimeStep
1329  << " " << dCurrTimeStep
1330  << " " << iStIter
1331  << " " << dTest
1332  << " " << dSolTest
1333  << " " << bSolConv
1334  << " " << bOut
1335  << std::endl;
1336  }
1337 
1338  bSolConv = false;
1339 
1341  dTime += dRefTimeStep;
1342 
1343  dTotErr += dTest;
1344  iTotIter += iStIter;
1345 
1346  if (EigAn.bAnalysis
1347  && EigAn.currAnalysis != EigAn.Analyses.end()
1348  && *EigAn.currAnalysis <= dTime)
1349  {
1350  std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
1351  EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
1352  if (i != EigAn.Analyses.end()) {
1353  EigAn.currAnalysis = --i;
1354  }
1355  Eig();
1356  ++EigAn.currAnalysis;
1357  }
1358 
1359  if (pRTSolver) {
1360  pRTSolver->Init();
1361  }
1362 
1363  /* Altri passi regolari */
1364  ASSERT(pRegularSteps != 0);
1365 
1366  /* Setup SolutionManager(s) */
1368  pCurrStepIntegrator = pRegularSteps;
1369 
1371 
1372  return true;
1373 }
1374 
1375 bool
1377 {
1378  DEBUGCOUTFNAME("Solver::Advance");
1379 
1380  // consistency check
1381  if (eStatus != SOLVER_STATUS_STARTED) {
1382  silent_cerr("Started() must be called first" << std::endl);
1384  }
1385 
1387 
1388  if (pDM->EndOfSimulation() || dTime >= dFinalTime) {
1389  if (pRTSolver) {
1391  }
1392  silent_cout(outputCounterPrefix
1393  << "End of simulation at time "
1394  << dTime << " after "
1395  << lStep << " steps;" << std::endl
1396  << "output in file \"" << sOutputFileName << "\"" << std::endl
1397  << "total iterations: " << iTotIter << std::endl
1398  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1399  << "total error: " << dTotErr << std::endl);
1400 
1401  if (pRTSolver) {
1402  pRTSolver->Log();
1403  }
1404 
1405  return false;
1406 
1407  } else if (pRTSolver && pRTSolver->IsStopCommanded()) {
1408  silent_cout(outputCounterPrefix
1409  << "Simulation is stopped by RTAI task" << std::endl
1410  << "Simulation ended at time "
1411  << dTime << " after "
1412  << lStep << " steps;" << std::endl
1413  << "total iterations: " << iTotIter << std::endl
1414  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1415  << "total error: " << dTotErr << std::endl);
1416  pRTSolver->Log();
1417  return false;
1418 
1419  } else if (mbdyn_stop_at_end_of_time_step()
1420 #ifdef USE_MPI
1421  || (MPI_Finalized(&mpi_finalize), mpi_finalize)
1422 #endif /* USE_MPI */
1423  )
1424  {
1425  if (pRTSolver) {
1427  }
1428 
1429  silent_cout(outputCounterPrefix
1430  << "Interrupted!" << std::endl
1431  << "Simulation ended at time "
1432  << dTime << " after "
1433  << lStep << " steps;" << std::endl
1434  << "total iterations: " << iTotIter << std::endl
1435  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1436  << "total error: " << dTotErr << std::endl);
1437 
1438  if (pRTSolver) {
1439  pRTSolver->Log();
1440  }
1441 
1443  }
1444 
1445  lStep++;
1446  pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
1447 
1448  Flip();
1449 
1450  if (pRTSolver) {
1451  pRTSolver->Wait();
1452  }
1453 
1454  int retries = -1;
1455 IfStepIsToBeRepeated:
1456  try {
1457  retries++;
1459  if (outputStep()) {
1460  if (outputCounter()) {
1461  silent_cout(std::endl);
1462  }
1463  silent_cout("Step(" << lStep << ':' << retries << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
1464  }
1467  qX, qXPrime, pX, pXPrime, iStIter,
1468  dTest, dSolTest);
1469  }
1471  if (dCurrTimeStep > dMinTimeStep) {
1472  /* Riduce il passo */
1474  doublereal dOldCurrTimeStep = dCurrTimeStep;
1476  if (dCurrTimeStep < dOldCurrTimeStep) {
1477  DEBUGCOUT("Changing time step"
1478  " from " << dOldCurrTimeStep
1479  << " to " << dCurrTimeStep
1480  << " during step "
1481  << lStep << " after "
1482  << iStIter << " iterations"
1483  << std::endl);
1484  goto IfStepIsToBeRepeated;
1485  }
1486  }
1487 
1488  silent_cerr(outputCounterPrefix
1489  << "Max iterations number "
1490  << std::abs(pRegularSteps->GetIntegratorMaxIters())
1491  << " has been reached during "
1492  "Step=" << lStep << ", "
1493  "Time=" << dTime + dCurrTimeStep << "; "
1494  "TimeStep=" << dCurrTimeStep
1495  << " cannot be reduced further; "
1496  "aborting..." << std::endl);
1498  }
1500  if (dCurrTimeStep > dMinTimeStep) {
1501  /* Riduce il passo */
1503  doublereal dOldCurrTimeStep = dCurrTimeStep;
1505  if (dCurrTimeStep < dOldCurrTimeStep) {
1506  DEBUGCOUT("Changing time step"
1507  " from " << dOldCurrTimeStep
1508  << " to " << dCurrTimeStep
1509  << " during step "
1510  << lStep << " after "
1511  << iStIter << " iterations"
1512  << std::endl);
1513  goto IfStepIsToBeRepeated;
1514  }
1515  }
1516 
1517  silent_cerr(outputCounterPrefix
1518  << "Simulation diverged after "
1519  << iStIter << " iterations, before "
1520  "reaching max iteration number "
1521  << std::abs(pRegularSteps->GetIntegratorMaxIters())
1522  << " during Step=" << lStep << ", "
1523  "Time=" << dTime + dCurrTimeStep << "; "
1524  "TimeStep=" << dCurrTimeStep
1525  << " cannot be reduced further; "
1526  "aborting..." << std::endl);
1528  }
1529  catch (LinearSolver::ErrFactor& err) {
1530  /*
1531  * Mettere qui eventuali azioni speciali
1532  * da intraprendere in caso di errore ...
1533  */
1534  silent_cerr(outputCounterPrefix
1535  << "Simulation failed because no pivot element "
1536  "could be found for column " << err.iCol
1537  << " (" << pDM->GetDofDescription(err.iCol) << ") "
1538  "after " << iStIter << " iterations "
1539  "during step " << lStep << "; "
1540  "aborting..." << std::endl);
1542  }
1544  bSolConv = true;
1545  }
1546  catch (EndOfSimulation& eos) {
1547  silent_cerr(outputCounterPrefix
1548  << "Simulation ended during a regular step:\n"
1549  << eos.what() << "\n");
1550 #ifdef USE_MPI
1551  MBDynComm.Abort(0);
1552 #endif /* USE_MPI */
1553  if (pRTSolver) {
1555  }
1556 
1557  silent_cout("Simulation ended at time "
1558  << dTime << " after "
1559  << lStep << " steps;" << std::endl
1560  << "total iterations: " << iTotIter << std::endl
1561  << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
1562  << "total error: " << dTotErr << std::endl);
1563 
1564  if (pRTSolver) {
1565  pRTSolver->Log();
1566  }
1567 
1568  return false;
1569  }
1570 
1571  dTotErr += dTest;
1572  iTotIter += iStIter;
1573 
1575 
1576  /* Si fa dare l'std::ostream al file di output per il log */
1577  std::ostream& Out = pDM->GetOutFile();
1578 
1579  if (outputMsg()) {
1580  Out << "Step " << lStep
1581  << " " << dTime + dCurrTimeStep
1582  << " " << dCurrTimeStep
1583  << " " << iStIter
1584  << " " << dTest
1585  << " " << dSolTest
1586  << " " << bSolConv
1587  << " " << bOut
1588  << std::endl;
1589  }
1590 
1591  if (bOutputCounter) {
1592  silent_cout("Step " << std::setw(5) << lStep
1593  << " " << std::setw(13) << dTime + dCurrTimeStep
1594  << " " << std::setw(13) << dCurrTimeStep
1595  << " " << std::setw(4) << iStIter
1596  << " " << std::setw(13) << dTest
1597  << " " << std::setw(13) << dSolTest
1598  << " " << bSolConv
1599  << " " << bOut
1601  }
1602 
1603  DEBUGCOUT("Step " << lStep
1604  << " has been successfully completed "
1605  "in " << iStIter << " iterations" << std::endl);
1606 
1608  dTime += dRefTimeStep;
1609 
1610  bSolConv = false;
1611 
1612  if (EigAn.bAnalysis
1613  && EigAn.currAnalysis != EigAn.Analyses.end()
1614  && *EigAn.currAnalysis <= dTime)
1615  {
1616  std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
1617  EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
1618  if (i != EigAn.Analyses.end()) {
1619  EigAn.currAnalysis = --i;
1620  }
1622  ++EigAn.currAnalysis;
1623  }
1624 
1625  /* Calcola il nuovo timestep */
1627  DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
1628 
1629  return true;
1630 }
1631 
1632 void
1634 {
1635  DEBUGCOUTFNAME("Solver::Run");
1636 
1637  if (Prepare()) {
1638  if (Start()) {
1639  while (Advance()) {
1640  NO_OP;
1641  }
1642  }
1643  }
1644 }
1645 
1646 /* Distruttore */
1648 {
1649  DEBUGCOUTFNAME("Solver::~Solver");
1650 
1651  if (!qX.empty()) {
1652  for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
1653  if (qX[ivec] != 0) {
1654  SAFEDELETE(qX[ivec]);
1655  SAFEDELETE(qXPrime[ivec]);
1656  }
1657  }
1658  }
1659 
1660  if (pX) {
1661  SAFEDELETE(pX);
1662  }
1663 
1664  if (pXPrime) {
1666  }
1667 
1668  if (pdWorkSpace) {
1670  }
1671 
1672  if (pDM) {
1673  SAFEDELETE(pDM);
1674  }
1675 
1676  if (pRTSolver) {
1678  }
1679 
1680  if (pDerivativeSteps) {
1682  }
1683 
1684  if (pFirstDummyStep) {
1686  }
1687 
1688  if (pDummySteps) {
1690  }
1691 
1692  if (pFirstRegularStep) {
1694  }
1695 
1696  if (pRegularSteps) {
1698  }
1699 
1700  if (pSM) {
1701  SAFEDELETE(pSM);
1702  }
1703 
1704  if (pNLS) {
1705  SAFEDELETE(pNLS);
1706  }
1707 
1708  if (pTSC) {
1709  delete pTSC;
1710  }
1711 
1713 }
1714 
1715 /*scrive il contributo al file di restart*/
1716 std::ostream &
1717 Solver::Restart(std::ostream& out,DataManager::eRestart type) const
1718 {
1719 
1720  out << "begin: initial value;" << std::endl;
1721  switch(type) {
1722  case DataManager::ATEND:
1723  out << " # initial time: " << pDM->dGetTime() << ";"
1724  << std::endl
1725  << " # final time: " << dFinalTime << ";"
1726  << std::endl
1727  << " # time step: " << dInitialTimeStep << ";"
1728  << std::endl;
1729  break;
1731  case DataManager::TIME:
1732  case DataManager::TIMES:
1733  out << " initial time: " << pDM->dGetTime()<< ";" << std::endl
1734  << " final time: " << dFinalTime << ";" << std::endl
1735  << " time step: " << dInitialTimeStep << ";"
1736  << std::endl;
1737  break;
1738  default:
1739  ASSERT(0);
1740  }
1741 
1742  out << " method: ";
1743  switch(RegularType) {
1744  case INT_CRANKNICOLSON:
1745  out << "Crank Nicolson; " << std::endl;
1746  break;
1747  case INT_MS2:
1748  out << "ms, ";
1749  pRhoRegular->Restart(out) << ", ";
1750  pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
1751  break;
1752  case INT_HOPE:
1753  out << "hope, ";
1754  pRhoRegular->Restart(out) << ", ";
1755  pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
1756  break;
1757 
1758  case INT_THIRDORDER:
1759  out << "thirdorder, ";
1760  if (!pRhoRegular)
1761  out << "ad hoc;" << std::endl;
1762  else
1763  pRhoRegular->Restart(out) << ";" << std::endl;
1764  break;
1765  case INT_IMPLICITEULER:
1766  out << "implicit euler;" << std::endl;
1767  break;
1768  default:
1769  ASSERT(0);
1770  }
1771 
1773  out << " max iterations: " << std::abs(iMI);
1774  if (iMI < 0) out << ", at most";
1775  out
1776  << ";" << std::endl
1777  << " tolerance: " << pRegularSteps->GetIntegratorDTol();
1778  switch(ResTest) {
1780  out << ", test, norm" ;
1781  break;
1783  out << ", test, minmax" ;
1784  break;
1786  NO_OP;
1787  default:
1788  silent_cerr("unhandled nonlinear solver test type" << std::endl);
1790  }
1791 
1792  if (bScale) {
1793  out << ", scale"
1794  << ", " << pRegularSteps->GetIntegratorDSolTol();
1795  }
1796 
1797  switch (SolTest) {
1799  out << ", test, norm" ;
1800  break;
1802  out << ", test, minmax" ;
1803  break;
1805  NO_OP;
1806  default:
1807  ASSERT(0);
1808  }
1809  out
1810  << ";" << std::endl;
1811 
1812  if ( pDerivativeSteps != 0 ) {
1813  out
1814  << " derivatives max iterations: " << pDerivativeSteps->GetIntegratorMaxIters() << ";" << std::endl
1815  << " derivatives tolerance: " << pDerivativeSteps->GetIntegratorDTol() << ";" << std::endl
1816  << " derivatives coefficient: " << dDerivativesCoef << ";" << std::endl;
1817  }
1818 
1819  if (iDummyStepsNumber) {
1820  out
1821  << " dummy steps max iterations: " << pDummySteps->GetIntegratorMaxIters() << ";" << std::endl
1822  << " dummy steps tolerance: " << pDummySteps->GetIntegratorDTol() << ";" << std::endl;
1823  }
1824  out
1825  << " dummy steps number: " << iDummyStepsNumber << ";" << std::endl
1826  << " dummy steps ratio: " << dDummyStepsRatio << ";" << std::endl;
1827  switch (NonlinearSolverType) {
1829  out << " # nonlinear solver: matrix free;" << std::endl;
1830  break;
1832  default :
1833  out << " nonlinear solver: newton raphson";
1834  if (!bTrueNewtonRaphson) {
1835  out << ", modified, " << iIterationsBeforeAssembly;
1836  if (bKeepJac) {
1837  out << ", keep jacobian matrix";
1838  }
1839  if (bHonorJacRequest) {
1840  out << ", honor element requests";
1841  }
1842  }
1843  out << ";" << std::endl;
1844  }
1845  out << " solver: ";
1847  out << "end: initial value;" << std::endl << std::endl;
1848  return out;
1849 }
1850 
1851 /* Dati dell'integratore */
1852 void
1854 {
1855  DEBUGCOUTFNAME("MultiStepIntegrator::ReadData");
1856 
1857  /* parole chiave */
1858  static const char*const sKeyWords[] = {
1859  "begin",
1860  "initial" "value",
1861  "multistep", /* deprecated */
1862  "end",
1863 
1864  "initial" "time",
1865  "final" "time",
1866  "time" "step",
1867  "min" "time" "step",
1868  "max" "time" "step",
1869  "tolerance",
1870  "max" "residual",
1871  "max" "iterations",
1872  "modify" "residual" "test",
1873  "enforce" "constraint" "equations",
1874  /* DEPRECATED */
1875  "fictitious" "steps" "number",
1876  "fictitious" "steps" "ratio",
1877  "fictitious" "steps" "tolerance",
1878  "fictitious" "steps" "max" "iterations",
1879  /* END OF DEPRECATED */
1880 
1881  "dummy" "steps" "number",
1882  "dummy" "steps" "ratio",
1883  "dummy" "steps" "tolerance",
1884  "dummy" "steps" "max" "iterations",
1885 
1886  "abort" "after",
1887  "input",
1888  "assembly",
1889  "derivatives",
1890 
1891  /* DEPRECATED */ "fictitious" "steps" /* END OF DEPRECATED */ ,
1892  "dummy" "steps",
1893 
1894  "output",
1895  "none",
1896  "iterations",
1897  "residual",
1898  "solution",
1899  /* DEPRECATED */ "jacobian" /* END OF DEPRECATED */ ,
1900  "jacobian" "matrix",
1901  "bailout",
1902  "messages",
1903  "counter",
1904  "matrix" "condition" "number",
1905  "solver" "condition" "number",
1906  "cpu" "time",
1907  "output" "meter",
1908 
1909  "method",
1910  /* DEPRECATED */ "fictitious" "steps" "method" /* END OF DEPRECATED */ ,
1911  "dummy" "steps" "method",
1912 
1913  "Crank" "Nicolson",
1914  /* DEPRECATED */ "Crank" "Nicholson" /* END OF DEPRECATED */ ,
1915  /* DEPRECATED */ "nostro" /* END OF DEPRECATED */ ,
1916  "ms",
1917  "hope",
1918  "bdf",
1919  "thirdorder",
1920  "implicit" "euler",
1921 
1922  "derivatives" "coefficient",
1923  "derivatives" "tolerance",
1924  "derivatives" "max" "iterations",
1925 
1926  /* DEPRECATED */
1927  "true",
1928  "modified",
1929  /* END OF DEPRECATED */
1930 
1931  "strategy",
1932  "factor",
1933  "no" "change",
1934  "change",
1935 
1936  "pod",
1937  "eigen" "analysis",
1938 
1939  /* DEPRECATED */
1940  "solver",
1941  "interface" "solver",
1942  /* END OF DEPRECATED */
1943  "linear" "solver",
1944  "interface" "linear" "solver",
1945 
1946  /* DEPRECATED */
1947  "preconditioner",
1948  /* END OF DEPRECATED */
1949 
1950  "nonlinear" "solver",
1951  "default",
1952  "newton" "raphson",
1953  "line" "search",
1954  "matrix" "free",
1955  "bicgstab",
1956  "gmres",
1957  /* DEPRECATED */ "full" "jacobian" /* END OF DEPRECATED */ ,
1958  "full" "jacobian" "matrix",
1959 
1960  /* RTAI stuff */
1961  "real" "time",
1962 
1963  /* multithread stuff */
1964  "threads",
1965 
1966  NULL
1967  };
1968 
1969  /* enum delle parole chiave */
1970  enum KeyWords {
1971  UNKNOWN = -1,
1972  BEGIN = 0,
1973  INITIAL_VALUE,
1974  MULTISTEP,
1975  END,
1976 
1977  INITIALTIME,
1978  FINALTIME,
1979  TIMESTEP,
1980  MINTIMESTEP,
1981  MAXTIMESTEP,
1982  TOLERANCE,
1983  MAXRESIDUAL,
1984  MAXITERATIONS,
1985  MODIFY_RES_TEST,
1986  ENFORCE_CONSTRAINT_EQUATIONS,
1987  FICTITIOUSSTEPSNUMBER,
1988  FICTITIOUSSTEPSRATIO,
1989  FICTITIOUSSTEPSTOLERANCE,
1990  FICTITIOUSSTEPSMAXITERATIONS,
1991 
1992  DUMMYSTEPSNUMBER,
1993  DUMMYSTEPSRATIO,
1994  DUMMYSTEPSTOLERANCE,
1995  DUMMYSTEPSMAXITERATIONS,
1996 
1997  ABORTAFTER,
1998  INPUT,
1999  ASSEMBLY,
2000  DERIVATIVES,
2001  FICTITIOUSSTEPS,
2002  DUMMYSTEPS,
2003 
2004  OUTPUT,
2005  NONE,
2006  ITERATIONS,
2007  RESIDUAL,
2008  SOLUTION,
2009  JACOBIAN,
2010  JACOBIANMATRIX,
2011  BAILOUT,
2012  MESSAGES,
2013  COUNTER,
2014  MATRIX_COND_NUM,
2015  SOLVER_COND_NUM,
2016  CPU_TIME,
2017  OUTPUTMETER,
2018 
2019  METHOD,
2020  FICTITIOUSSTEPSMETHOD,
2021  DUMMYSTEPSMETHOD,
2022  CRANKNICOLSON,
2023  CRANKNICHOLSON,
2024  NOSTRO,
2025  MS,
2026  HOPE,
2027  BDF,
2028  THIRDORDER,
2029  IMPLICITEULER,
2030 
2031  DERIVATIVESCOEFFICIENT,
2032  DERIVATIVESTOLERANCE,
2033  DERIVATIVESMAXITERATIONS,
2034 
2035  /* DEPRECATED */
2036  NR_TRUE,
2037  MODIFIED,
2038  /* END OF DEPRECATED */
2039 
2040  STRATEGY,
2041  STRATEGYFACTOR,
2042  STRATEGYNOCHANGE,
2043  STRATEGYCHANGE,
2044 
2045  POD,
2046  EIGENANALYSIS,
2047 
2048  /* DEPRECATED */
2049  SOLVER,
2050  INTERFACESOLVER,
2051  /* END OF DEPRECATED */
2052  LINEARSOLVER,
2053  INTERFACELINEARSOLVER,
2054 
2055  /* DEPRECATED */
2056  PRECONDITIONER,
2057  /* END OF DEPRECATED */
2058 
2059  NONLINEARSOLVER,
2060  DEFAULT,
2061  NEWTONRAPHSON,
2062  LINESEARCH,
2063  MATRIXFREE,
2064  BICGSTAB,
2065  GMRES,
2066  FULLJACOBIAN,
2067  FULLJACOBIANMATRIX,
2068 
2069  /* RTAI stuff */
2070  REALTIME,
2071 
2072  THREADS,
2073 
2074  LASTKEYWORD
2075  };
2076 
2077  /* tabella delle parole chiave */
2078  KeyTable K(HP, sKeyWords);
2079 
2080  /* legge i dati della simulazione */
2081  if (KeyWords(HP.GetDescription()) != BEGIN) {
2082  silent_cerr("Error: <begin> expected at line "
2083  << HP.GetLineData() << "; aborting..." << std::endl);
2085  }
2086 
2087  switch (KeyWords(HP.GetWord())) {
2088  case MULTISTEP:
2089  pedantic_cout("warning: \"begin: multistep\" is deprecated; "
2090  "use \"begin: initial value;\" instead." << std::endl);
2091  case INITIAL_VALUE:
2092  break;
2093 
2094  default:
2095  silent_cerr("Error: \"begin: initial value;\" expected at line "
2096  << HP.GetLineData() << "; aborting..." << std::endl);
2098  }
2099 
2100  bool bMethod(false);
2101  bool bDummyStepsMethod(false);
2102 
2103  /* dati letti qui ma da passare alle classi
2104  * StepIntegration e NonlinearSolver
2105  */
2106 
2107  doublereal dTol = ::dDefaultTol;
2108  doublereal dSolutionTol = 0.;
2110  bool bModResTest = false;
2111  bool bSetScaleAlgebraic = false;
2112 
2113  /* Dati dei passi fittizi di trimmaggio iniziale */
2114  doublereal dDummyStepsTolerance = ::dDefaultDummyStepsTolerance;
2115  integer iDummyStepsMaxIterations = ::iDefaultMaxIterations;
2116 
2117  /* Dati del passo iniziale di calcolo delle derivate */
2118 
2119  doublereal dDerivativesTol = ::dDefaultTol;
2120  integer iDerivativesMaxIterations = ::iDefaultMaxIterations;
2121  integer iDerivativesCoefMaxIter = 0;
2122  doublereal dDerivativesCoefFactor = 10.;
2123 
2124 #ifdef USE_MULTITHREAD
2125  bool bSolverThreads(false);
2126  unsigned nSolverThreads = 0;
2127 #endif /* USE_MULTITHREAD */
2128 
2129 
2130  /* Ciclo infinito */
2131  while (true) {
2132  KeyWords CurrKeyWord = KeyWords(HP.GetDescription());
2133 
2134  switch (CurrKeyWord) {
2135  case INITIALTIME:
2136  dInitialTime = HP.GetReal();
2137  DEBUGLCOUT(MYDEBUG_INPUT, "Initial time is "
2138  << dInitialTime << std::endl);
2139  break;
2140 
2141  case FINALTIME:
2142  if (HP.IsKeyWord("forever")) {
2143  dFinalTime = std::numeric_limits<doublereal>::max();
2144  } else {
2145  dFinalTime = HP.GetReal();
2146  }
2147  DEBUGLCOUT(MYDEBUG_INPUT, "Final time is "
2148  << dFinalTime << std::endl);
2149 
2150  if(dFinalTime <= dInitialTime) {
2151  silent_cerr("warning: final time " << dFinalTime
2152  << " is less than initial time "
2153  << dInitialTime << ';' << std::endl
2154  << "this will cause the simulation"
2155  " to abort" << std::endl);
2156  }
2157  break;
2158 
2159  case TIMESTEP:
2160  dInitialTimeStep = HP.GetReal();
2161  DEBUGLCOUT(MYDEBUG_INPUT, "Initial time step is "
2162  << dInitialTimeStep << std::endl);
2163 
2164  if (dInitialTimeStep == 0.) {
2165  silent_cerr("warning, null initial time step"
2166  " is not allowed" << std::endl);
2167  } else if (dInitialTimeStep < 0.) {
2169  silent_cerr("warning, negative initial time step"
2170  " is not allowed;" << std::endl
2171  << "its modulus " << dInitialTimeStep
2172  << " will be considered" << std::endl);
2173  }
2174  break;
2175 
2176  case MINTIMESTEP:
2177  try {
2178  dMinTimeStep = HP.GetReal(std::numeric_limits<doublereal>::min(), HighParser::range_gt<doublereal>(0.));
2179 
2181  silent_cerr("invalid min time step " << e.Get() << " (must be positive) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
2182  throw e;
2183  }
2184  break;
2185 
2186  case MAXTIMESTEP:
2187  if (HP.IsKeyWord("unlimited")) {
2188  DriveCaller *pDC = 0;
2190  MaxTimeStep.Set(pDC);
2191 
2192  } else {
2194  }
2195 
2196 #ifdef DEBUG
2197  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
2198  DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is postponed" << std::endl);
2199 
2200  } else {
2201  DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is " << MaxTimeStep.dGet() << std::endl);
2202  }
2203 #endif // DEBUG
2204 
2205  eTimeStepLimit = (HP.IsKeyWord("hard" "limit")
2206  && HP.GetYesNoOrBool())
2207  ? TS_HARD_LIMIT
2208  : TS_SOFT_LIMIT;
2209 
2210  if (dGetInitialMaxTimeStep() == 0.) {
2211  silent_cout("no max time step limit will be"
2212  " considered" << std::endl);
2213 
2214  } else if (dGetInitialMaxTimeStep() < 0.) {
2215  silent_cerr("negative max time step"
2216  " is not allowed" << std::endl);
2218  }
2219  break;
2220 
2221  case FICTITIOUSSTEPSNUMBER:
2222  case DUMMYSTEPSNUMBER:
2223  iDummyStepsNumber = HP.GetInt();
2224  if (iDummyStepsNumber < 0) {
2227  silent_cerr("negative dummy steps number"
2228  " is illegal" << std::endl);
2230 
2231  } else if (iDummyStepsNumber == 1) {
2232  silent_cerr("warning, a single dummy step"
2233  " may be useless" << std::endl);
2234  }
2235 
2236  DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps number: "
2237  << iDummyStepsNumber << std::endl);
2238  break;
2239 
2240  case FICTITIOUSSTEPSRATIO:
2241  case DUMMYSTEPSRATIO:
2242  dDummyStepsRatio = HP.GetReal();
2243  if (dDummyStepsRatio < 0.) {
2244  silent_cerr("negative dummy steps ratio"
2245  " is illegal" << std::endl);
2247  }
2248 
2249  if (dDummyStepsRatio > 1.) {
2250  silent_cerr("warning, dummy steps ratio"
2251  " is larger than one." << std::endl
2252  << "Something like 1.e-3 should"
2253  " be safer ..." << std::endl);
2254  }
2255 
2256  DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps ratio: "
2257  << dDummyStepsRatio << std::endl);
2258  break;
2259 
2260  case FICTITIOUSSTEPSTOLERANCE:
2261  case DUMMYSTEPSTOLERANCE:
2262  dDummyStepsTolerance = HP.GetReal();
2263  if (dDummyStepsTolerance <= 0.) {
2264  silent_cerr("negative dummy steps"
2265  " tolerance is illegal" << std::endl);
2267  }
2269  "Dummy steps tolerance: "
2270  << dDummyStepsTolerance << std::endl);
2271  break;
2272 
2273  case ABORTAFTER: {
2274  KeyWords WhenToAbort(KeyWords(HP.GetWord()));
2275  switch (WhenToAbort) {
2276  case INPUT:
2279  "Simulation will abort after"
2280  " data input" << std::endl);
2281  break;
2282 
2283  case ASSEMBLY:
2286  "Simulation will abort after"
2287  " initial assembly" << std::endl);
2288  break;
2289 
2290  case DERIVATIVES:
2293  "Simulation will abort after"
2294  " derivatives solution" << std::endl);
2295  break;
2296 
2297  case FICTITIOUSSTEPS:
2298  case DUMMYSTEPS:
2301  "Simulation will abort after"
2302  " dummy steps solution" << std::endl);
2303  break;
2304 
2305  default:
2306  silent_cerr("Don't know when to abort,"
2307  " so I'm going to abort now" << std::endl);
2309  }
2310  } break;
2311 
2312  case OUTPUT: {
2313  unsigned OF = OUTPUT_DEFAULT;
2314  bool setOutput = false;
2315 
2316  while (HP.IsArg()) {
2317  KeyWords OutputFlag(KeyWords(HP.GetWord()));
2318  switch (OutputFlag) {
2319  case NONE:
2320  OF = OUTPUT_NONE;
2321  setOutput = true;
2322  break;
2323 
2324  case ITERATIONS:
2325  OF |= OUTPUT_ITERS;
2326  break;
2327 
2328  case RESIDUAL:
2329  OF |= OUTPUT_RES;
2330  break;
2331 
2332  case SOLUTION:
2333  OF |= OUTPUT_SOL;
2334  break;
2335 
2336  case JACOBIAN:
2337  case JACOBIANMATRIX:
2338  OF |= OUTPUT_JAC;
2339  break;
2340 
2341  case BAILOUT:
2342  OF |= OUTPUT_BAILOUT;
2343  break;
2344 
2345  case MESSAGES:
2346  OF |= OUTPUT_MSG;
2347  break;
2348 
2349  case COUNTER:
2350  OF |= OUTPUT_COUNTER;
2351  break;
2352 
2353  case MATRIX_COND_NUM:
2355  if (HP.IsKeyWord("norm")) {
2356  if (HP.IsKeyWord("inf")) {
2358  } else {
2359  const double P = HP.GetReal();
2360  if (P != 1.) {
2361  silent_cerr("Only inf or 1 norm are supported for condition numbers at line " << HP.GetLineData() << std::endl);
2363  }
2364  OF |= OUTPUT_MAT_COND_NUM_1;
2365  }
2366  } else {
2367  OF |= OUTPUT_MAT_COND_NUM_1;
2368  }
2369  break;
2370 
2371  case SOLVER_COND_NUM:
2372  OF |= OUTPUT_SOLVER_COND_NUM;
2373  if (HP.IsKeyWord("stat")) {
2374  if (HP.GetYesNoOrBool()) {
2376  } else {
2378  }
2379  }
2380  break;
2381  case CPU_TIME:
2382  OF |= OUTPUT_CPU_TIME;
2383  break;
2384 
2385  default:
2386  silent_cerr("Warning, unknown output flag "
2387  "at line " << HP.GetLineData()
2388  << "; ignored" << std::endl);
2389  break;
2390  }
2391  }
2392 
2393  if (setOutput) {
2394  SetOutputFlags(OF);
2395 
2396  } else {
2397  AddOutputFlags(OF);
2398  }
2399  } break;
2400 
2401  case OUTPUTMETER:
2402  SetOutputMeter(HP.GetDriveCaller(true));
2403  break;
2404 
2405  case METHOD: {
2406  if (bMethod) {
2407  silent_cerr("error: multiple definition"
2408  " of integration method at line "
2409  << HP.GetLineData() << std::endl);
2411  }
2412  bMethod = true;
2413 
2414  KeyWords KMethod = KeyWords(HP.GetWord());
2415  switch (KMethod) {
2416  case CRANKNICHOLSON:
2417  silent_cout("warning: \"crank nicolson\" is the correct spelling "
2418  "at line " << HP.GetLineData() << std::endl);
2419  case CRANKNICOLSON:
2421  break;
2422 
2423  case BDF:
2424  /* default (order 2) */
2425  RegularType = INT_MS2;
2426 
2427  if (HP.IsKeyWord("order")) {
2428  int iOrder = HP.GetInt();
2429 
2430  switch (iOrder) {
2431  case 1:
2433  break;
2434 
2435  case 2:
2436  break;
2437 
2438  default:
2439  silent_cerr("unhandled BDF order " << iOrder << std::endl);
2441  }
2442  }
2443 
2444  if (RegularType == INT_MS2) {
2447  }
2448  break;
2449 
2450  case NOSTRO:
2451  silent_cerr("integration method \"nostro\" "
2452  "is deprecated; use \"ms\" "
2453  "instead at line "
2454  << HP.GetLineData()
2455  << std::endl);
2456  case MS:
2457  case HOPE:
2458  pRhoRegular = HP.GetDriveCaller(true);
2459 
2461  if (HP.IsArg()) {
2463  } else {
2465  }
2466 
2467  switch (KMethod) {
2468  case NOSTRO:
2469  case MS:
2470  RegularType = INT_MS2;
2471  break;
2472 
2473  case HOPE:
2475  break;
2476 
2477  default:
2479  }
2480  break;
2481 
2482  case THIRDORDER:
2483  if (HP.IsKeyWord("ad" "hoc")) {
2484  /* do nothing */ ;
2485  } else {
2486  pRhoRegular = HP.GetDriveCaller(true);
2487  }
2489  break;
2490 
2491  case IMPLICITEULER:
2493  break;
2494 
2495  default:
2496  silent_cerr("Unknown integration method at line "
2497  << HP.GetLineData() << std::endl);
2499  }
2500  break;
2501  }
2502 
2503  case FICTITIOUSSTEPSMETHOD:
2504  case DUMMYSTEPSMETHOD: {
2505  if (bDummyStepsMethod) {
2506  silent_cerr("error: multiple definition "
2507  "of dummy steps integration method "
2508  "at line " << HP.GetLineData()
2509  << std::endl);
2511  }
2512  bDummyStepsMethod = true;
2513 
2514  KeyWords KMethod = KeyWords(HP.GetWord());
2515  switch (KMethod) {
2516  case CRANKNICHOLSON:
2517  silent_cout("warning: \"crank nicolson\" is the correct spelling" << std::endl);
2518  case CRANKNICOLSON:
2520  break;
2521 
2522  case BDF:
2523  /* default (order 2) */
2524  DummyType = INT_MS2;
2525 
2526  if (HP.IsKeyWord("order")) {
2527  int iOrder = HP.GetInt();
2528 
2529  switch (iOrder) {
2530  case 1:
2532  break;
2533 
2534  case 2:
2535  break;
2536 
2537  default:
2538  silent_cerr("unhandled BDF order " << iOrder << std::endl);
2540  }
2541  }
2542 
2543  if (DummyType == INT_MS2) {
2546  }
2547  break;
2548 
2549  case NOSTRO:
2550  case MS:
2551  case HOPE:
2552  pRhoDummy = HP.GetDriveCaller(true);
2553 
2554  if (HP.IsArg()) {
2556  } else {
2558  }
2559 
2560  switch (KMethod) {
2561  case NOSTRO:
2562  case MS:
2563  DummyType = INT_MS2;
2564  break;
2565 
2566  case HOPE:
2567  DummyType = INT_HOPE;
2568  break;
2569 
2570  default:
2572  }
2573  break;
2574 
2575  case THIRDORDER:
2576  if (HP.IsKeyWord("ad" "hoc")) {
2577  /* do nothing */ ;
2578  } else {
2579  pRhoDummy = HP.GetDriveCaller(true);
2580  }
2582  break;
2583 
2584  case IMPLICITEULER:
2586  break;
2587 
2588  default:
2589  silent_cerr("Unknown integration method at line " << HP.GetLineData() << std::endl);
2591  }
2592  break;
2593  }
2594 
2595  case TOLERANCE: {
2596  /*
2597  * residual tolerance; can be the keyword "null",
2598  * which means that the convergence test
2599  * will be computed on the solution, or a number
2600  */
2601  if (HP.IsKeyWord("null")) {
2602  dTol = 0.;
2603 
2604  } else {
2605  dTol = HP.GetReal();
2606  if (dTol < 0.) {
2607  dTol = ::dDefaultTol;
2608  silent_cerr("warning, residual tolerance "
2609  "< 0. is illegal; "
2610  "using default value " << dTol
2611  << std::endl);
2612  }
2613  }
2614 
2615  /* safe default */
2616  if (dTol == 0.) {
2618  }
2619 
2620  if (HP.IsArg()) {
2621  if (HP.IsKeyWord("test")) {
2622  if (HP.IsKeyWord("norm")) {
2624  } else if (HP.IsKeyWord("minmax")) {
2626  } else if (HP.IsKeyWord("none")) {
2628  } else {
2629  silent_cerr("unknown test "
2630  "method at line "
2631  << HP.GetLineData()
2632  << std::endl);
2634  }
2635 
2636  if (HP.IsKeyWord("scale")) {
2638  silent_cerr("it's a nonsense "
2639  "to scale a disabled test; "
2640  "\"scale\" ignored"
2641  << std::endl);
2642  bScale = false;
2643  } else {
2644  bScale = true;
2645  }
2646  }
2647  }
2648  }
2649 
2650  if (HP.IsArg()) {
2651  if (!HP.IsKeyWord("null")) {
2652  dSolutionTol = HP.GetReal();
2653  }
2654 
2655  /* safe default */
2656  if (dSolutionTol != 0.) {
2658  }
2659 
2660  if (HP.IsArg()) {
2661  if (HP.IsKeyWord("test")) {
2662  if (HP.IsKeyWord("norm")) {
2664  } else if (HP.IsKeyWord("minmax")) {
2666  } else if (HP.IsKeyWord("none")) {
2668  } else {
2669  silent_cerr("unknown test "
2670  "method at line "
2671  << HP.GetLineData()
2672  << std::endl);
2674  }
2675  }
2676  }
2677 
2678  } else if (dTol == 0.) {
2679  silent_cerr("need solution tolerance "
2680  "with null residual tolerance"
2681  << std::endl);
2683  }
2684 
2685  if (dSolutionTol < 0.) {
2686  dSolutionTol = 0.;
2687  silent_cerr("warning, solution tolerance "
2688  "< 0. is illegal; "
2689  "solution test is disabled"
2690  << std::endl);
2691  }
2692 
2693  if (dTol == 0. && dSolutionTol == 0.) {
2694  silent_cerr("both residual and solution "
2695  "tolerances are zero" << std::endl);
2697  }
2698 
2699  DEBUGLCOUT(MYDEBUG_INPUT, "tolerance = " << dTol
2700  << ", " << dSolutionTol << std::endl);
2701  break;
2702  }
2703 
2704  case MAXRESIDUAL: {
2705  if (HP.IsKeyWord("differential" "equations")) {
2706  dMaxResidualDiff = HP.GetReal();
2707 
2708  if (dMaxResidualDiff <= 0.) {
2709  silent_cerr("error: max residual for differential equations "
2710  "must be greater than zero at line "
2711  << HP.GetLineData() << std::endl);
2713  }
2714  }
2715 
2716  if (HP.IsKeyWord("all" "equations") || HP.IsArg()) {
2717  dMaxResidual = HP.GetReal();
2718 
2719  if (dMaxResidual <= 0.) {
2720  silent_cerr("error: max residual must be greater than zero at line "
2721  << HP.GetLineData() << std::endl);
2723  }
2724  }
2725  break;
2726  }
2727 
2728  case DERIVATIVESTOLERANCE: {
2729  dDerivativesTol = HP.GetReal();
2730  if (dDerivativesTol <= 0.) {
2731  dDerivativesTol = ::dDefaultTol;
2732  silent_cerr("warning, derivatives "
2733  "tolerance <= 0.0 is illegal; "
2734  "using default value "
2735  << dDerivativesTol
2736  << std::endl);
2737  }
2739  "Derivatives tolerance = "
2740  << dDerivativesTol
2741  << std::endl);
2742  break;
2743  }
2744 
2745  case MAXITERATIONS: {
2746  iMaxIterations = HP.GetInt();
2747  if (iMaxIterations < 1) {
2749  silent_cerr("warning, max iterations "
2750  "< 1 is illegal; using default value "
2751  << iMaxIterations
2752  << std::endl);
2753  }
2755  "Max iterations = "
2756  << iMaxIterations << std::endl);
2757  if (HP.IsKeyWord("at" "most")) {
2759  }
2760  break;
2761  }
2762 
2763  case MODIFY_RES_TEST:
2764  if (bParallel) {
2765  silent_cerr("\"modify residual test\" "
2766  "not supported by schur data manager "
2767  "at line " << HP.GetLineData()
2768  << "; ignored" << std::endl);
2769  } else {
2770  bModResTest = true;
2772  "Modify residual test" << std::endl);
2773  }
2774  break;
2775 
2776  case ENFORCE_CONSTRAINT_EQUATIONS:
2777  if (HP.IsKeyWord("constraint" "violations")) {
2779 
2780  bSetScaleAlgebraic = !HP.IsKeyWord("scale" "factor");
2781 
2782  if (!bSetScaleAlgebraic) {
2783  dScaleAlgebraic = HP.GetReal();
2784  }
2785  } else if (HP.IsKeyWord("constraint" "violation" "rates")) {
2787  } else {
2788  silent_cerr("Keyword \"constraint violations\" or "
2789  "\"constraint violation rates\" expected at line "
2790  << HP.GetLineData() << std::endl);
2791 
2793  }
2794  break;
2795 
2796  case DERIVATIVESMAXITERATIONS: {
2797  iDerivativesMaxIterations = HP.GetInt();
2798  if (iDerivativesMaxIterations < 1) {
2799  iDerivativesMaxIterations = ::iDefaultMaxIterations;
2800  silent_cerr("warning, derivatives "
2801  "max iterations < 1 is illegal; "
2802  "using default value "
2803  << iDerivativesMaxIterations
2804  << std::endl);
2805  }
2806  DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives "
2807  "max iterations = "
2808  << iDerivativesMaxIterations
2809  << std::endl);
2810  break;
2811  }
2812 
2813  case FICTITIOUSSTEPSMAXITERATIONS:
2814  case DUMMYSTEPSMAXITERATIONS: {
2815  iDummyStepsMaxIterations = HP.GetInt();
2816  if (iDummyStepsMaxIterations < 1) {
2817  iDummyStepsMaxIterations = ::iDefaultMaxIterations;
2818  silent_cerr("warning, dummy steps "
2819  "max iterations < 1 is illegal; "
2820  "using default value "
2821  << iDummyStepsMaxIterations
2822  << std::endl);
2823  }
2824  DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps "
2825  "max iterations = "
2826  << iDummyStepsMaxIterations
2827  << std::endl);
2828  break;
2829  }
2830 
2831  case DERIVATIVESCOEFFICIENT: {
2832  const bool bAutoDerivCoef = HP.IsKeyWord("auto");
2833 
2834  if (!bAutoDerivCoef) {
2835  dDerivativesCoef = HP.GetReal();
2836  if (dDerivativesCoef <= 0.) {
2838  silent_cerr("warning, derivatives "
2839  "coefficient <= 0. is illegal; "
2840  "using default value "
2841  << dDerivativesCoef
2842  << std::endl);
2843  }
2844  DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives coefficient = "
2845  << dDerivativesCoef << std::endl);
2846  }
2847 
2848  if (bAutoDerivCoef || HP.IsKeyWord("auto")) {
2849  if (HP.IsKeyWord("max" "iterations")) {
2850  iDerivativesCoefMaxIter = HP.GetInt();
2851  if (iDerivativesCoefMaxIter < 0) {
2852  silent_cerr("number of iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
2854  }
2855 
2856  } else {
2857  iDerivativesCoefMaxIter = 6;
2858  }
2859 
2860  if (HP.IsKeyWord("factor")) {
2861  dDerivativesCoefFactor = HP.GetReal();
2862 
2863  if (dDerivativesCoefFactor <= 1.) {
2864  silent_cerr("factor must be greater than one at line " << HP.GetLineData() << std::endl);
2866  }
2867  }
2868  }
2869 
2870  break;
2871  }
2872 
2873  case NEWTONRAPHSON: {
2874  pedantic_cout("Newton Raphson is deprecated; use "
2875  "\"nonlinear solver: newton raphson "
2876  "[ , modified, <steps> ]\" instead"
2877  << std::endl);
2878  KeyWords NewRaph(KeyWords(HP.GetWord()));
2879  switch(NewRaph) {
2880  case MODIFIED:
2881  bTrueNewtonRaphson = 0;
2882  if (HP.IsArg()) {
2884  } else {
2886  }
2887  DEBUGLCOUT(MYDEBUG_INPUT, "Modified "
2888  "Newton-Raphson will be used; "
2889  "matrix will be assembled "
2890  "at most after "
2892  << " iterations" << std::endl);
2893  break;
2894 
2895  default:
2896  silent_cerr("warning: unknown case; "
2897  "using default" << std::endl);
2898 
2899  /* no break: fall-thru to next case */
2900  case NR_TRUE:
2901  bTrueNewtonRaphson = 1;
2903  break;
2904  }
2905  break;
2906  }
2907 
2908  case END:
2909  switch (KeyWords(HP.GetWord())) {
2910  case MULTISTEP:
2911  pedantic_cout("\"end: multistep;\" is deprecated; "
2912  "use \"end: initial value;\" instead." << std::endl);
2913  case INITIAL_VALUE:
2914  break;
2915 
2916  default:
2917  silent_cerr("\"end: initial value;\" expected "
2918  "at line " << HP.GetLineData()
2919  << "; aborting..." << std::endl);
2921  }
2922  goto EndOfCycle;
2923 
2924  case STRATEGY:
2925  pTSC = ReadTimeStepData(this, HP);
2926  break;
2927 
2928  case POD:
2929  silent_cerr("line " << HP.GetLineData()
2930  << ": POD analysis not supported (ignored)"
2931  << std::endl);
2932  for (; HP.IsArg();) {
2933  (void)HP.GetReal();
2934  }
2935  break;
2936 
2937  case EIGENANALYSIS:
2938  // initialize output precision: 0 means use default precision
2939  EigAn.iMatrixPrecision = 0;
2941 
2942  // read eigenanalysis time (to be changed)
2943  if (HP.IsKeyWord("list")) {
2944  int iNumTimes = HP.GetInt();
2945  if (iNumTimes <= 0) {
2946  silent_cerr("invalid number of eigenanalysis times "
2947  "at line " << HP.GetLineData()
2948  << std::endl);
2950  }
2951 
2952  EigAn.Analyses.resize(iNumTimes);
2953  for (std::vector<doublereal>::iterator i = EigAn.Analyses.begin();
2954  i != EigAn.Analyses.end(); ++i)
2955  {
2956  *i = HP.GetReal();
2957  if (i > EigAn.Analyses.begin() && *i <= *(i-1)) {
2958  silent_cerr("eigenanalysis times must be in strict ascending order "
2959  "at line " << HP.GetLineData()
2960  << std::endl);
2962  }
2963  }
2964 
2965  } else {
2966  EigAn.Analyses.resize(1);
2967  EigAn.Analyses[0] = HP.GetReal();
2968  }
2969 
2970  ASSERT(EigAn.Analyses.size() > 0);
2971  // initialize EigAn
2972  EigAn.currAnalysis = EigAn.Analyses.begin();
2973  EigAn.bAnalysis = true;
2974 
2975  // permute is the default; use "balance, no" to disable
2977 
2978  while (HP.IsArg()) {
2979  if (HP.IsKeyWord("parameter")) {
2980  EigAn.dParam = HP.GetReal();
2981 
2982  } else if (HP.IsKeyWord("output" "matrices")) {
2984 
2985  } else if (HP.IsKeyWord("output" "full" "matrices")) {
2986 #ifndef USE_EIG
2987  silent_cerr("\"output full matrices\" needs eigenanalysis support" << std::endl);
2989 #endif // !USE_EIG
2991 
2992  } else if (HP.IsKeyWord("output" "sparse" "matrices")) {
2994 
2995  } else if (HP.IsKeyWord("output" "eigenvectors")) {
2996 #ifndef USE_EIG
2997  silent_cerr("\"output eigenvectors\" needs eigenanalysis support" << std::endl);
2999 #endif // !USE_EIG
3001 
3002  } else if (HP.IsKeyWord("output" "geometry")) {
3004 
3005  } else if (HP.IsKeyWord("matrix" "output" "precision")) {
3006  EigAn.iMatrixPrecision = HP.GetInt();
3007  if (EigAn.iMatrixPrecision <= 0) {
3008  silent_cerr("negative or null \"matrix output precision\" "
3009  "parameter at line " << HP.GetLineData()
3010  << std::endl);
3012  }
3013 
3014  } else if (HP.IsKeyWord("results" "output" "precision")) {
3016  if (EigAn.iResultsPrecision <= 0) {
3017  silent_cerr("negative or null \"results output precision\" "
3018  "parameter at line " << HP.GetLineData()
3019  << std::endl);
3021  }
3022 
3023  } else if (HP.IsKeyWord("upper" "frequency" "limit")) {
3024  EigAn.dUpperFreq = HP.GetReal();
3025  if (EigAn.dUpperFreq < 0.) {
3026  silent_cerr("invalid \"upper frequency limit\" "
3027  "at line " << HP.GetLineData()
3028  << std::endl);
3030  }
3031 
3032  } else if (HP.IsKeyWord("lower" "frequency" "limit")) {
3033  EigAn.dLowerFreq = HP.GetReal();
3034  if (EigAn.dLowerFreq < 0.) {
3035  silent_cerr("invalid \"lower frequency limit\" "
3036  "at line " << HP.GetLineData()
3037  << std::endl);
3039  }
3040 
3041  } else if (HP.IsKeyWord("use" "lapack")) {
3043  silent_cerr("eigenanalysis routine already selected "
3044  "at line " << HP.GetLineData()
3045  << std::endl);
3047  }
3048 #ifdef USE_LAPACK
3050 #else // !USE_LAPACK
3051  silent_cerr("\"use lapack\" "
3052  "needs to configure --with-lapack "
3053  "at line " << HP.GetLineData()
3054  << std::endl);
3056 #endif // !USE_LAPACK
3057 
3058  } else if (HP.IsKeyWord("use" "arpack")) {
3060  silent_cerr("eigenanalysis routine already selected "
3061  "at line " << HP.GetLineData()
3062  << std::endl);
3064  }
3065 #ifdef USE_ARPACK
3067 
3068  EigAn.arpack.iNEV = HP.GetInt();
3069  if (EigAn.arpack.iNEV <= 0) {
3070  silent_cerr("invalid number of eigenvalues "
3071  "at line " << HP.GetLineData()
3072  << std::endl);
3074  }
3075 
3076  EigAn.arpack.iNCV = HP.GetInt();
3077  if (EigAn.arpack.iNCV <= 0
3078  || EigAn.arpack.iNCV <= EigAn.arpack.iNEV + 2)
3079  {
3080  silent_cerr("invalid number of Arnoldi vectors "
3081  "(must be > NEV+2) "
3082  "at line " << HP.GetLineData()
3083  << std::endl);
3085  }
3086 
3087  if (EigAn.arpack.iNCV <= 2*EigAn.arpack.iNEV) {
3088  silent_cerr("warning, possibly incorrect number of Arnoldi vectors "
3089  "(should be > 2*NEV) "
3090  "at line " << HP.GetLineData()
3091  << std::endl);
3092  }
3093 
3094  EigAn.arpack.dTOL = HP.GetReal();
3095  if (EigAn.arpack.dTOL < 0.) {
3096  silent_cerr("tolerance must be non-negative "
3097  "at line " << HP.GetLineData()
3098  << std::endl);
3099  EigAn.arpack.dTOL = 0.;
3100  }
3101 #else // !USE_ARPACK
3102  silent_cerr("\"use arpack\" "
3103  "needs to configure --with-arpack "
3104  "at line " << HP.GetLineData()
3105  << std::endl);
3107 #endif // !USE_ARPACK
3108 
3109  } else if (HP.IsKeyWord("use" "jdqz")) {
3111  silent_cerr("eigenanalysis routine already selected "
3112  "at line " << HP.GetLineData()
3113  << std::endl);
3115  }
3116 #ifdef USE_JDQZ
3118 
3119  EigAn.jdqz.kmax = HP.GetInt();
3120  if (EigAn.jdqz.kmax <= 0) {
3121  silent_cerr("invalid number of eigenvalues "
3122  "at line " << HP.GetLineData()
3123  << std::endl);
3125  }
3126 
3127  EigAn.jdqz.jmax = HP.GetInt();
3128  if (EigAn.jdqz.jmax < 20
3129  || EigAn.jdqz.jmax < 2*EigAn.jdqz.kmax)
3130  {
3131  silent_cerr("invalid size of the search space "
3132  "(must be >= 20 && >= 2*kmax) "
3133  "at line " << HP.GetLineData()
3134  << std::endl);
3136  }
3137 
3138  EigAn.jdqz.jmin = 2*EigAn.jdqz.kmax;
3139 
3140  EigAn.jdqz.eps = HP.GetReal();
3141  if (EigAn.jdqz.eps <= 0.) {
3142  silent_cerr("tolerance must be non-negative "
3143  "at line " << HP.GetLineData()
3144  << std::endl);
3145  EigAn.jdqz.eps = std::numeric_limits<doublereal>::epsilon();
3146  }
3147 #else // !USE_JDQZ
3148  silent_cerr("\"use jdqz\" "
3149  "needs to configure --with-jdqz "
3150  "at line " << HP.GetLineData()
3151  << std::endl);
3153 #endif // !USE_JDQZ
3154 
3155  } else if (HP.IsKeyWord("balance")) {
3156  if (HP.IsKeyWord("no")) {
3157  EigAn.uFlags &= ~EigenAnalysis::EIG_BALANCE;
3158 
3159  } else if (HP.IsKeyWord("permute")) {
3161 
3162  } else if (HP.IsKeyWord("scale")) {
3164 
3165  } else if (HP.IsKeyWord("all")) {
3167 
3168  } else {
3169  silent_cerr("unknown balance option "
3170  "at line " << HP.GetLineData()
3171  << std::endl);
3173  }
3174 
3175  } else if (HP.IsKeyWord("suffix" "width")) {
3176  if (HP.IsKeyWord("compute")) {
3178 
3179  } else {
3181  }
3182 
3183  } else if (HP.IsKeyWord("suffix" "format")) {
3185  // check?
3186 
3187  } else {
3188  silent_cerr("unknown option "
3189  "at line " << HP.GetLineData()
3190  << std::endl);
3192  }
3193  }
3194 
3195  // lower must be less than upper
3196  if (EigAn.dLowerFreq > EigAn.dUpperFreq) {
3197  silent_cerr("upper frequency limit " << EigAn.dUpperFreq
3198  << " less than lower frequency limit " << EigAn.dLowerFreq
3199  << std::endl);
3201  }
3202 
3203  // if only upper is defined, make lower equal to -upper
3204  if (EigAn.dLowerFreq == -1.) {
3206  }
3207 
3211  silent_cerr("sparse matrices output "
3212  "incompatible with lapack "
3213  "at line " << HP.GetLineData()
3214  << std::endl);
3216  }
3219  }
3220  break;
3221 
3224  silent_cerr("full matrices output "
3225  "incompatible with arpack "
3226  "at line " << HP.GetLineData()
3227  << std::endl);
3229  }
3232  }
3233  break;
3234 
3237  silent_cerr("full matrices output "
3238  "incompatible with jdqz "
3239  "at line " << HP.GetLineData()
3240  << std::endl);
3242  }
3245  }
3246  break;
3247 
3248  default:
3249  break;
3250  }
3251 
3252 #ifdef USE_EIG
3253  // if an eigenanalysis routine is selected
3254  // or sparse matrix output is not requested,
3255  // force direct eigensolution
3258  {
3260  }
3261 
3262  // if no eigenanalysis routine is selected,
3263  // force the use of LAPACK's
3266  {
3270  }
3271  }
3272 #else // !USE_EIG
3275  }
3276  silent_cerr("warning: \"eigenanalysis\" not supported; ignored" << std::endl);
3277 #endif // !USE_EIG
3278  break;
3279 
3280  case SOLVER:
3281  silent_cerr("\"solver\" keyword at line "
3282  << HP.GetLineData()
3283  << " is deprecated; "
3284  "use \"linear solver\" instead"
3285  << std::endl);
3286  case LINEARSOLVER:
3288  break;
3289 
3290  case INTERFACESOLVER:
3291  silent_cerr("\"interface solver\" keyword at line "
3292  << HP.GetLineData()
3293  << " is deprecated; "
3294  "use \"interface linear solver\" "
3295  "instead" << std::endl);
3296  case INTERFACELINEARSOLVER:
3297  ReadLinSol(CurrIntSolver, HP, true);
3298 
3299 #ifndef USE_MPI
3300  silent_cerr("Interface solver only allowed "
3301  "when compiled with MPI support" << std::endl);
3302 #endif /* ! USE_MPI */
3303  break;
3304 
3305  case NONLINEARSOLVER:
3306  switch (KeyWords(HP.GetWord())) {
3307  case DEFAULT:
3309  break;
3310 
3311  case NEWTONRAPHSON:
3313  break;
3314 
3315  case LINESEARCH:
3317  break;
3318 
3319  case MATRIXFREE:
3321  break;
3322 
3323  default:
3324  silent_cerr("unknown nonlinear solver "
3325  "at line " << HP.GetLineData()
3326  << std::endl);
3328  }
3329 
3330  switch (NonlinearSolverType) {
3333  bTrueNewtonRaphson = true;
3334  bKeepJac = false;
3336 
3338  break;
3339  }
3340 
3341  if (HP.IsKeyWord("modified")) {
3342  bTrueNewtonRaphson = false;
3344 
3345  if (HP.IsKeyWord("keep" "jacobian")) {
3346  pedantic_cout("Use of deprecated \"keep jacobian\" "
3347  "at line " << HP.GetLineData() << std::endl);
3348  bKeepJac = true;
3349 
3350  } else if (HP.IsKeyWord("keep" "jacobian" "matrix")) {
3351  bKeepJac = true;
3352  }
3353 
3354  DEBUGLCOUT(MYDEBUG_INPUT, "modified "
3355  "Newton-Raphson "
3356  "will be used; "
3357  "matrix will be "
3358  "assembled at most "
3359  "after "
3361  << " iterations"
3362  << std::endl);
3363  if (HP.IsKeyWord("honor" "element" "requests")) {
3364  bHonorJacRequest = true;
3366  "honor elements' "
3367  "request to update "
3368  "the preconditioner"
3369  << std::endl);
3370  }
3371  }
3372 
3374  while (HP.IsArg()) {
3375  if (HP.IsKeyWord("tolerance" "x")) {
3376  LineSearch.dTolX = HP.GetReal();
3377  if (LineSearch.dTolX < 0.) {
3378  silent_cerr("tolerance x must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3380  }
3381  } else if (HP.IsKeyWord("tolerance" "min")) {
3382  LineSearch.dTolMin = HP.GetReal();
3383  if (LineSearch.dTolMin < 0.) {
3384  silent_cerr("tolerance min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3386  }
3387  } else if (HP.IsKeyWord("max" "iterations")) {
3389  if (LineSearch.iMaxIterations < 0) {
3390  silent_cerr("max iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3392  }
3393  } else if (HP.IsKeyWord("alpha")) {
3394  LineSearch.dAlpha = HP.GetReal();
3395  if (LineSearch.dAlpha < 0.) {
3396  silent_cerr("alpha must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3398  }
3399  } else if (HP.IsKeyWord("lambda" "min")) {
3400  LineSearch.dLambdaMin = HP.GetReal();
3401  if (LineSearch.dLambdaMin < 0.) {
3402  silent_cerr("lambda min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
3404  }
3405  if (HP.IsKeyWord("relative")) {
3406  if (HP.GetYesNoOrBool()) {
3408  } else {
3409  LineSearch.uFlags &= ~LineSearchParameters::RELATIVE_LAMBDA_MIN;
3410  }
3411  }
3412  } else if (HP.IsKeyWord("lambda" "factor" "min")) {
3414  if (LineSearch.dLambdaFactMin <= 0. || LineSearch.dLambdaFactMin >= 1.) {
3415  silent_cerr("lambda factor min must be in between zero and one at line" << HP.GetLineData() << std::endl);
3417  }
3418  } else if (HP.IsKeyWord("max" "step")) {
3419  LineSearch.dMaxStep = HP.GetReal();
3420  if (LineSearch.dMaxStep <= 0.) {
3421  silent_cerr("max step must be greater than zero at line " << HP.GetLineData() << std::endl);
3423  }
3424  } else if (HP.IsKeyWord("zero" "gradient")) {
3425  if (HP.IsKeyWord("continue")) {
3426  if (HP.GetYesNoOrBool()) {
3428  } else {
3429  LineSearch.uFlags &= ~LineSearchParameters::ZERO_GRADIENT_CONTINUE;
3430  }
3431  }
3432  else {
3433  silent_cerr("Keyword \"continue\" expected at line " << HP.GetLineData() << std::endl);
3435  }
3436  } else if (HP.IsKeyWord("divergence" "check")) {
3437  if (HP.GetYesNoOrBool()) {
3439  } else {
3440  LineSearch.uFlags &= ~LineSearchParameters::DIVERGENCE_CHECK;
3441  }
3442  if (HP.IsKeyWord("factor")) {
3444  if (LineSearch.dDivergenceCheck <= 0.) {
3445  silent_cerr("divergence check factor must be greater than zero at line " << HP.GetLineData() << std::endl);
3447  }
3448  }
3449  } else if (HP.IsKeyWord("algorithm")) {
3450  LineSearch.uFlags &= ~LineSearchParameters::ALGORITHM;
3451  if (HP.IsKeyWord("cubic")) {
3453  } else if (HP.IsKeyWord("factor")) {
3455  } else {
3456  silent_cerr("Keyword \"cubic\" or \"factor\" expected at line " << HP.GetLineData() << std::endl);
3458  }
3459  } else if (HP.IsKeyWord("scale" "newton" "step")) {
3460  if (HP.GetYesNoOrBool()) {
3462  } else {
3463  LineSearch.uFlags &= ~LineSearchParameters::SCALE_NEWTON_STEP;
3464  }
3465  if (HP.IsKeyWord("min" "scale")) {
3467  if (LineSearch.dMinStepScale < 0. || LineSearch.dMinStepScale > 1.) {
3468  silent_cerr("min scale must be in range [0-1] at line " << HP.GetLineData() << std::endl);
3470  }
3471  }
3472  } else if (HP.IsKeyWord("print" "convergence" "info")) {
3473  if (HP.GetYesNoOrBool()) {
3475  } else {
3476  LineSearch.uFlags &= ~LineSearchParameters::PRINT_CONVERGENCE_INFO;
3477  }
3478  } else if (HP.IsKeyWord("verbose")) {
3479  if (HP.GetYesNoOrBool()) {
3481  } else {
3482  LineSearch.uFlags &= ~LineSearchParameters::VERBOSE_MODE;
3483  }
3484  } else if (HP.IsKeyWord("abort" "at" "lambda" "min")) {
3485  if (HP.GetYesNoOrBool()) {
3487  } else {
3488  LineSearch.uFlags &= ~LineSearchParameters::ABORT_AT_LAMBDA_MIN;
3489  }
3490  } else if (HP.IsKeyWord("non" "negative" "slope" "continue")) {
3491  if (HP.GetYesNoOrBool()) {
3493  } else {
3494  LineSearch.uFlags &= ~LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE;
3495  }
3496  } else {
3497  silent_cerr("Keyword \"tolerance x\", "
3498  "\"tolerance min\", \"max iterations\", \"alpha\", "
3499  "\"lambda min\" \"lambda factor min\", \"max step\" "
3500  "\"divergence check\", \"algorithm\", \"scale newton step\" "
3501  "\"print convergence info\", \"verbose\", "
3502  "\"abort at lambda min\" "
3503  "or \"zero gradient\" expected at line " << HP.GetLineData() << std::endl);
3505  }
3506  }
3507  }
3508  break;
3509 
3511  switch (KeyWords(HP.GetWord())) {
3512  case DEFAULT:
3514  break;
3515 
3516 
3517  case BICGSTAB:
3519  break;
3520 
3521  case GMRES:
3523  break;
3524 
3525  default:
3526  silent_cerr("unknown iterative "
3527  "solver at line "
3528  << HP.GetLineData()
3529  << std::endl);
3531  }
3532 
3533  if (HP.IsKeyWord("tolerance")) {
3534  dIterTol = HP.GetReal();
3535  DEBUGLCOUT(MYDEBUG_INPUT,"inner "
3536  "iterative solver "
3537  "tolerance: "
3538  << dIterTol
3539  << std::endl);
3540  }
3541 
3542  if (HP.IsKeyWord("steps")) {
3543  iIterativeMaxSteps = HP.GetInt();
3544  DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
3545  "number of inner "
3546  "steps for iterative "
3547  "solver: "
3549  << std::endl);
3550  }
3551 
3552  if (HP.IsKeyWord("tau")) {
3553  dIterertiveTau = HP.GetReal();
3555  "tau scaling "
3556  "coefficient "
3557  "for iterative "
3558  "solver: "
3559  << dIterertiveTau
3560  << std::endl);
3561  }
3562 
3563  if (HP.IsKeyWord("eta")) {
3564  dIterertiveEtaMax = HP.GetReal();
3565  DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
3566  "eta coefficient "
3567  "for iterative "
3568  "solver: "
3570  << std::endl);
3571  }
3572 
3573  if (HP.IsKeyWord("preconditioner")) {
3574  KeyWords KPrecond = KeyWords(HP.GetWord());
3575  switch (KPrecond) {
3576  case FULLJACOBIAN:
3577  case FULLJACOBIANMATRIX:
3579  if (HP.IsKeyWord("steps")) {
3580  iPrecondSteps = HP.GetInt();
3582  "number of steps "
3583  "before recomputing "
3584  "the preconditioner: "
3585  << iPrecondSteps
3586  << std::endl);
3587  }
3588  if (HP.IsKeyWord("honor" "element" "requests")) {
3589  bHonorJacRequest = true;
3591  "honor elements' "
3592  "request to update "
3593  "the preconditioner"
3594  << std::endl);
3595  }
3596  break;
3597 
3598  /* add other preconditioners
3599  * here */
3600 
3601  default:
3602  silent_cerr("unknown "
3603  "preconditioner "
3604  "at line "
3605  << HP.GetLineData()
3606  << std::endl);
3608  }
3609  break;
3610  }
3611  break;
3612 
3613  default:
3614  ASSERT(0);
3616  }
3617  break;
3618 
3619  case REALTIME:
3620  pRTSolver = ReadRTSolver(this, HP);
3621  break;
3622 
3623  case THREADS:
3624  if (HP.IsKeyWord("auto")) {
3625 #ifdef USE_MULTITHREAD
3626  int n = get_nprocs();
3627  /* sanity checks ... */
3628  if (n <= 0) {
3629  silent_cerr("got " << n << " CPUs "
3630  "at line "
3631  << HP.GetLineData()
3632  << std::endl);
3633  nThreads = 1;
3634  } else {
3635  nThreads = n;
3636  }
3637 #else /* ! USE_MULTITHREAD */
3638  silent_cerr("configure with "
3639  "--enable-multithread "
3640  "for multithreaded assembly"
3641  << std::endl);
3642 #endif /* ! USE_MULTITHREAD */
3643 
3644  } else if (HP.IsKeyWord("disable")) {
3645 #ifdef USE_MULTITHREAD
3646  nThreads = 1;
3647 #endif /* USE_MULTITHREAD */
3648 
3649  } else {
3650 #ifdef USE_MULTITHREAD
3651  bool bAssembly = false;
3652  bool bSolver = false;
3653  bool bAll = true;
3654  unsigned nt;
3655 #endif // USE_MULTITHREAD
3656 
3657  if (HP.IsKeyWord("assembly")) {
3658 #ifdef USE_MULTITHREAD
3659  bAll = false;
3660  bAssembly = true;
3661 #endif // USE_MULTITHREAD
3662 
3663  } else if (HP.IsKeyWord("solver")) {
3664 #ifdef USE_MULTITHREAD
3665  bAll = false;
3666  bSolver = true;
3667 #endif // USE_MULTITHREAD
3668  }
3669 
3670 #ifdef USE_MULTITHREAD
3671  nt =
3672 #endif // USE_MULTITHREAD
3673  HP.GetInt();
3674 
3675 #ifdef USE_MULTITHREAD
3676  if (bAll || bAssembly) {
3677  nThreads = nt;
3678  }
3679 
3680  if (bAll || bSolver) {
3681  bSolverThreads = true;
3682  nSolverThreads = nt;
3683  }
3684 #else /* ! USE_MULTITHREAD */
3685  silent_cerr("configure with "
3686  "--enable-multithread "
3687  "for multithreaded assembly"
3688  << std::endl);
3689 #endif /* ! USE_MULTITHREAD */
3690  }
3691  break;
3692 
3693 
3694  default:
3695  silent_cerr("unknown description at line "
3696  << HP.GetLineData() << "; aborting..."
3697  << std::endl);
3699  }
3700  }
3701 
3702 EndOfCycle: /* esce dal ciclo di lettura */
3703 
3704  if (pTSC == 0) {
3705  pedantic_cout("No time step control strategy defined; defaulting to NoChange time step control" );
3706  pTSC = new NoChange(this);
3707  }
3708 
3709  if (dFinalTime < dInitialTime) {
3711  }
3712 
3713  if (dFinalTime == dInitialTime) {
3715  }
3716 
3717  /* Metodo di integrazione di default */
3718  if (!bMethod) {
3720 
3721  /* FIXME: maybe we should use a better value
3722  * like 0.6; however, BDF should be conservative */
3724 
3725  /* DriveCaller per Rho asintotico per variabili algebriche */
3727 
3728  RegularType = INT_MS2;
3729  }
3730 
3731  /* Metodo di integrazione di default */
3732  if (iDummyStepsNumber && !bDummyStepsMethod) {
3734 
3736 
3737  /* DriveCaller per Rho asintotico per variabili algebriche */
3739 
3740  DummyType = INT_MS2;
3741  }
3742 
3743  /* costruzione dello step solver derivative */
3746  DerivativeSolver(dDerivativesTol,
3747  0.,
3749  iDerivativesMaxIterations,
3750  bModResTest,
3751  iDerivativesCoefMaxIter,
3752  dDerivativesCoefFactor));
3753 
3754  /* First step prediction must always be Crank-Nicolson for accuracy */
3755  if (iDummyStepsNumber) {
3758  CrankNicolsonIntegrator(dDummyStepsTolerance,
3759  0.,
3760  iDummyStepsMaxIterations,
3761  bModResTest));
3762 
3763  /* costruzione dello step solver dummy */
3764  switch (DummyType) {
3765  case INT_CRANKNICOLSON:
3768  CrankNicolsonIntegrator(dDummyStepsTolerance,
3769  0.,
3770  iDummyStepsMaxIterations,
3771  bModResTest));
3772  break;
3773 
3774  case INT_MS2:
3777  MultistepSolver(dDummyStepsTolerance,
3778  0.,
3779  iDummyStepsMaxIterations,
3780  pRhoDummy,
3782  bModResTest));
3783  break;
3784 
3785  case INT_HOPE:
3787  HopeSolver,
3788  HopeSolver(dDummyStepsTolerance,
3789  dSolutionTol,
3790  iDummyStepsMaxIterations,
3791  pRhoDummy,
3793  bModResTest));
3794  break;
3795 
3796  case INT_THIRDORDER:
3797  if (pRhoDummy == 0) {
3800  AdHocThirdOrderIntegrator(dDummyStepsTolerance,
3801  dSolutionTol,
3802  iDummyStepsMaxIterations,
3803  bModResTest));
3804  } else {
3807  TunableThirdOrderIntegrator(dDummyStepsTolerance,
3808  dSolutionTol,
3809  iDummyStepsMaxIterations,
3810  pRhoDummy,
3811  bModResTest));
3812  }
3813  break;
3814 
3815  case INT_IMPLICITEULER:
3818  ImplicitEulerIntegrator(dDummyStepsTolerance,
3819  dSolutionTol, iDummyStepsMaxIterations,
3820  bModResTest));
3821  break;
3822 
3823  default:
3824  silent_cerr("unknown dummy steps integration method"
3825  << std::endl);
3827  break;
3828  }
3829  }
3830 
3834  dSolutionTol,
3836  bModResTest));
3837 
3838  /* costruzione dello step solver per i passi normali */
3839  switch (RegularType) {
3840  case INT_CRANKNICOLSON:
3844  dSolutionTol,
3846  bModResTest));
3847  break;
3848 
3849  case INT_MS2:
3852  MultistepSolver(dTol,
3853  dSolutionTol,
3855  pRhoRegular,
3857  bModResTest));
3858  break;
3859 
3860  case INT_HOPE:
3862  HopeSolver,
3863  HopeSolver(dTol,
3864  dSolutionTol,
3866  pRhoRegular,
3868  bModResTest));
3869  break;
3870 
3871  case INT_THIRDORDER:
3872  if (pRhoRegular == 0) {
3876  dSolutionTol,
3878  bModResTest));
3879  } else {
3883  dSolutionTol,
3885  pRhoRegular,
3886  bModResTest));
3887  }
3888  break;
3889 
3890  case INT_IMPLICITEULER:
3894  dSolutionTol,
3896  bModResTest));
3897  break;
3898 
3899  default:
3900  silent_cerr("Unknown integration method" << std::endl);
3902  break;
3903  }
3904 
3905  if (bSetScaleAlgebraic) {
3907  }
3908 
3909 #ifdef USE_MULTITHREAD
3910  if (bSolverThreads) {
3911  if (CurrLinearSolver.SetNumThreads(nSolverThreads)) {
3912  silent_cerr("linear solver "
3914  << " does not support "
3915  "threaded solution" << std::endl);
3916  }
3917  }
3918 #endif /* USE_MULTITHREAD */
3919 }
3920 
3921 static int
3922 do_eig(const doublereal& b, const doublereal& re,
3923  const doublereal& im, const doublereal& h,
3924  doublereal& sigma, doublereal& omega,
3925  doublereal& csi, doublereal& freq)
3926 {
3927  // denominator
3928  doublereal d = re + b;
3929  d *= d;
3930  d += im*im;
3931  d *= h/2.;
3932 
3933  // real & imag
3934  sigma = (re*re - b*b + im*im)/d;
3935  omega = 2.*b*im/d;
3936 
3937  // frequency and damping factor
3938  if (im != 0.) {
3939  d = sigma*sigma + omega*omega;
3940  if (d > std::numeric_limits<doublereal>::epsilon()) {
3941  csi = -100*sigma/sqrt(d);
3942 
3943  } else {
3944  csi = 0.;
3945  }
3946 
3947  freq = omega/(2*M_PI);
3948 
3949  } else {
3950  if (std::abs(sigma) < std::numeric_limits<doublereal>::epsilon()) {
3951  csi = 0.;
3952 
3953  } else {
3954  csi = -100.*copysign(1, sigma);
3955  }
3956 
3957  freq = 0.;
3958  }
3959 
3960  return 0;
3961 }
3962 
3963 // Writes eigenvalues to the .out file in human-readable form
3964 static void
3966  const VectorHandler& R, const VectorHandler& I,
3967  const doublereal& dShiftR,
3968  DataManager* pDM,
3969  const Solver::EigenAnalysis *pEA,
3970  integer iLow, integer iHigh,
3971  std::vector<bool>& vOut)
3972 {
3973  std::ostream& Out = pDM->GetOutFile();
3974 
3975  /* Output? */
3976  Out << "Mode n. " " " " Real " " " " Imag " " " " " " Damp % " " Freq Hz" << std::endl;
3977 
3978  integer iNVec = R.iGetSize();
3979 
3980  for (int iCnt = 1; iCnt <= iNVec; iCnt++) {
3981  doublereal b = pBeta ? (*pBeta)(iCnt) : 1.;
3982  doublereal re = R(iCnt) + dShiftR;
3983  doublereal im = I(iCnt);
3984  doublereal sigma;
3985  doublereal omega;
3986  doublereal csi;
3987  doublereal freq;
3988 
3989  if (iCnt < iLow || iCnt > iHigh) {
3990  vOut[iCnt - 1] = false;
3991  continue;
3992  }
3993 
3994  const doublereal& h = pEA->dParam;
3995  do_eig(b, re, im, h, sigma, omega, csi, freq);
3996 
3997  if (freq < pEA->dLowerFreq || freq > pEA->dUpperFreq) {
3998  vOut[iCnt - 1] = false;
3999  continue;
4000  }
4001 
4002  vOut[iCnt - 1] = true;
4003 
4004  Out << std::setw(8) << iCnt << ": "
4005  << std::setw(12) << sigma << " + " << std::setw(12) << omega << " j";
4006 
4007  if (fabs(csi) > std::numeric_limits<doublereal>::epsilon()) {
4008  Out << " " << std::setw(12) << csi;
4009  } else {
4010  Out << " " << std::setw(12) << 0.;
4011  }
4012 
4013  Out << " " << std::setw(12) << freq;
4014 
4015  Out << std::endl;
4016  }
4017 }
4018 
4019 #ifdef USE_LAPACK
4020 // Computes eigenvalues and eigenvectors using LAPACK's
4021 // generalized non-symmetric eigenanalysis
4022 static void
4023 eig_lapack(const MatrixHandler* pMatA, const MatrixHandler* pMatB,
4025  bool bNewLine, const unsigned uCurr)
4026 {
4027  const FullMatrixHandler& MatA = dynamic_cast<const FullMatrixHandler &>(*pMatA);
4028  const FullMatrixHandler& MatB = dynamic_cast<const FullMatrixHandler &>(*pMatB);
4029 
4030  char sB[2] = "N";
4033  {
4034  sB[0] = 'B';
4035 
4036  } else if (pEA->uFlags & Solver::EigenAnalysis::EIG_PERMUTE) {
4037  sB[0] = 'P';
4038 
4039  } else if (pEA->uFlags & Solver::EigenAnalysis::EIG_SCALE) {
4040  sB[0] = 'S';
4041  }
4042 
4043  char sL[2] = "V";
4044  char sR[2] = "V";
4045  char sS[2] = "N";
4046 
4047  // iNumDof is a member, set after dataman constr.
4048  integer iSize = MatA.iGetNumRows();
4049 
4050  // Minimum workspace size. To be improved.
4051  // NOTE: optimal iWorkSize is computed by dggev() and dggevx()
4052  // when called with iWorkSize = -1.
4053  // The computed value is stored in WorkVec[0].
4054  integer iWorkSize = -1;
4055  integer iMinWorkSize = -1;
4056  integer iInfo = 0;
4057 
4058  integer iILO = 1, iIHI = iSize;
4059  doublereal dABNRM = -1., dBBNRM = -1.;
4060 
4061  doublereal dDmy;
4062  integer iDmy;
4063  logical lDmy;
4064  doublereal dWV;
4065  if (sB[0] == 'N') {
4066  iMinWorkSize = 8*iSize;
4067 
4068  __FC_DECL__(dggev)(
4069  sL, // JOBVL
4070  sR, // JOBVR
4071  &iSize, // N
4072  &dDmy, // A
4073  &iSize, // LDA
4074  &dDmy, // B
4075  &iSize, // LDB
4076  &dDmy, // ALPHAR
4077  &dDmy, // ALPHAI
4078  &dDmy, // BETA
4079  &dDmy, // VL
4080  &iSize, // LDVL
4081  &dDmy, // VR
4082  &iSize, // LDVR
4083  &dWV, // WORK
4084  &iWorkSize, // LWORK
4085  &iInfo);
4086 
4087  } else {
4088  iMinWorkSize = 6*iSize;
4089 
4090  __FC_DECL__(dggevx)(
4091  sB, // BALANC
4092  sL, // JOBVL
4093  sR, // JOBVR
4094  sS, // SENSE
4095  &iSize, // N
4096  &dDmy, // A
4097  &iSize, // LDA
4098  &dDmy, // B
4099  &iSize, // LDB
4100  &dDmy, // ALPHAR
4101  &dDmy, // ALPHAI
4102  &dDmy, // BETA
4103  &dDmy, // VL
4104  &iSize, // LDVL
4105  &dDmy, // VR
4106  &iSize, // LDVR
4107  &iILO, // ILO
4108  &iIHI, // IHI
4109  &dDmy, // LSCALE
4110  &dDmy, // RSCALE
4111  &dABNRM, // ABNRM
4112  &dBBNRM, // BBNRM
4113  &dDmy, // RCONDE
4114  &dDmy, // RCONDV
4115  &dWV, // WORK
4116  &iWorkSize, // LWORK
4117  &iDmy, // IWORK
4118  &lDmy, // BWORK
4119  &iInfo);
4120  }
4121 
4122  if (iInfo != 0) {
4123  if (bNewLine && silent_err) {
4124  silent_cerr(std::endl);
4125  bNewLine = false;
4126  }
4127  silent_cerr("dggev[x]() query for worksize failed "
4128  "INFO=" << iInfo << std::endl);
4129  iInfo = 0;
4130  }
4131 
4132  iWorkSize = (integer)dWV;
4133  if (iWorkSize < iMinWorkSize) {
4134  if (bNewLine && silent_err) {
4135  silent_cerr(std::endl);
4136  bNewLine = false;
4137  }
4138  silent_cerr("dggev[x]() asked for a worksize " << iWorkSize
4139  << " less than the minimum, " << iMinWorkSize
4140  << "; using the minimum" << std::endl);
4141  iWorkSize = iMinWorkSize;
4142  }
4143 
4144  // Workspaces
4145  // 2 matrices iSize x iSize
4146  // 5 vectors iSize x 1
4147  // 1 vector iWorkSize x 1
4148  doublereal* pd = 0;
4149  int iTmpSize = 2*(iSize*iSize) + 3*iSize + iWorkSize;
4150  if (sB[0] != 'N') {
4151  iTmpSize += 2*iSize;
4152  }
4153  SAFENEWARR(pd, doublereal, iTmpSize);
4154 #if defined HAVE_MEMSET
4155  memset(pd, 0, iTmpSize*sizeof(doublereal));
4156 #else // !HAVE_MEMSET
4157  for (int iCnt = iTmpSize; iCnt-- > 0; ) {
4158  pd[iCnt] = 0.;
4159  }
4160 #endif // !HAVE_MEMSET
4161 
4162  // 2 pointer arrays iSize x 1 for the matrices
4163  doublereal** ppd = 0;
4164  SAFENEWARR(ppd, doublereal*, 2*iSize);
4165 
4166  // Data Handlers
4167  doublereal* pdTmp = pd;
4168  doublereal** ppdTmp = ppd;
4169 
4170  FullMatrixHandler MatVL(pdTmp, ppdTmp, iSize*iSize, iSize, iSize);
4171  pdTmp += iSize*iSize;
4172  ppdTmp += iSize;
4173 
4174  FullMatrixHandler MatVR(pdTmp, ppdTmp, iSize*iSize, iSize, iSize);
4175  pdTmp += iSize*iSize;
4176 
4177  MyVectorHandler AlphaR(iSize, pdTmp);
4178  pdTmp += iSize;
4179 
4180  MyVectorHandler AlphaI(iSize, pdTmp);
4181  pdTmp += iSize;
4182 
4183  MyVectorHandler Beta(iSize, pdTmp);
4184  pdTmp += iSize;
4185 
4186  MyVectorHandler LScale;
4187  MyVectorHandler RScale;
4188  if (sB[0] != 'N') {
4189  LScale.Attach(iSize, pdTmp);
4190  pdTmp += iSize;
4191 
4192  RScale.Attach(iSize, pdTmp);
4193  pdTmp += iSize;
4194  }
4195 
4196  MyVectorHandler WorkVec(iWorkSize, pdTmp);
4197 
4198  // Eigenanalysis
4199  // NOTE: according to lapack's documentation, dgegv() is deprecated
4200  // in favour of dggev()... I find dgegv() a little bit faster (10%?)
4201  // for typical problems (N ~ 1000).
4202  if (sB[0] == 'N') {
4203  __FC_DECL__(dggev)(
4204  sL, // JOBVL
4205  sR, // JOBVR
4206  &iSize, // N
4207  const_cast<doublereal *>(MatA.pdGetMat()), // A; remove const'ness (hack)
4208  &iSize, // LDA
4209  const_cast<doublereal *>(MatB.pdGetMat()), // B; remove const'ness (hack)
4210  &iSize, // LDB
4211  AlphaR.pdGetVec(), // ALPHAR
4212  AlphaI.pdGetVec(), // ALPHAI
4213  Beta.pdGetVec(), // BETA
4214  MatVL.pdGetMat(), // VL
4215  &iSize, // LDVL
4216  MatVR.pdGetMat(), // VR
4217  &iSize, // LDVR
4218  WorkVec.pdGetVec(), // WORK
4219  &iWorkSize, // LWORK
4220  &iInfo); // INFO
4221 
4222  } else {
4223  std::vector<integer> iIWORK(iSize + 6);
4224 
4225  __FC_DECL__(dggevx)(
4226  sB, // BALANCE
4227  sL, // JOBVL
4228  sR, // JOBVR
4229  sS, // SENSE
4230  &iSize, // N
4231  const_cast<doublereal *>(MatA.pdGetMat()), // A; remove const'ness (hack)
4232  &iSize, // LDA
4233  const_cast<doublereal *>(MatB.pdGetMat()), // B; remove const'ness (hack)
4234  &iSize, // LDB
4235  AlphaR.pdGetVec(), // ALPHAR
4236  AlphaI.pdGetVec(), // ALPHAI
4237  Beta.pdGetVec(), // BETA
4238  MatVL.pdGetMat(), // VL
4239  &iSize, // LDVL
4240  MatVR.pdGetMat(), // VR
4241  &iSize, // LDVR
4242  &iILO, // ILO
4243  &iIHI, // IHI
4244  LScale.pdGetVec(), // LSCALE
4245  RScale.pdGetVec(), // RSCALE
4246  &dABNRM, // ABNRM
4247  &dBBNRM, // BBNRM
4248  &dDmy, // RCONDE
4249  &dDmy, // RCONDV
4250  WorkVec.pdGetVec(), // WORK
4251  &iWorkSize, // LWORK
4252  &iIWORK[0], // IWORK
4253  &lDmy, // BWORK
4254  &iInfo); // INFO
4255  }
4256 
4257  std::ostream& Out = pDM->GetOutFile();
4258  Out << "Info: " << iInfo << ", ";
4259 
4260  if (iInfo == 0) {
4261  // = 0: successful exit
4262  Out << "success"
4263  << " BALANC=\"" << sB << "\""
4264  << " ILO=" << iILO
4265  << " IHI=" << iIHI
4266  << " ABNRM=" << dABNRM
4267  << " BBNRM=" << dBBNRM
4268  << std::endl;
4269 
4270  } else if (iInfo < 0) {
4271  const char *arg[] = {
4272  "JOBVL",
4273  "JOBVR",
4274  "N",
4275  "A",
4276  "LDA",
4277  "B",
4278  "LDB",
4279  "ALPHAR",
4280  "ALPHAI",
4281  "BETA",
4282  "VL",
4283  "LDVL",
4284  "VR",
4285  "LDVR",
4286  "WORK",
4287  "LWORK",
4288  "INFO",
4289  NULL
4290  };
4291 
4292  const char *argx[] = {
4293  "BALANCE",
4294  "JOBVL",
4295  "JOBVR",
4296  "SENSE",
4297  "N",
4298  "A",
4299  "LDA",
4300  "B",
4301  "LDB",
4302  "ALPHAR",
4303  "ALPHAI",
4304  "BETA",
4305  "VL",
4306  "LDVL",
4307  "VR",
4308  "LDVR",
4309  "ILO",
4310  "IHI",
4311  "LSCALE",
4312  "RSCALE",
4313  "ABNRM",
4314  "BBNRM",
4315  "RCONDE",
4316  "RCONDV",
4317  "WORK",
4318  "LWORK",
4319  "IWORK",
4320  "BWORK",
4321  "INFO",
4322  NULL
4323  };
4324 
4325  const char **argv = (sB[0] == 'N' ? arg : argx );
4326 
4327  // < 0: if INFO = -i, the i-th argument had an illegal value.
4328  Out << "argument #" << -iInfo
4329  << " (" << argv[-iInfo - 1] << ") "
4330  << "was passed an illegal value" << std::endl;
4331 
4332  } else if (iInfo > 0 && iInfo <= iSize) {
4333  /* = 1,...,N:
4334  * The QZ iteration failed. No eigenvectors have been
4335  * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
4336  * should be correct for j=INFO+1,...,N. */
4337  Out << "the QZ iteration failed, but eigenvalues "
4338  << iInfo + 1 << "->" << iSize << " should be correct"
4339  << std::endl;
4340 
4341  } else if (iInfo > iSize) {
4342  const char* const sErrs[] = {
4343  "DHGEQZ (other than QZ iteration failed in DHGEQZ)",
4344  "DTGEVC",
4345  NULL
4346  };
4347 
4348  Out << "error return from " << sErrs[iInfo - iSize - 1]
4349  << std::endl;
4350  }
4351 
4352  std::vector<bool> vOut(iSize);
4353  output_eigenvalues(&Beta, AlphaR, AlphaI, 0., pDM, pEA, iILO, iIHI, vOut);
4354 
4356  pDM->OutputEigGeometry(uCurr, pEA->iResultsPrecision);
4357  }
4358 
4360  pDM->OutputEigenvectors(&Beta, AlphaR, AlphaI, 0.,
4361  &MatVL, MatVR, vOut, uCurr, pEA->iResultsPrecision);
4362  }
4363 
4364  SAFEDELETEARR(pd);
4365  SAFEDELETEARR(ppd);
4366 }
4367 #endif // USE_LAPACK
4368 
4369 #ifdef USE_ARPACK
4370 // Computes eigenvalues and eigenvectors using ARPACK's
4371 // canonical non-symmetric eigenanalysis
4372 static void
4373 eig_arpack(const MatrixHandler* pMatA, SolutionManager* pSM,
4375  bool bNewLine, const unsigned uCurr)
4376 {
4378  = dynamic_cast<NaiveSparsePermSolutionManager<Colamd_ordering> &>(*pSM);
4379  const NaiveMatrixHandler& MatA = dynamic_cast<const NaiveMatrixHandler &>(*pMatA);
4380 
4381  // shift
4382  doublereal SIGMAR = 0.;
4383  doublereal SIGMAI = 0.;
4384 
4385  // arpack-related vars
4386  integer IDO; // 0 at first iteration; then set by dnaupd
4387  const char *BMAT; // 'I' for standard problem
4388  integer N; // size of problem
4389  const char *WHICH; // "SM" to request smallest eigenvalues
4390  integer NEV; // number of eigenvalues
4391  doublereal TOL; // -1 to use machine precision
4392  std::vector<doublereal> RESID; // residual vector (ignored if IDO==0)
4393  integer NCV; // number of vectors in subspace
4394  std::vector<doublereal> V; // Schur basis
4395  integer LDV; // leading dimension of V (==N!)
4396  integer IPARAM[11] = { 0 };
4397  integer IPNTR[14] = { 0 };
4398  std::vector<doublereal> WORKD;
4399  std::vector<doublereal> WORKL;
4400  integer LWORKL;
4401  integer INFO;
4402 
4403  IDO = 0;
4404  BMAT = "I";
4405  N = MatA.iGetNumRows();
4406  WHICH = "SM";
4407  NEV = pEA->arpack.iNEV;
4408  if (NEV > N) {
4409  silent_cerr("eig_arpack: invalid NEV=" << NEV << " > size of problem (=" << N << ")" << std::endl);
4411  }
4412 
4413  TOL = pEA->arpack.dTOL;
4414  RESID.resize(N, 0.);
4415  NCV = pEA->arpack.iNCV;
4416  if (NCV > N) {
4417  silent_cerr("eig_arpack: invalid NCV=" << NCV << " > size of problem (=" << N << ")" << std::endl);
4419  }
4420  // NOTE: we accommodate NCV + 1 vectors to handle the case
4421  // of a complex last eigenvalue which requires a vector
4422  // for the real part and one for the imaginary part
4423  V.resize(N*(NCV + 1), 0.);
4424  LDV = N;
4425  IPARAM[0] = 1;
4426  IPARAM[2] = 300; // configurable?
4427  IPARAM[3] = 1;
4428  IPARAM[6] = 1; // mode 1: canonical problem
4429  WORKD.resize(3*N, 0.);
4430  LWORKL = 3*NCV*NCV + 6*NCV;
4431  WORKL.resize(LWORKL, 0.);
4432  INFO = 0;
4433 
4434  int cnt = 0;
4435 
4436  const bool bOutputStatus = isatty(fileno(stderr));
4437 
4438  do {
4439  __FC_DECL__(dnaupd)(&IDO, &BMAT[0], &N, &WHICH[0], &NEV,
4440  &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
4441  &WORKD[0], &WORKL[0], &LWORKL, &INFO);
4442 
4443 #if 0
4444  std::cout << "cnt=" << cnt << ": IDO=" << IDO << ", INFO=" << INFO << std::endl;
4445 #endif
4446 
4447  // compute Y = OP*X
4448  ASSERT(IPNTR[0] - 1 <= 2*N);
4449  ASSERT(IPNTR[1] - 1 <= 2*N);
4450  MyVectorHandler X(N, &WORKD[IPNTR[0] - 1]);
4451  MyVectorHandler Y(N, &WORKD[IPNTR[1] - 1]);
4452 
4453  /*
4454  * NOTE: we are solving the problem
4455 
4456  MatB * X * Lambda = MatA * X
4457 
4458  * and we want to focus on Ritz parameters Lambda
4459  * as close as possible to (1., 0.), which maps
4460  * to (0., 0.) in continuous time.
4461  *
4462  * We are casting the problem in the form
4463 
4464  X * Alpha = A * X
4465 
4466  * by putting the problem in canonical form
4467 
4468  X * Lambda = MatB \ MatA * X
4469 
4470  * and then subtracting a shift Sigma = (1., 0) after :
4471 
4472  X * Lambda - X = MatB \ MatA * X - X
4473 
4474  X * (Lambda - 1.) = (MatB \ MatA - I) * X
4475 
4476  * so
4477 
4478  Alpha = Lambda - 1.
4479 
4480  A = MatB \ MatA - I
4481 
4482  * and the sequence of operations for Y = A * X is
4483 
4484  X' = MatA * X
4485  X'' = MatB \ X'
4486  Y = X'' - X
4487 
4488  * the eigenvalues need to be modified by adding 1.
4489  */
4490 
4491 
4492  MatA.MatVecMul(*sm.pResHdl(), X);
4493  sm.Solve();
4494  *sm.pSolHdl() -= X;
4495 
4496  Y = *sm.pSolHdl();
4497 
4498  static const int CNT = 100;
4499  cnt++;
4500  if (bOutputStatus && !(cnt % CNT)) {
4501  if (bNewLine && silent_err) {
4502  silent_cerr(std::endl);
4503  bNewLine = false;
4504  }
4505  silent_cerr("\r" "cnt=" << cnt);
4506  }
4507 
4509  if (bNewLine && silent_err) {
4510  silent_cerr(std::endl);
4511  bNewLine = false;
4512  }
4513  silent_cerr((cnt >= CNT ? "\n" : "")
4514  << "ARPACK: interrupted" << std::endl);
4515  return;
4516  }
4517  } while (IDO == 1 || IDO == -1);
4518 
4519  if (bOutputStatus) {
4520  if (bNewLine && silent_err) {
4521  silent_cerr(std::endl);
4522  bNewLine = false;
4523  }
4524  silent_cerr("\r" "cnt=" << cnt << std::endl);
4525  }
4526 
4527  // NOTE: improve diagnostics
4528  if (INFO < 0) {
4529  if (bNewLine && silent_err) {
4530  silent_cerr(std::endl);
4531  bNewLine = false;
4532  }
4533  silent_cerr("ARPACK error after " << cnt << " iterations; "
4534  "IDO=" << IDO << ", INFO=" << INFO << std::endl);
4535  return;
4536  }
4537 
4538  switch (INFO) {
4539  case 0:
4540  break;
4541 
4542  case 1:
4543  if (bNewLine && silent_err) {
4544  silent_cerr(std::endl);
4545  bNewLine = false;
4546  }
4547  silent_cerr("INFO=1: Maximum number of iterations taken. "
4548  "All possible eigenvalues of OP have been found. IPARAM(5) "
4549  "returns the number of wanted converged Ritz values "
4550  "(currently = " << IPARAM[4] << "; requested NEV = " << NEV << ")."
4551  << std::endl);
4552  break;
4553 
4554  case 2:
4555  if (bNewLine && silent_err) {
4556  silent_cerr(std::endl);
4557  bNewLine = false;
4558  }
4559  silent_cerr("INFO=2: No longer an informational error. Deprecated starting "
4560  "with release 2 of ARPACK."
4561  << std::endl);
4562  break;
4563 
4564  case 3:
4565  if (bNewLine && silent_err) {
4566  silent_cerr(std::endl);
4567  bNewLine = false;
4568  }
4569  silent_cerr("INFO=3: No shifts could be applied during a cycle of the "
4570  "implicitly restarted Arnoldi iteration. One possibility "
4571  "is to increase the size of NCV (currently = " << NCV << ") "
4572  "relative to NEV (currently = " << NEV << "). "
4573  "See remark 4 in dnaupd(3)."
4574  << std::endl);
4575  break;
4576 
4577  default:
4578  if (bNewLine && silent_err) {
4579  silent_cerr(std::endl);
4580  bNewLine = false;
4581  }
4582  silent_cerr("INFO=" << INFO << ": undocumented value." << std::endl);
4583  break;
4584  }
4585 
4586  std::ostream& Out = pDM->GetOutFile();
4587  Out << "INFO: " << INFO << std::endl;
4588 
4589 #if 0
4590  for (int ii = 0; ii < 11; ii++) {
4591  std::cerr << "IPARAM(" << ii + 1 << ")=" << IPARAM[ii] << std::endl;
4592  }
4593 #endif
4594 
4595  logical RVEC = true;
4596  const char *HOWMNY = "A";
4597  std::vector<logical> SELECT(NCV);
4598  std::vector<doublereal> D(2*NCV);
4599  doublereal *DR = &D[0], *DI = &D[NCV];
4600  std::vector<doublereal> Z(N*(NCV + 1));
4601  integer LDZ = N;
4602  std::vector<doublereal> WORKEV(3*NCV);
4603 
4604 #if 0
4605  std::cerr << "dneupd:" << std::endl
4606  << "RVEC = " << RVEC << std::endl
4607  << "HOWMNY = " << HOWMNY << std::endl
4608  << "SELECT(" << NCV << ")" << std::endl
4609  << "DR(" << NEV + 1 << ")" << std::endl
4610  << "DI(" << NEV + 1 << ")" << std::endl
4611  << "Z(" << N << ", " << NEV + 1 << ")" << std::endl
4612  << "LDZ = " << LDZ << std::endl
4613  << "SIGMAR = " << SIGMAR << std::endl
4614  << "SIGMAI = " << SIGMAI << std::endl
4615  << "WORKEV(" << 3*NCV << ")" << std::endl
4616  << "BMAT = " << BMAT << std::endl
4617  << "N = " << N << std::endl
4618  << "WHICH = " << WHICH << std::endl
4619  << "NEV = " << NEV << std::endl
4620  << "TOL = " << TOL << std::endl
4621  << "RESID(" << N << ")" << std::endl
4622  << "NCV = " << NCV << std::endl
4623  << "V(" << N << ", " << NCV << ")" << std::endl
4624  << "LDV = " << LDV << std::endl
4625  << "IPARAM(" << sizeof(IPARAM)/sizeof(IPARAM[0]) << ")" << std::endl
4626  << "IPNTR(" << sizeof(IPNTR)/sizeof(IPNTR[0]) << ")" << std::endl
4627  << "WORKD(" << 3*N << ")" << std::endl
4628  << "WORKL(" << LWORKL << ")" << std::endl
4629  << "LWORKL = " << LWORKL << std::endl
4630  << "INFO = " << INFO << std::endl;
4631 #endif
4632 
4633  __FC_DECL__(dneupd)(&RVEC, &HOWMNY[0], &SELECT[0], &DR[0], &DI[0],
4634  &Z[0], &LDZ, &SIGMAR, &SIGMAI, &WORKEV[0],
4635  &BMAT[0], &N, &WHICH[0], &NEV,
4636  &TOL, &RESID[0], &NCV, &V[0], &LDV, &IPARAM[0], &IPNTR[0],
4637  &WORKD[0], &WORKL[0], &LWORKL, &INFO);
4638 
4639  int nconv = IPARAM[4];
4640  if (nconv > 0) {
4641  ASSERT(nconv <= NCV);
4642  MyVectorHandler AlphaR(nconv, DR);
4643  MyVectorHandler AlphaI(nconv, DI);
4644  std::vector<bool> vOut(nconv);
4645  output_eigenvalues(0, AlphaR, AlphaI, 1., pDM, pEA, 1, nconv, vOut);
4646 
4648  pDM->OutputEigGeometry(uCurr, pEA->iResultsPrecision);
4649  }
4650 
4652  std::vector<doublereal *> ZC(nconv);
4653  FullMatrixHandler VR(&Z[0], &ZC[0], N*nconv, N, nconv);
4654  pDM->OutputEigenvectors(0, AlphaR, AlphaI, 1.,
4655  0, VR, vOut, uCurr, pEA->iResultsPrecision);
4656  }
4657 
4658  } else {
4659  if (bNewLine && silent_err) {
4660  silent_cerr(std::endl);
4661  bNewLine = false;
4662  }
4663  silent_cerr("no converged Ritz coefficients" << std::endl);
4664  }
4665 }
4666 #endif // USE_ARPACK
4667 
4668 #ifdef USE_JDQZ
4669 // Computes eigenvalues and eigenvectors using ARPACK's
4670 // canonical non-symmetric eigenanalysis
4671 static void
4672 eig_jdqz(const MatrixHandler *pMatA, const MatrixHandler *pMatB,
4674  bool bNewLine, const unsigned uCurr)
4675 {
4676  const NaiveMatrixHandler& MatA = dynamic_cast<const NaiveMatrixHandler &>(*pMatA);
4677  const NaiveMatrixHandler& MatB = dynamic_cast<const NaiveMatrixHandler &>(*pMatB);
4678 
4679  MBJDQZ mbjdqz(MatA, MatB);
4680  mbjdqzp = &mbjdqz;
4681 
4682  /*
4683 
4684 alpha, beta Obvious from equation (1)
4685 wanted Compute the converged eigenvectors (if wanted =
4686  .true.)
4687 eivec Converged eigenvectors if wanted = .true., else con-
4688  verged Schur vectors
4689 n The size of the problem
4690 target The value near which the eigenvalues are sought
4691 eps Tolerance of the eigensolutions, Ax-Bx /|/| <
4692 kmax Number of wanted eigensolutions, on output: number of
4693  converged eigenpairs
4694 jmax Maximum size of the search space
4695 jmin Minimum size of the search space
4696 method Linear equation solver:
4697  1: GMRESm , [2]
4698  2: BiCGstab(), [3]
4699 m Maximum dimension of searchspace of GMRESm
4700 l Degree of GMRES-polynomial in Bi-CGstab()
4701 mxmv Maximum number of matrix-vector multiplications in
4702  GMRESm or BiCGstab()
4703 maxstep Maximum number of Jacobi-Davidson iterations
4704 lock Tracking parameter (section 2.5.1)
4705 order Selection criterion for Ritz values:
4706  0: nearest to target
4707  -1: smallest real part
4708  1: largest real part
4709  -2: smallest imaginary part
4710  2: largest imaginary part
4711 testspace Determines how to expand the testspace W
4712  1: w = "Standard Petrov" ×v (Section 3.1.1)
4713  2: w = "Standard 'variable' Petrov" ×v (Section 3.1.2)
4714  3: w = "Harmonic Petrov" ×v (Section 3.5.1)
4715 zwork Workspace
4716 lwork Size of workspace, >= 4+m+5jmax+3kmax if GMRESm
4717  is used, >= 10 + 6 + 5jmax + 3kmax if Bi-CGstab() is
4718  used.
4719 
4720 */
4721 
4722  std::vector<doublecomplex> alpha;
4723  std::vector<doublecomplex> beta;
4724  std::vector<doublecomplex> eivec;
4725  logical wanted = 1;
4726  integer n = pMatA->iGetNumRows();
4727  doublecomplex target = { 1., 0. };
4728  doublereal eps = pEA->jdqz.eps;
4729  integer kmax = pEA->jdqz.kmax;
4730  integer jmax = pEA->jdqz.jmax;
4731  integer jmin = pEA->jdqz.jmin;
4732  integer method = pEA->jdqz.method;
4733  integer m = pEA->jdqz.m;
4734  integer l = pEA->jdqz.l;
4735  integer mxmv = pEA->jdqz.mxmv;
4736  integer maxstep = pEA->jdqz.maxstep;
4737  doublereal lock = pEA->jdqz.lock;
4738  integer order = pEA->jdqz.order;
4739  integer testspace = pEA->jdqz.testspace;
4740  std::vector<doublecomplex> zwork;
4741  integer lwork;
4742 
4743  switch (method) {
4745  lwork = 4 + m + 5*jmax + 3*kmax;
4746  break;
4747 
4749  lwork = 10 + 6*l + 5*jmax + 3*kmax;
4750  break;
4751 
4752  default:
4754  }
4755 
4756  alpha.resize(jmax);
4757  beta.resize(jmax);
4758  eivec.resize(n*kmax);
4759  zwork.resize(n*lwork);
4760 
4761  __FC_DECL__(jdqz)(
4762  &alpha[0],
4763  &beta[0],
4764  &eivec[0],
4765  &wanted,
4766  &n,
4767  &target,
4768  &eps,
4769  &kmax,
4770  &jmax,
4771  &jmin,
4772  &method,
4773  &m,
4774  &l,
4775  &mxmv,
4776  &maxstep,
4777  &lock,
4778  &order,
4779  &testspace,
4780  &zwork[0],
4781  &lwork);
4782 
4783  if (bNewLine && silent_err) {
4784  silent_cerr(std::endl);
4785  bNewLine = false;
4786  }
4787  silent_cerr("\r" "cnt=" << mbjdqz.Cnt() << std::endl);
4788 
4789  if (kmax > 0) {
4790  int nconv = kmax;
4791  MyVectorHandler AlphaR(nconv);
4792  MyVectorHandler AlphaI(nconv);
4793  MyVectorHandler Beta(nconv);
4794  for (integer c = 0; c < nconv; c++) {
4795  if (beta[c].i != 0.) {
4796  doublereal d = std::sqrt(beta[c].r*beta[c].r + beta[c].i*beta[c].i);
4797  Beta(c + 1) = d;
4798  AlphaR(c + 1) = (alpha[c].r*beta[c].r + alpha[c].i*beta[c].i)/d;
4799  AlphaI(c + 1) = (alpha[c].i*beta[c].r - alpha[c].r*beta[c].i)/d;
4800 
4801  } else {
4802  Beta(c + 1) = beta[c].r;
4803  AlphaR(c + 1) = alpha[c].r;
4804  AlphaI(c + 1) = alpha[c].i;
4805  }
4806  }
4807  std::vector<bool> vOut(nconv);
4808  output_eigenvalues(0, AlphaR, AlphaI, 1., pDM, pEA, 1, nconv, vOut);
4809 
4811  pDM->OutputGeometry(pEA->iResultsPrecision);
4812  }
4813 
4815  FullMatrixHandler VR(n, nconv);
4816  doublecomplex *p = &eivec[0] - 1;
4817  for (integer c = 1; c <= nconv; c++) {
4818 
4819  if (AlphaI(c) == 0.) {
4820  for (integer r = 1; r <= n; r++) {
4821  VR(r, c) = p[r].r;
4822  }
4823 
4824  } else {
4825  for (integer r = 1; r <= n; r++) {
4826  VR(r, c) = p[r].r;
4827  VR(r, c + 1) = p[n + r].i;
4828  }
4829 
4830  p += n;
4831  c++;
4832  }
4833 
4834  p += n;
4835  }
4836 
4837  pDM->OutputEigenvectors(&Beta, AlphaR, AlphaI, 1.,
4838  0, VR, vOut, uCurr, pEA->iResultsPrecision);
4839  }
4840 
4841  } else {
4842  if (bNewLine && silent_err) {
4843  silent_cerr(std::endl);
4844  bNewLine = false;
4845  }
4846  silent_cerr("no converged eigenpairs" << std::endl);
4847  }
4848 }
4849 #endif // USE_JDQZ
4850 
4851 // Driver for eigenanalysis
4852 void
4853 Solver::Eig(bool bNewLine)
4854 {
4855  DEBUGCOUTFNAME("Solver::Eig");
4856 
4857  /*
4858  * MatA, MatB: MatrixHandlers to eigenanalysis matrices
4859  * MatVL, MatVR: MatrixHandlers to eigenvectors, if required
4860  * AlphaR, AlphaI Beta: eigenvalues
4861  * WorkVec: Workspace
4862  * iWorkSize: Size of the workspace
4863  */
4864 
4865  DEBUGCOUT("Solver::Eig(): performing eigenanalysis" << std::endl);
4866 
4867  integer iSize = iNumDofs;
4868 
4869  SolutionManager *pSM = 0;
4870  MatrixHandler *pMatA = 0;
4871  MatrixHandler *pMatB = 0;
4872 
4875  FullMatrixHandler(iSize));
4876 
4878  FullMatrixHandler(iSize));
4879 
4884  NaiveMatrixHandler(iSize));
4885  pMatB = pSM->pMatHdl();
4886 
4887  } else if (EigAn.uFlags & EigenAnalysis::EIG_USE_JDQZ) {
4889  NaiveMatrixHandler(iSize));
4891  NaiveMatrixHandler(iSize));
4892 
4895  SpMapMatrixHandler(iSize));
4897  SpMapMatrixHandler(iSize));
4898  }
4899 
4900  pMatA->Reset();
4901  pMatB->Reset();
4902 
4903  // Matrices assembly (see eig.ps)
4904  doublereal h = EigAn.dParam;
4905  pDM->AssJac(*pMatA, -h/2.);
4906  pDM->AssJac(*pMatB, h/2.);
4907 
4908 #ifdef DEBUG
4909  DEBUGCOUT(std::endl
4910  << "Matrix A:" << std::endl << *pMatA << std::endl
4911  << "Matrix B:" << std::endl << *pMatB << std::endl);
4912 #endif /* DEBUG */
4913 
4914  unsigned uCurr = EigAn.currAnalysis - EigAn.Analyses.begin();
4916  unsigned uSize = EigAn.Analyses.size();
4917  if (uSize >= 1) {
4918  std::stringstream postfix_ss;
4919 
4920  postfix_ss << '_';
4921 
4922  if (EigAn.iFNameWidth > 0) {
4923  postfix_ss << std::setw(EigAn.iFNameWidth) << std::setfill('0') << uCurr;
4924 
4925  } else if (!EigAn.iFNameFormat.empty()) {
4926  char buf[BUFSIZ];
4927 
4928  snprintf(buf, sizeof(buf), EigAn.iFNameFormat.c_str(), (int)uCurr);
4929  postfix_ss << buf;
4930 
4931  } else {
4932  postfix_ss << uCurr;
4933  }
4934 
4935  pDM->OutputEigOpen(postfix_ss.str());
4936  pDM->OutputEigParams(dTime, h/2., uCurr, EigAn.iResultsPrecision);
4937  }
4938  }
4940  pDM->OutputEigFullMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4941  }
4943  if (dynamic_cast<const NaiveMatrixHandler *>(pMatB)) {
4944  pDM->OutputEigNaiveMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4945  } else {
4946  pDM->OutputEigSparseMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
4947  }
4948  }
4949 
4951 #ifdef USE_LAPACK
4953  eig_lapack(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
4954  break;
4955 #endif // USE_LAPACK
4956 
4957 #ifdef USE_ARPACK
4959  eig_arpack(pMatA, pSM, pDM, &EigAn, bNewLine, uCurr);
4960  break;
4961 #endif // USE_ARPACK
4962 
4963 #ifdef USE_JDQZ
4965  eig_jdqz(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
4966  break;
4967 #endif // USE_JDQZ
4968 
4969  default:
4970  // only output matrices, use external eigenanalysis
4971  break;
4972  }
4973 
4974  pDM->OutputEigClose();
4975 
4976  if (pSM) {
4977  pMatB = 0;
4978  SAFEDELETE(pSM);
4979  }
4980 
4981  if (pMatA) {
4982  SAFEDELETE(pMatA);
4983  }
4984 
4985  if (pMatB) {
4986  SAFEDELETE(pMatB);
4987  }
4988 }
4989 
4991 {
4992  if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
4994  }
4995 
4996  // The same behavior like in previous releases
4997  return MaxTimeStep.dGet();
4998 }
4999 
5000 SolutionManager *const
5002 {
5003  SolutionManager *pCurrSM = CurrLinearSolver.GetSolutionManager(iNLD, iLWS);
5004 
5005  /* special extra parameters if required */
5006  switch (CurrLinearSolver.GetSolver()) {
5008 #ifdef HAVE_UMFPACK_TIC_DISABLE
5009  if (pRTSolver) {
5010  /* disable profiling, to avoid times() system call
5011  *
5012  * This function has been introduced in Umfpack 4.1
5013  * by our patch at
5014  *
5015  * http://mbdyn.aero.polimi.it/~masarati/Download/\
5016  * mbdyn/umfpack-4.1-nosyscalls.patch
5017  *
5018  * but since Umfpack 4.3 is no longer required,
5019  * provided the library is compiled with -DNO_TIMER
5020  * to disable run-time syscalls to timing routines.
5021  */
5022  umfpack_tic_disable();
5023  }
5024 #endif // HAVE_UMFPACK_TIC_DISABLE
5025  break;
5026 
5027  default:
5028  break;
5029  }
5030 
5031  return pCurrSM;
5032 };
5033 
5034 
5035 SolutionManager *const
5037 {
5038  SolutionManager *pSSM(0);
5039 
5040 #ifdef USE_MPI
5041  switch (CurrIntSolver.GetSolver()) {
5042  case LinSol::LAPACK_SOLVER:
5044  case LinSol::NAIVE_SOLVER:
5046  case LinSol::Y12_SOLVER:
5047  break;
5048 
5049  default:
5050  silent_cerr("apparently solver "
5052  << " is not allowed as interface solver "
5053  "for SchurSolutionManager" << std::endl);
5055  }
5056 
5060  iNumLocDofs,
5063 
5064 #else /* !USE_MPI */
5065  silent_cerr("Configure --with-mpi to enable Schur solver" << std::endl);
5067 #endif /* !USE_MPI */
5068 
5069  return pSSM;
5070 };
5071 
5072 NonlinearSolver *const
5074 {
5075  NonlinearSolver *pNLS = 0;
5076 
5077  switch (NonlinearSolverType) {
5079  switch (MFSolverType) {
5082  BiCGStab,
5083  BiCGStab(PcType,
5084  iPrecondSteps,
5085  dIterTol,
5089  *this));
5090  break;
5091 
5092  default:
5093  pedantic_cout("unknown matrix free solver type; "
5094  "using default" << std::endl);
5095  /* warning: should be unreachable */
5096 
5099  Gmres,
5100  Gmres(PcType,
5101  iPrecondSteps,
5102  dIterTol,
5106  *this));
5107  break;
5108  }
5109  break;
5110 
5111  default:
5112  pedantic_cout("unknown nonlinear solver type; using default"
5113  << std::endl);
5114 
5119  bKeepJac,
5121  *this));
5122  break;
5126  LineSearchSolver(pDM,
5127  *this,
5128  LineSearch));
5129  break;
5130  }
5131  return pNLS;
5132 }
5133 
5134 void
5135 Solver::SetupSolmans(integer iStates, bool bCanBeParallel)
5136 {
5137  DEBUGLCOUT(MYDEBUG_MEM, "creating SolutionManager\n\tsize = "
5138  << iNumDofs*iUnkStates <<
5139  "\n\tnumdofs = " << iNumDofs
5140  << "\n\tnumstates = " << iStates << std::endl);
5141 
5142  /* delete previous solmans */
5143  if (pSM != 0) {
5144  SAFEDELETE(pSM);
5145  pSM = 0;
5146  }
5147  if (pLocalSM != 0) {
5149  pLocalSM = 0;
5150  }
5151 
5152  integer iWorkSpaceSize = CurrLinearSolver.iGetWorkSpaceSize();
5153  integer iLWS = iWorkSpaceSize;
5154  integer iNLD = iNumDofs*iStates;
5155  if (bCanBeParallel && bParallel) {
5156  /* FIXME BEPPE! */
5157  iLWS = iWorkSpaceSize*iNumLocDofs/(iNumDofs*iNumDofs);
5158  /* FIXME: GIUSTO QUESTO? */
5159  iNLD = iNumLocDofs*iStates;
5160  }
5161 
5162  SolutionManager *pCurrSM = AllocateSolman(iNLD, iLWS);
5163 
5164  /*
5165  * This is the LOCAL solver if instantiating a parallel
5166  * integrator; otherwise it is the MAIN solver
5167  */
5168  if (bCanBeParallel && bParallel) {
5169  pLocalSM = pCurrSM;
5170 
5171  /* Crea il solutore di Schur globale */
5172  pSM = AllocateSchurSolman(iStates);
5173 
5174  } else {
5175  pSM = pCurrSM;
5176  }
5177  /*
5178  * FIXME: at present there MUST be a pSM
5179  * (even for matrix-free nonlinear solvers)
5180  */
5181  if (pSM == 0) {
5182  silent_cerr("No linear solver defined" << std::endl);
5184  }
5185 }
5186 
5187 clock_t
5189 {
5190  return pDM->GetCPUTime();
5191 }
5192 
5193 void
5194 Solver::PrintResidual(const VectorHandler& Res, integer iIterCnt) const
5195 {
5196  pDM->PrintResidual(Res, iIterCnt);
5197 }
5198 
5199 void
5200 Solver::PrintSolution(const VectorHandler& Sol, integer iIterCnt) const
5201 {
5202  pDM->PrintSolution(Sol, iIterCnt);
5203 }
5204 
5206 {
5207  if (pDerivativeSteps) {
5208  // Time step cannot be reduced
5209  return;
5210  }
5211 
5212  if (dErr > dMaxResidual) {
5213  if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
5214  if (outputIters()) {
5215 #ifdef USE_MPI
5216  if (!bParallel || MBDynComm.Get_rank() == 0)
5217 #endif /* USE_MPI */
5218  {
5219  silent_cerr("warning: current residual = " << dErr
5220  << " > maximum residual = " << dMaxResidual
5221  << std::endl);
5222  }
5223  }
5224 
5226  }
5227  }
5228 
5229  if (dErrDiff > dMaxResidualDiff) {
5230  if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
5231  if (outputIters()) {
5232 #ifdef USE_MPI
5233  if (!bParallel || MBDynComm.Get_rank() == 0)
5234 #endif /* USE_MPI */
5235  {
5236  silent_cerr("warning: current residual for differential equations = " << dErrDiff
5237  << " > maximum residual = " << dMaxResidualDiff
5238  << std::endl);
5239  }
5240  }
5241 
5243  }
5244  }
5245 
5246  switch (eTimeStepLimit) {
5247  case TS_SOFT_LIMIT:
5248  break;
5249 
5250  case TS_HARD_LIMIT: {
5251  const doublereal dMaxTS = MaxTimeStep.dGet();
5252 
5253  if (dCurrTimeStep > dMaxTS && dCurrTimeStep > dMinTimeStep) {
5254  if (outputIters()) {
5255 #ifdef USE_MPI
5256  if (!bParallel || MBDynComm.Get_rank() == 0)
5257 #endif /* USE_MPI */
5258  {
5259  silent_cerr("warning: current time step = "
5260  << dCurrTimeStep
5261  << " > hard limit of the maximum time step = "
5262  << dMaxTS << std::endl);
5263  }
5264  }
5265