#include <solver.h>

Inheritance diagram for Solver:

Collaboration diagram for Solver:

Classes
struct	EigenAnalysis

class	EndOfSimulation

class	ErrGeneric

class	ErrMaxIterations

class	SimulationDiverged

Public Member Functions
	Solver (MBDynParser &HP, const std::string &sInputFileName, const std::string &sOutputFileName, unsigned int nThreads, bool bParallel=false)

virtual	~Solver (void)

virtual bool	Prepare (void)

virtual bool	Start (void)

virtual bool	Advance (void)

void	Run (void)

std::ostream &	Restart (std::ostream &out, DataManager::eRestart type) const

virtual DataManager *	pGetDataManager (void) const

virtual SolutionManager *	pGetSolutionManager (void) const

virtual const LinSol &	GetLinearSolver (void) const

virtual NonlinearSolver *	pGetNonlinearSolver (void) const

virtual StepIntegrator *	pGetStepIntegrator (void) const

virtual doublereal	dGetFinalTime (void) const

virtual doublereal	dGetInitialTimeStep (void) const

virtual clock_t	GetCPUTime (void) const

virtual void	PrintResidual (const VectorHandler &Res, integer iIterCnt) const

virtual void	PrintSolution (const VectorHandler &Sol, integer iIterCnt) const

virtual void	CheckTimeStepLimit (doublereal dErr, doublereal dErrDiff) const throw (NonlinearSolver::TimeStepLimitExceeded, NonlinearSolver::MaxResidualExceeded)

Public Member Functions inherited from SolverDiagnostics
	SolverDiagnostics (unsigned OF=OUTPUT_DEFAULT, DriveCaller *pOM=0)

virtual	~SolverDiagnostics (void)

void	SetNoOutput (void)

void	SetOutputMeter (DriveCaller *pOM)

void	SetOutputDriveHandler (const DriveHandler *pDH)

void	SetOutputFlags (unsigned OF)

void	AddOutputFlags (unsigned OF)

void	DelOutputFlags (unsigned OF)

MatrixHandler::Norm_t	GetCondMatNorm (void) const

bool	outputMeter (void) const

bool	outputIters (void) const

bool	outputRes (void) const

bool	outputSol (void) const

bool	outputJac (void) const

bool	outputStep (void) const

bool	outputBailout (void) const

bool	outputCounter (void) const

bool	outputMatrixConditionNumber (void) const

bool	outputSolverConditionNumber (void) const

bool	outputSolverConditionStat (void) const

bool	outputCPUTime (void) const

bool	outputMsg (void) const

Protected Types
enum	TimeStepFlags { TS_SOFT_LIMIT = 0, TS_HARD_LIMIT = 1 }

enum	AbortAfter { AFTER_UNKNOWN, AFTER_INPUT, AFTER_ASSEMBLY, AFTER_DERIVATIVES, AFTER_DUMMY_STEPS }

enum	StepIntegratorType { INT_CRANKNICOLSON, INT_MODCRANKNICOLSON, INT_MS2, INT_HOPE, INT_THIRDORDER, INT_IMPLICITEULER, INT_UNKNOWN }

enum	{ SOLVER_STATUS_UNINITIALIZED, SOLVER_STATUS_PREPARED, SOLVER_STATUS_STARTED }

Protected Types inherited from SolverDiagnostics
enum	{ OUTPUT_NONE = 0x0000, OUTPUT_ITERS = 0x0001, OUTPUT_RES = 0x0002, OUTPUT_SOL = 0x0004, OUTPUT_JAC = 0x0008, OUTPUT_BAILOUT = 0x0010, OUTPUT_MSG = 0x0020, OUTPUT_COUNTER = 0x0040, OUTPUT_MAT_COND_NUM_1 = 0x0080, OUTPUT_MAT_COND_NUM_INF = 0x0100, OUTPUT_SOLVER_COND_NUM = 0x0200, OUTPUT_SOLVER_COND_STAT = 0x400, OUTPUT_CPU_TIME = 0x800, OUTPUT_MAT_COND_NUM = OUTPUT_MAT_COND_NUM_1 \| OUTPUT_MAT_COND_NUM_INF, OUTPUT_DEFAULT = OUTPUT_MSG, OUTPUT_STEP = (OUTPUT_ITERS \| OUTPUT_RES \| OUTPUT_SOL \| OUTPUT_JAC \| OUTPUT_MAT_COND_NUM \| OUTPUT_SOLVER_COND_NUM), OUTPUT_MASK = 0x07FF }

Protected Types inherited from NonlinearSolverOptions
enum	ScaleFlags { SCALE_ALGEBRAIC_EQUATIONS_NO = 0, SCALE_ALGEBRAIC_EQUATIONS_YES = 1 }

Protected Member Functions
void	Eig (bool bNewLine=false)

void	Flip (void)

void	ReadData (MBDynParser &HP)

SolutionManager *const	AllocateSolman (integer iNLD, integer iLWS=0)

SolutionManager *const	AllocateSchurSolman (integer iStates)

NonlinearSolver *const	AllocateNonlinearSolver ()

void	SetupSolmans (integer iStates, bool bCanBeParallel=false)

doublereal	dGetInitialMaxTimeStep () const

Protected Member Functions inherited from NonlinearSolverOptions
	NonlinearSolverOptions (bool bHonorJacRequest=false, enum ScaleFlags eScaleFlags=SCALE_ALGEBRAIC_EQUATIONS_NO, doublereal dScaleAlgebraic=1.)

Protected Attributes
TimeStepControl *	pTSC

doublereal	dCurrTimeStep

integer	iStIter

doublereal	dTime

DriveOwner	MaxTimeStep

doublereal	dMinTimeStep

StepIntegrator::StepChange	CurrStep

std::string	sInputFileName

std::string	sOutputFileName

MBDynParser &	HP

integer	iMaxIterations

struct EigenAnalysis	EigAn

RTSolverBase *	pRTSolver

integer	iNumPreviousVectors

integer	iUnkStates

doublereal *	pdWorkSpace

std::deque< MyVectorHandler * >	qX

std::deque< MyVectorHandler * >	qXPrime

MyVectorHandler *	pX

MyVectorHandler *	pXPrime

doublereal	dInitialTime

doublereal	dFinalTime

doublereal	dRefTimeStep

doublereal	dInitialTimeStep

doublereal	dMaxResidual

doublereal	dMaxResidualDiff

enum Solver::TimeStepFlags	eTimeStepLimit

integer	iDummyStepsNumber

doublereal	dDummyStepsRatio

AbortAfter	eAbortAfter

StepIntegratorType	RegularType

StepIntegratorType	DummyType

StepIntegrator *	pDerivativeSteps

StepIntegrator *	pFirstDummyStep

StepIntegrator *	pDummySteps

StepIntegrator *	pFirstRegularStep

StepIntegrator *	pRegularSteps

StepIntegrator *	pCurrStepIntegrator

DriveCaller *	pRhoRegular

DriveCaller *	pRhoAlgebraicRegular

DriveCaller *	pRhoDummy

DriveCaller *	pRhoAlgebraicDummy

doublereal	dDerivativesCoef

LinSol	CurrLinearSolver

NonlinearSolverTest::Type	ResTest

NonlinearSolverTest::Type	SolTest

bool	bScale

MyVectorHandler	Scale

bool	bTrueNewtonRaphson

bool	bKeepJac

integer	iIterationsBeforeAssembly

NonlinearSolver::Type	NonlinearSolverType

MatrixFreeSolver::SolverType	MFSolverType

doublereal	dIterTol

Preconditioner::PrecondType	PcType

integer	iPrecondSteps

integer	iIterativeMaxSteps

doublereal	dIterertiveEtaMax

doublereal	dIterertiveTau

struct LineSearchParameters	LineSearch

bool	bParallel

SchurDataManager *	pSDM

LinSol	CurrIntSolver

integer	iNumLocDofs

integer	iNumIntDofs

integer *	pLocDofs

integer *	pIntDofs

Dof *	pDofs

SolutionManager *	pLocalSM

DataManager *	pDM

integer	iNumDofs

SolutionManager *	pSM

NonlinearSolver *	pNLS

enum Solver:: { ... }	eStatus

bool	bOutputCounter

std::string	outputCounterPrefix

std::string	outputCounterPostfix

integer	iTotIter

doublereal	dTotErr

doublereal	dTest

doublereal	dSolTest

bool	bSolConv

bool	bOut

long	lStep

Protected Attributes inherited from SolverDiagnostics
unsigned	OutputFlags

DriveCaller *	pOutputMeter

Protected Attributes inherited from NonlinearSolverOptions
bool	bHonorJacRequest

enum NonlinearSolverOptions::ScaleFlags	eScaleFlags

doublereal	dScaleAlgebraic

Detailed Description

Definition at line 78 of file solver.h.

Member Enumeration Documentation

anonymous enum

protected

Enumerator
SOLVER_STATUS_UNINITIALIZED
SOLVER_STATUS_PREPARED
SOLVER_STATUS_STARTED

Definition at line 355 of file solver.h.

              {
                 SOLVER_STATUS_UNINITIALIZED,
                 SOLVER_STATUS_PREPARED,
                 SOLVER_STATUS_STARTED
         } eStatus;

enum Solver::AbortAfter

protected

Enumerator
AFTER_UNKNOWN
AFTER_INPUT
AFTER_ASSEMBLY
AFTER_DERIVATIVES
AFTER_DUMMY_STEPS

Definition at line 257 of file solver.h.

                         {
                 AFTER_UNKNOWN,
                 AFTER_INPUT,
                 AFTER_ASSEMBLY,
                 AFTER_DERIVATIVES,
                 AFTER_DUMMY_STEPS
         };

enum Solver::StepIntegratorType

protected

Enumerator
INT_CRANKNICOLSON
INT_MODCRANKNICOLSON
INT_MS2
INT_HOPE
INT_THIRDORDER
INT_IMPLICITEULER
INT_UNKNOWN

Definition at line 267 of file solver.h.

                                 {
                         INT_CRANKNICOLSON,
                         INT_MODCRANKNICOLSON,
                         INT_MS2,
                         INT_HOPE,
                         INT_THIRDORDER,
                         INT_IMPLICITEULER,
                         INT_UNKNOWN
         };

enum Solver::TimeStepFlags

protected

Enumerator
TS_SOFT_LIMIT
TS_HARD_LIMIT

Definition at line 247 of file solver.h.

                            {
                 TS_SOFT_LIMIT = 0,
                 TS_HARD_LIMIT = 1
         } eTimeStepLimit;

Constructor & Destructor Documentation

Solver::Solver	(	MBDynParser &	HP,
		const std::string &	sInputFileName,
		const std::string &	sOutputFileName,
		unsigned int	nThreads,
		bool	bParallel = `false`
	)

Definition at line 274 of file solver.cc.

References ASSERT, DEBUGCOUTFNAME, and InitTimeStepData().

 :
 #ifdef USE_MULTITHREAD
 nThreads(nThreads),
 #endif /* USE_MULTITHREAD */
 pTSC(0),
 dCurrTimeStep(0.),
 iStIter(0),
 dTime(0.),
 MaxTimeStep(new ConstDriveCaller(::dDefaultMaxTimeStep)),
 dMinTimeStep(::dDefaultMinTimeStep),
 CurrStep(StepIntegrator::NEWSTEP),
 sInputFileName(sInFName),
 sOutputFileName(sOutFName),
 HP(HPar),
 iMaxIterations(::iDefaultMaxIterations),
 EigAn(),
 pRTSolver(0),
 iNumPreviousVectors(2),
 iUnkStates(1),
 pdWorkSpace(0),
 qX(),
 qXPrime(),
 pX(0),
 pXPrime(0),
 dInitialTime(0.),
 dFinalTime(0.),
 dRefTimeStep(0.),
 dInitialTimeStep(1.),
 dMaxResidual(std::numeric_limits<doublereal>::max()),
 dMaxResidualDiff(std::numeric_limits<doublereal>::max()),
 eTimeStepLimit(TS_SOFT_LIMIT),
 iDummyStepsNumber(::iDefaultDummyStepsNumber),
 dDummyStepsRatio(::dDefaultDummyStepsRatio),
 eAbortAfter(AFTER_UNKNOWN),
 RegularType(INT_UNKNOWN),
 DummyType(INT_UNKNOWN),
 pDerivativeSteps(0),
 pFirstDummyStep(0),
 pDummySteps(0),
 pFirstRegularStep(0),
 pRegularSteps(0),
 pCurrStepIntegrator(0),
 pRhoRegular(0),
 pRhoAlgebraicRegular(0),
 pRhoDummy(0),
 pRhoAlgebraicDummy(0),
 dDerivativesCoef(::dDefaultDerivativesCoefficient),
 CurrLinearSolver(),
 ResTest(NonlinearSolverTest::NORM),
 SolTest(NonlinearSolverTest::NONE),
 bScale(false),
 bTrueNewtonRaphson(true),
 bKeepJac(false),
 iIterationsBeforeAssembly(0),
 NonlinearSolverType(NonlinearSolver::UNKNOWN),
 /* for matrix-free solvers */
 MFSolverType(MatrixFreeSolver::UNKNOWN),
 dIterTol(::dDefaultTol),
 PcType(Preconditioner::FULLJACOBIANMATRIX),
 iPrecondSteps(::iDefaultPreconditionerSteps),
 iIterativeMaxSteps(::iDefaultPreconditionerSteps),
 dIterertiveEtaMax(defaultIterativeEtaMax),
 dIterertiveTau(defaultIterativeTau),
 /* end of matrix-free solvers */
 /* for line search solver */
 LineSearch(),
 /* end of line search solver */
 /* for parallel solvers */
 bParallel(bPar),
 pSDM(0),
 iNumLocDofs(0),
 iNumIntDofs(0),
 pLocDofs(0),
 pIntDofs(0),
 pDofs(0),
 pLocalSM(0),
 /* end of parallel solvers */
 pDM(0),
 iNumDofs(0),
 pSM(0),
 pNLS(0),
 eStatus(SOLVER_STATUS_UNINITIALIZED),
 bOutputCounter(false),
 outputCounterPrefix(),
 outputCounterPostfix(),
 iTotIter(0),
 dTotErr(0.),
 dTest(std::numeric_limits<double>::max()),
 dSolTest(std::numeric_limits<double>::max()),
 bSolConv(false),
 bOut(false),
 lStep(0)
 {
         DEBUGCOUTFNAME("Solver::Solver");
         ::InitTimeStepData();
         ASSERT(!sInFName.empty());
 }

Here is the call graph for this function:

Solver::~Solver ( void )

virtual

Definition at line 1647 of file solver.cc.

References DEBUGCOUTFNAME, DestroyTimeStepData(), iNumPreviousVectors, pDerivativeSteps, pDM, pDummySteps, pdWorkSpace, pFirstDummyStep, pFirstRegularStep, pNLS, pRegularSteps, pRTSolver, pSM, pTSC, pX, pXPrime, qX, qXPrime, SAFEDELETE, and SAFEDELETEARR.

 {
         DEBUGCOUTFNAME("Solver::~Solver");
 
         if (!qX.empty()) {
                 for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
                         if (qX[ivec] != 0) {
                                 SAFEDELETE(qX[ivec]);
                                 SAFEDELETE(qXPrime[ivec]);
                         }
                 }
         }
 
         if (pX) {
                 SAFEDELETE(pX);
         }
 
         if (pXPrime) {
                 SAFEDELETE(pXPrime);
         }
 
         if (pdWorkSpace) {
                 SAFEDELETEARR(pdWorkSpace);
         }
 
         if (pDM) {
                 SAFEDELETE(pDM);
         }
 
         if (pRTSolver) {
                 SAFEDELETE(pRTSolver);
         }
 
         if (pDerivativeSteps) {
                 SAFEDELETE(pDerivativeSteps);
         }
 
         if (pFirstDummyStep) {
                 SAFEDELETE(pFirstDummyStep);
         }
 
         if (pDummySteps) {
                 SAFEDELETE(pDummySteps);
         }
 
         if (pFirstRegularStep) {
                 SAFEDELETE(pFirstRegularStep);
         }
 
         if (pRegularSteps) {
                 SAFEDELETE(pRegularSteps);
         }
 
         if (pSM) {
                 SAFEDELETE(pSM);
         }
 
         if (pNLS) {
                 SAFEDELETE(pNLS);
         }
 
         if (pTSC) {
                 delete pTSC;
         }
 
         DestroyTimeStepData();
 }

Here is the call graph for this function:

Member Function Documentation

bool Solver::Advance ( void )

virtual

Reimplemented in InverseSolver.

Definition at line 1376 of file solver.cc.

References StepIntegrator::Advance(), Solver::EigenAnalysis::Analyses, Solver::EigenAnalysis::bAnalysis, DataManager::BeforePredict(), bOut, bOutputCounter, bSolConv, Solver::EigenAnalysis::currAnalysis, CurrStep, dCurrTimeStep, DEBUGCOUT, DEBUGCOUTFNAME, dFinalTime, TimeStepControl::dGetNewStepTime(), dMinTimeStep, dRefTimeStep, dSolTest, dTest, dTime, dTotErr, Eig(), EigAn, DataManager::EndOfSimulation(), eStatus, Flip(), DataManager::GetDofDescription(), StepIntegrator::GetIntegratorMaxIters(), DataManager::GetOutFile(), LinearSolver::ErrFactor::iCol, RTSolverBase::IsStopCommanded(), iStIter, iTotIter, RTSolverBase::Log(), lStep, MBDYN_EXCEPT_ARGS, mbdyn_stop_at_end_of_time_step(), StepIntegrator::NEWSTEP, DataManager::Output(), SolverDiagnostics::outputCounter(), outputCounterPostfix, outputCounterPrefix, SolverDiagnostics::outputMsg(), SolverDiagnostics::outputStep(), pDM, pNLS, pRegularSteps, pRTSolver, pTSC, pX, pXPrime, qX, qXPrime, StepIntegrator::REPEATSTEP, DataManager::SetTime(), SOLVER_STATUS_STARTED, sOutputFileName, RTSolverBase::StopCommanded(), NonlinearSolver::TotalAssembledJacobian(), RTSolverBase::Wait(), and MBDynErrBase::what().

Referenced by mb_sol_advance(), and Run().

 {
         DEBUGCOUTFNAME("Solver::Advance");
 
         // consistency check
         if (eStatus != SOLVER_STATUS_STARTED) {
                 silent_cerr("Started() must be called first" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         CurrStep = StepIntegrator::NEWSTEP;
 
         if (pDM->EndOfSimulation() || dTime >= dFinalTime) {
                 if (pRTSolver) {
                         pRTSolver->StopCommanded();
                 }
                 silent_cout(outputCounterPrefix
                         << "End of simulation at time "
                         << dTime << " after "
                         << lStep << " steps;" << std::endl
                         << "output in file \"" << sOutputFileName << "\"" << std::endl
                         << "total iterations: " << iTotIter << std::endl
                         << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
                         << "total error: " << dTotErr << std::endl);
 
                 if (pRTSolver) {
                         pRTSolver->Log();
                 }
 
                 return false;
 
         } else if (pRTSolver && pRTSolver->IsStopCommanded()) {
                 silent_cout(outputCounterPrefix
                         << "Simulation is stopped by RTAI task" << std::endl
                         << "Simulation ended at time "
                         << dTime << " after "
                         << lStep << " steps;" << std::endl
                         << "total iterations: " << iTotIter << std::endl
                         << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
                         << "total error: " << dTotErr << std::endl);
                 pRTSolver->Log();
                 return false;
 
         } else if (mbdyn_stop_at_end_of_time_step()
 #ifdef USE_MPI
                 || (MPI_Finalized(&mpi_finalize), mpi_finalize)
 #endif /* USE_MPI */
                         )
         {
                 if (pRTSolver) {
                         pRTSolver->StopCommanded();
                 }
 
                 silent_cout(outputCounterPrefix
                         << "Interrupted!" << std::endl
                         << "Simulation ended at time "
                         << dTime << " after "
                         << lStep << " steps;" << std::endl
                         << "total iterations: " << iTotIter << std::endl
                         << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
                         << "total error: " << dTotErr << std::endl);
 
                 if (pRTSolver) {
                         pRTSolver->Log();
                 }
 
                 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
         }
 
         lStep++;
         pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
 
         Flip();
 
         if (pRTSolver) {
                 pRTSolver->Wait();
         }
 
         int retries = -1;
 IfStepIsToBeRepeated:
         try {
                 retries++;
                 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, lStep);
                 if (outputStep()) {
                         if (outputCounter()) {
                                 silent_cout(std::endl);
                         }
                         silent_cout("Step(" << lStep << ':' << retries << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
                 }
                 dTest = pRegularSteps->Advance(this, dRefTimeStep,
                                 dCurrTimeStep/dRefTimeStep, CurrStep,
                                 qX, qXPrime, pX, pXPrime, iStIter,
                                 dTest, dSolTest);
         }
         catch (NonlinearSolver::NoConvergence) {
                 if (dCurrTimeStep > dMinTimeStep) {
                         /* Riduce il passo */
                         CurrStep = StepIntegrator::REPEATSTEP;
                         doublereal dOldCurrTimeStep = dCurrTimeStep;
                         dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
                         if (dCurrTimeStep < dOldCurrTimeStep) {
                                 DEBUGCOUT("Changing time step"
                                         " from " << dOldCurrTimeStep
                                         << " to " << dCurrTimeStep
                                         << " during step "
                                         << lStep << " after "
                                         << iStIter << " iterations"
                                         << std::endl);
                                 goto IfStepIsToBeRepeated;
                         }
                 }
 
                 silent_cerr(outputCounterPrefix
                         << "Max iterations number "
                         << std::abs(pRegularSteps->GetIntegratorMaxIters())
                         << " has been reached during "
                         "Step=" << lStep << ", "
                         "Time=" << dTime + dCurrTimeStep << "; "
                         "TimeStep=" << dCurrTimeStep
                         << " cannot be reduced further; "
                         "aborting..." << std::endl);
                 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ErrSimulationDiverged) {
                 if (dCurrTimeStep > dMinTimeStep) {
                         /* Riduce il passo */
                         CurrStep = StepIntegrator::REPEATSTEP;
                         doublereal dOldCurrTimeStep = dCurrTimeStep;
                         dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
                         if (dCurrTimeStep < dOldCurrTimeStep) {
                                 DEBUGCOUT("Changing time step"
                                         " from " << dOldCurrTimeStep
                                         << " to " << dCurrTimeStep
                                         << " during step "
                                         << lStep << " after "
                                         << iStIter << " iterations"
                                         << std::endl);
                                 goto IfStepIsToBeRepeated;
                         }
                 }
 
                 silent_cerr(outputCounterPrefix
                         << "Simulation diverged after "
                         << iStIter << " iterations, before "
                         "reaching max iteration number "
                         << std::abs(pRegularSteps->GetIntegratorMaxIters())
                         << " during Step=" << lStep << ", "
                         "Time=" << dTime + dCurrTimeStep << "; "
                         "TimeStep=" << dCurrTimeStep
                         << " cannot be reduced further; "
                         "aborting..." << std::endl);
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (LinearSolver::ErrFactor& err) {
                 /*
                  * Mettere qui eventuali azioni speciali
                  * da intraprendere in caso di errore ...
                  */
                 silent_cerr(outputCounterPrefix
                         << "Simulation failed because no pivot element "
                         "could be found for column " << err.iCol
                         << " (" << pDM->GetDofDescription(err.iCol) << ") "
                         "after " << iStIter << " iterations "
                         "during step " << lStep << "; "
                         "aborting..." << std::endl);
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ConvergenceOnSolution) {
                 bSolConv = true;
         }
         catch (EndOfSimulation& eos) {
                 silent_cerr(outputCounterPrefix
                         << "Simulation ended during a regular step:\n"
                         << eos.what() << "\n");
 #ifdef USE_MPI
                 MBDynComm.Abort(0);
 #endif /* USE_MPI */
                 if (pRTSolver) {
                         pRTSolver->StopCommanded();
                 }
 
                 silent_cout("Simulation ended at time "
                         << dTime << " after "
                         << lStep << " steps;" << std::endl
                         << "total iterations: " << iTotIter << std::endl
                         << "total Jacobian matrices: " << pNLS->TotalAssembledJacobian() << std::endl
                         << "total error: " << dTotErr << std::endl);
 
                 if (pRTSolver) {
                         pRTSolver->Log();
                 }
 
                 return false;
         }
 
         dTotErr += dTest;
         iTotIter += iStIter;
 
         bOut = pDM->Output(lStep, dTime + dCurrTimeStep, dCurrTimeStep);
 
         /* Si fa dare l'std::ostream al file di output per il log */
         std::ostream& Out = pDM->GetOutFile();
 
         if (outputMsg()) {
                 Out << "Step " << lStep
                         << " " << dTime + dCurrTimeStep
                         << " " << dCurrTimeStep
                         << " " << iStIter
                         << " " << dTest
                         << " " << dSolTest
                         << " " << bSolConv
                         << " " << bOut
                         << std::endl;
         }
 
         if (bOutputCounter) {
                 silent_cout("Step " << std::setw(5) << lStep
                         << " " << std::setw(13) << dTime + dCurrTimeStep
                         << " " << std::setw(13) << dCurrTimeStep
                         << " " << std::setw(4) << iStIter
                         << " " << std::setw(13) << dTest
                         << " " << std::setw(13) << dSolTest
                         << " " << bSolConv
                         << " " << bOut
                         << outputCounterPostfix);
         }
 
         DEBUGCOUT("Step " << lStep
                 << " has been successfully completed "
                 "in " << iStIter << " iterations" << std::endl);
 
         dRefTimeStep = dCurrTimeStep;
         dTime += dRefTimeStep;
 
         bSolConv = false;
 
         if (EigAn.bAnalysis
                 && EigAn.currAnalysis != EigAn.Analyses.end()
                 && *EigAn.currAnalysis <= dTime)
         {
                 std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
                         EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
                 if (i != EigAn.Analyses.end()) {
                         EigAn.currAnalysis = --i;
                 }
                 Eig(bOutputCounter);
                 ++EigAn.currAnalysis;
         }
 
         /* Calcola il nuovo timestep */
         dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
         DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
 
         return true;
 }

Here is the call graph for this function:

NonlinearSolver *const Solver::AllocateNonlinearSolver ( )

protected

Definition at line 5073 of file solver.cc.

References MatrixFreeSolver::BICGSTAB, bKeepJac, bTrueNewtonRaphson, dIterertiveEtaMax, dIterertiveTau, dIterTol, MatrixFreeSolver::GMRES, iIterationsBeforeAssembly, iIterativeMaxSteps, iPrecondSteps, NonlinearSolver::LINESEARCH, LineSearch, NonlinearSolver::MATRIXFREE, MFSolverType, NonlinearSolver::NEWTONRAPHSON, NonlinearSolverType, PcType, pNLS, and SAFENEWWITHCONSTRUCTOR.

Referenced by InverseSolver::Prepare(), and Prepare().

 {
         NonlinearSolver *pNLS = 0;
 
         switch (NonlinearSolverType) {
         case NonlinearSolver::MATRIXFREE:
                 switch (MFSolverType) {
                 case MatrixFreeSolver::BICGSTAB:
                         SAFENEWWITHCONSTRUCTOR(pNLS,
                                         BiCGStab,
                                         BiCGStab(PcType,
                                                 iPrecondSteps,
                                                 dIterTol,
                                                 iIterativeMaxSteps,
                                                 dIterertiveEtaMax,
                                                 dIterertiveTau,
                                                 *this));
                         break;
 
                 default:
                         pedantic_cout("unknown matrix free solver type; "
                                         "using default" << std::endl);
                         /* warning: should be unreachable */
 
                 case MatrixFreeSolver::GMRES:
                         SAFENEWWITHCONSTRUCTOR(pNLS,
                                         Gmres,
                                         Gmres(PcType,
                                                 iPrecondSteps,
                                                 dIterTol,
                                                 iIterativeMaxSteps,
                                                 dIterertiveEtaMax,
                                                 dIterertiveTau,
                                                 *this));
                         break;
                 }
                 break;
 
         default:
                 pedantic_cout("unknown nonlinear solver type; using default"
                                 << std::endl);
 
         case NonlinearSolver::NEWTONRAPHSON:
                 SAFENEWWITHCONSTRUCTOR(pNLS,
                                 NewtonRaphsonSolver,
                                 NewtonRaphsonSolver(bTrueNewtonRaphson,
                                         bKeepJac,
                                         iIterationsBeforeAssembly,
                                         *this));
                 break;
         case NonlinearSolver::LINESEARCH:
                 SAFENEWWITHCONSTRUCTOR(pNLS,
                                 LineSearchSolver,
                                 LineSearchSolver(pDM,
                                  *this,
                                  LineSearch));
                 break;
         }
         return pNLS;
 }

SolutionManager *const Solver::AllocateSchurSolman ( integer iStates )

protected

Definition at line 5036 of file solver.cc.

References CurrIntSolver, LinSol::GetSolver(), LinSol::GetSolverName(), iNumDofs, iNumIntDofs, iNumLocDofs, LinSol::LAPACK_SOLVER, MBDYN_EXCEPT_ARGS, LinSol::MESCHACH_SOLVER, LinSol::NAIVE_SOLVER, pIntDofs, pLocalSM, pLocDofs, SAFENEWWITHCONSTRUCTOR, LinSol::UMFPACK_SOLVER, and LinSol::Y12_SOLVER.

Referenced by SetupSolmans().

 {
         SolutionManager *pSSM(0);
 
 #ifdef USE_MPI
         switch (CurrIntSolver.GetSolver()) {
         case LinSol::LAPACK_SOLVER:
         case LinSol::MESCHACH_SOLVER:
         case LinSol::NAIVE_SOLVER:
         case LinSol::UMFPACK_SOLVER:
         case LinSol::Y12_SOLVER:
                 break;
 
         default:
                 silent_cerr("apparently solver "
                                 << CurrIntSolver.GetSolverName()
                                 << " is not allowed as interface solver "
                                 "for SchurSolutionManager" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         SAFENEWWITHCONSTRUCTOR(pSSM,
                         SchurSolutionManager,
                         SchurSolutionManager(iNumDofs, iStates, pLocDofs,
                                 iNumLocDofs,
                                 pIntDofs, iNumIntDofs,
                                 pLocalSM, CurrIntSolver));
 
 #else /* !USE_MPI */
         silent_cerr("Configure --with-mpi to enable Schur solver" << std::endl);
         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif /* !USE_MPI */
 
         return pSSM;
 };

Here is the call graph for this function:

SolutionManager *const Solver::AllocateSolman	(	integer	iNLD,
		integer	iLWS = `0`
	)

protected

Definition at line 5001 of file solver.cc.

References CurrLinearSolver, LinSol::GetSolutionManager(), LinSol::GetSolver(), pRTSolver, and LinSol::UMFPACK_SOLVER.

Referenced by SetupSolmans().

 {
         SolutionManager *pCurrSM = CurrLinearSolver.GetSolutionManager(iNLD, iLWS);
 
         /* special extra parameters if required */
         switch (CurrLinearSolver.GetSolver()) {
         case LinSol::UMFPACK_SOLVER:
 #ifdef HAVE_UMFPACK_TIC_DISABLE
                 if (pRTSolver) {
                         /* disable profiling, to avoid times() system call
                          *
                          * This function has been introduced in Umfpack 4.1
                          * by our patch at
                          *
                          * http://mbdyn.aero.polimi.it/~masarati/Download/\
                          *      mbdyn/umfpack-4.1-nosyscalls.patch
                          *
                          * but since Umfpack 4.3 is no longer required,
                          * provided the library is compiled with -DNO_TIMER
                          * to disable run-time syscalls to timing routines.
                          */
                         umfpack_tic_disable();
                 }
 #endif // HAVE_UMFPACK_TIC_DISABLE
                 break;
 
         default:
                 break;
         }
 
         return pCurrSM;
 };

Here is the call graph for this function:

void Solver::CheckTimeStepLimit	(	doublereal	dErr,
		doublereal	dErrDiff
	)		const
throw	(	NonlinearSolver::TimeStepLimitExceeded,
		NonlinearSolver::MaxResidualExceeded
	)

virtual

Definition at line 5205 of file solver.cc.

References ASSERT, and MBDYN_EXCEPT_ARGS.

Referenced by LineSearchSolver::LineSearch(), BiCGStab::Solve(), Gmres::Solve(), NewtonRaphsonSolver::Solve(), and LineSearchSolver::Solve().

 {
         if (pDerivativeSteps) {
                 // Time step cannot be reduced
                 return;
         }
 
         if (dErr > dMaxResidual) {
             if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
                 if (outputIters()) {
 #ifdef USE_MPI
                         if (!bParallel || MBDynComm.Get_rank() == 0)
 #endif /* USE_MPI */
                         {
                                 silent_cerr("warning: current residual = " << dErr
                                                 << " > maximum residual = " << dMaxResidual
                                                 << std::endl);
                         }
                 }
 
                 throw NonlinearSolver::MaxResidualExceeded(MBDYN_EXCEPT_ARGS);
         }
         }
 
         if (dErrDiff > dMaxResidualDiff) {
             if (dCurrTimeStep > 2 * dMinTimeStep) { // FIXME
                 if (outputIters()) {
 #ifdef USE_MPI
                         if (!bParallel || MBDynComm.Get_rank() == 0)
 #endif /* USE_MPI */
                         {
                                 silent_cerr("warning: current residual for differential equations = " << dErrDiff
                                                 << " > maximum residual = " << dMaxResidualDiff
                                                 << std::endl);
                         }
                 }
 
                 throw NonlinearSolver::MaxResidualExceeded(MBDYN_EXCEPT_ARGS);
         }
         }
 
         switch (eTimeStepLimit) {
         case TS_SOFT_LIMIT:
                 break;
 
         case TS_HARD_LIMIT: {
                 const doublereal dMaxTS = MaxTimeStep.dGet();
 
                 if (dCurrTimeStep > dMaxTS && dCurrTimeStep > dMinTimeStep) {
                         if (outputIters()) {
 #ifdef USE_MPI
                                 if (!bParallel || MBDynComm.Get_rank() == 0)
 #endif /* USE_MPI */
                                 {
                                         silent_cerr("warning: current time step = "
                                                 << dCurrTimeStep
                                                 << " > hard limit of the maximum time step = "
                                                 << dMaxTS << std::endl);
                                 }
                         }
 
                         throw NonlinearSolver::TimeStepLimitExceeded(MBDYN_EXCEPT_ARGS);
                 }
                 } break;
 
         default:
                 ASSERT(0);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 }

virtual doublereal Solver::dGetFinalTime ( void ) const

inlinevirtual

Definition at line 414 of file solver.h.

References dFinalTime.

                                                      {
                 return dFinalTime;
         };

doublereal Solver::dGetInitialMaxTimeStep ( ) const

protected

Definition at line 4990 of file solver.cc.

References dDefaultMaxTimeStep, DriveOwner::dGet(), MaxTimeStep, and DriveOwner::pGetDriveCaller().

Referenced by InverseSolver::ReadData(), and ReadData().

 {
         if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
                 return ::dDefaultMaxTimeStep;
         }
 
         // The same behavior like in previous releases
         return MaxTimeStep.dGet();
 }

Here is the call graph for this function:

virtual doublereal Solver::dGetInitialTimeStep ( void ) const

inlinevirtual

Definition at line 417 of file solver.h.

References dInitialTimeStep.

Referenced by DataManager::AfterConvergence(), and ModuleFMU::ModuleFMU().

                                                            {
                 return dInitialTimeStep;
         };

void Solver::Eig ( bool bNewLine = false )

protected

Definition at line 4853 of file solver.cc.

Referenced by Advance(), Prepare(), and Start().

 {
         DEBUGCOUTFNAME("Solver::Eig");
 
         /*
          * MatA, MatB: MatrixHandlers to eigenanalysis matrices
          * MatVL, MatVR: MatrixHandlers to eigenvectors, if required
          * AlphaR, AlphaI Beta: eigenvalues
          * WorkVec:    Workspace
          * iWorkSize:  Size of the workspace
          */
 
         DEBUGCOUT("Solver::Eig(): performing eigenanalysis" << std::endl);
 
         integer iSize = iNumDofs;
 
         SolutionManager *pSM = 0;
         MatrixHandler *pMatA = 0;
         MatrixHandler *pMatB = 0;
 
         if (EigAn.uFlags & EigenAnalysis::EIG_USE_LAPACK) {
                 SAFENEWWITHCONSTRUCTOR(pMatA, FullMatrixHandler,
                         FullMatrixHandler(iSize));
 
                 SAFENEWWITHCONSTRUCTOR(pMatB, FullMatrixHandler,
                         FullMatrixHandler(iSize));
 
         } else if (EigAn.uFlags & EigenAnalysis::EIG_USE_ARPACK) {
                 SAFENEWWITHCONSTRUCTOR(pSM, NaiveSparsePermSolutionManager<Colamd_ordering>,
                         NaiveSparsePermSolutionManager<Colamd_ordering>(iSize));
                 SAFENEWWITHCONSTRUCTOR(pMatA, NaiveMatrixHandler,
                         NaiveMatrixHandler(iSize));
                 pMatB = pSM->pMatHdl();
 
         } else if (EigAn.uFlags & EigenAnalysis::EIG_USE_JDQZ) {
                 SAFENEWWITHCONSTRUCTOR(pMatA, NaiveMatrixHandler,
                         NaiveMatrixHandler(iSize));
                 SAFENEWWITHCONSTRUCTOR(pMatB, NaiveMatrixHandler,
                         NaiveMatrixHandler(iSize));
 
         } else if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
                 SAFENEWWITHCONSTRUCTOR(pMatA, SpMapMatrixHandler,
                         SpMapMatrixHandler(iSize));
                 SAFENEWWITHCONSTRUCTOR(pMatB, SpMapMatrixHandler,
                         SpMapMatrixHandler(iSize));
         }
 
         pMatA->Reset();
         pMatB->Reset();
 
         // Matrices assembly (see eig.ps)
         doublereal h = EigAn.dParam;
         pDM->AssJac(*pMatA, -h/2.);
         pDM->AssJac(*pMatB, h/2.);
 
 #ifdef DEBUG
         DEBUGCOUT(std::endl
                 << "Matrix A:" << std::endl << *pMatA << std::endl
                 << "Matrix B:" << std::endl << *pMatB << std::endl);
 #endif /* DEBUG */
 
         unsigned uCurr = EigAn.currAnalysis - EigAn.Analyses.begin();
         if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT) {
                 unsigned uSize = EigAn.Analyses.size();
                 if (uSize >= 1) {
                         std::stringstream postfix_ss;
 
                         postfix_ss << '_';
 
                         if (EigAn.iFNameWidth > 0) {
                                 postfix_ss << std::setw(EigAn.iFNameWidth) << std::setfill('0') << uCurr;
 
                         } else if (!EigAn.iFNameFormat.empty()) {
                                 char buf[BUFSIZ];
 
                                 snprintf(buf, sizeof(buf), EigAn.iFNameFormat.c_str(), (int)uCurr);
                                 postfix_ss << buf;
 
                         } else {
                                 postfix_ss << uCurr;
                         }
                         
                         pDM->OutputEigOpen(postfix_ss.str());
                         pDM->OutputEigParams(dTime, h/2., uCurr, EigAn.iResultsPrecision);
                 }
         }
         if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
                 pDM->OutputEigFullMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
         }
         else if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
                 if (dynamic_cast<const NaiveMatrixHandler *>(pMatB)) {
                         pDM->OutputEigNaiveMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
                 } else {
                         pDM->OutputEigSparseMatrices(pMatA, pMatB, uCurr, EigAn.iMatrixPrecision);
                 }
         }
 
         switch (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
 #ifdef USE_LAPACK
         case EigenAnalysis::EIG_USE_LAPACK:
                 eig_lapack(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
                 break;
 #endif // USE_LAPACK
 
 #ifdef USE_ARPACK
         case EigenAnalysis::EIG_USE_ARPACK:
                 eig_arpack(pMatA, pSM, pDM, &EigAn, bNewLine, uCurr);
                 break;
 #endif // USE_ARPACK
 
 #ifdef USE_JDQZ
         case EigenAnalysis::EIG_USE_JDQZ:
                 eig_jdqz(pMatA, pMatB, pDM, &EigAn, bNewLine, uCurr);
                 break;
 #endif // USE_JDQZ
 
         default:
                 // only output matrices, use external eigenanalysis
                 break;
         }
 
         pDM->OutputEigClose();
 
         if (pSM) {
                 pMatB = 0;
                 SAFEDELETE(pSM);
         }
 
         if (pMatA) {
                 SAFEDELETE(pMatA);
         }
 
         if (pMatB) {
                 SAFEDELETE(pMatB);
         }
 }

Here is the call graph for this function:

void Solver::Flip ( void )

inlineprotected

Definition at line 429 of file solver.h.

References iNumDofs, MyVectorHandler::PutCoef(), pX, pXPrime, qX, and qXPrime.

Referenced by Advance(), Prepare(), and Start().

 {
         /*
          * switcha i puntatori; in questo modo non e' necessario
          * copiare i vettori per cambiare passo
          */
         qX.push_front(qX.back());
         qX.pop_back();
         qXPrime.push_front(qXPrime.back());
         qXPrime.pop_back();
 
         /* copy from pX, pXPrime to qx[0], qxPrime[0] */
         MyVectorHandler* x = qX[0];
         MyVectorHandler* xp = qXPrime[0];
         for (integer i = 1; i <= iNumDofs; i++) {
                 x->PutCoef(i, pX->operator()(i));
                 xp->PutCoef(i, pXPrime->operator()(i));
         }
 }

Here is the call graph for this function:

clock_t Solver::GetCPUTime ( void ) const

virtual

Definition at line 5188 of file solver.cc.

References DataManager::GetCPUTime().

Referenced by mbdyn_program().

 {
         return pDM->GetCPUTime();
 }

Here is the call graph for this function:

virtual const LinSol& Solver::GetLinearSolver ( void ) const

inlinevirtual

Definition at line 401 of file solver.h.

References CurrLinearSolver.

                                                           {
                 return CurrLinearSolver;
         };

virtual DataManager* Solver::pGetDataManager ( void ) const

inlinevirtual

Definition at line 395 of file solver.h.

References pDM.

Referenced by NonlinearSolver::MakeResTest(), NonlinearSolver::MakeSolTest(), NonlinearSolverTest::MakeTest(), mb_sol_setbufin(), and mb_sol_setbufout().

                                                          {
                 return pDM;
         };

virtual NonlinearSolver* Solver::pGetNonlinearSolver ( void ) const

inlinevirtual

Definition at line 404 of file solver.h.

References pNLS.

Referenced by ThirdOrderIntegrator::Advance(), DerivativeSolver::Advance(), Step1Integrator::Advance(), Step2Integrator::Advance(), and InverseDynamicsStepSolver::Advance().

                                                                  {
                 return pNLS;
         };

virtual SolutionManager* Solver::pGetSolutionManager ( void ) const

inlinevirtual

Definition at line 398 of file solver.h.

References pSM.

Referenced by InverseDynamicsStepSolver::Advance(), LineSearchSolver::Jacobian(), NonlinearSolverTest::MakeTest(), BiCGStab::Solve(), NewtonRaphsonSolver::Solve(), Gmres::Solve(), and LineSearchSolver::Solve().

                                                                  {
                 return pSM;
         };

virtual StepIntegrator* Solver::pGetStepIntegrator ( void ) const

inlinevirtual

Definition at line 407 of file solver.h.

References ASSERT, and pCurrStepIntegrator.

Referenced by DataManager::AssConstrJac().

                                                                {
                 // if (pCurrStepIntegrator == NULL)
                 //      return pRegularSteps;
                 ASSERT(pCurrStepIntegrator != 0);
 
                 return pCurrStepIntegrator;
         };

bool Solver::Prepare ( void )

virtual

Reimplemented in InverseSolver.

Definition at line 399 of file solver.cc.

Referenced by mb_sol_prepare(), and Run().

 {
         DEBUGCOUTFNAME("Solver::Prepare");
 
         // consistency check
         if (eStatus != SOLVER_STATUS_UNINITIALIZED) {
                 silent_cerr("Prepare() must be called first" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         mbdyn_signal_init(1);
 
         /* Legge i dati relativi al metodo di integrazione */
         ReadData(HP);
 
 #ifdef USE_MULTITHREAD
         ThreadPrepare();
 #endif /* USE_MULTITHREAD */
 
         if (pRTSolver) {
                 pRTSolver->Setup();
         }
 
 #ifdef USE_MPI
         int mpi_finalize = 0;
 #endif
 
 #ifdef USE_SCHUR
         if (bParallel) {
                 DEBUGLCOUT(MYDEBUG_MEM, "creating parallel SchurDataManager"
                                 << std::endl);
 
                 SAFENEWWITHCONSTRUCTOR(pSDM,
                         SchurDataManager,
                         SchurDataManager(HP,
                                 OutputFlags,
                                 this,
                                 dInitialTime,
                                 sOutputFileName.c_str(),
                                 sInputFileName.c_str(),
                                 eAbortAfter == AFTER_INPUT));
 
                 pDM = pSDM;
 
         } else
 #endif // USE_SCHUR
         {
                 /* chiama il gestore dei dati generali della simulazione */
 #ifdef USE_MULTITHREAD
                 if (nThreads > 1) {
                         if (!(CurrLinearSolver.GetSolverFlags() & LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS)) {
                                 /* conservative: dir may use too much memory */
                                 if (!CurrLinearSolver.AddSolverFlags(LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS)) {
                                         bool b;
 
 #if defined(USE_UMFPACK)
                                         b = CurrLinearSolver.SetSolver(LinSol::UMFPACK_SOLVER,
                                                         LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS);
 #elif defined(USE_Y12)
                                         b = CurrLinearSolver.SetSolver(LinSol::Y12_SOLVER,
                                                         LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS);
 #else
                                         b = false;
 #endif
                                         if (!b) {
                                                 silent_cerr("unable to select a CC-capable solver"
                                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                 }
                         }
 
                         silent_cout("Creating multithread solver "
                                         "with " << nThreads << " threads "
                                         "and "
                                         << CurrLinearSolver.GetSolverName()
                                         << " linear solver"
                                         << std::endl);
 
                         SAFENEWWITHCONSTRUCTOR(pDM,
                                         MultiThreadDataManager,
                                         MultiThreadDataManager(HP,
                                                 OutputFlags,
                                                 this,
                                                 dInitialTime,
                                                 sOutputFileName.c_str(),
                                                 sInputFileName.c_str(),
                                                 eAbortAfter == AFTER_INPUT,
                                                 nThreads));
 
                 } else
 #endif /* USE_MULTITHREAD */
                 {
                         DEBUGLCOUT(MYDEBUG_MEM, "creating DataManager"
                                         << std::endl);
 
                         silent_cout("Creating scalar solver "
                                         "with "
                                         << CurrLinearSolver.GetSolverName()
                                         << " linear solver"
                                         << std::endl);
 
                         SAFENEWWITHCONSTRUCTOR(pDM,
                                         DataManager,
                                         DataManager(HP,
                                                 OutputFlags,
                                                 this,
                                                 dInitialTime,
                                                 sOutputFileName.c_str(),
                                                 sInputFileName.c_str(),
                                                 eAbortAfter == AFTER_INPUT));
                 }
         }
 
         // log symbol table
         std::ostream& log = pDM->GetLogFile();
         log << "Symbol table:" << std::endl;
         log << HP.GetMathParser().GetSymbolTable();
 
 #ifdef HAVE_ENVIRON
         // log environment
         log << "Environment:" << std::endl;
         for (int i = 0; environ[i] != NULL; i++) {
                 log << "  " << environ[i] << std::endl;
         }
 #endif // HAVE_ENVIRON
 
         // close input stream
         HP.Close();
 
         /* Si fa dare il DriveHandler e linka i drivers di rho ecc. */
         const DriveHandler* pDH = pDM->pGetDrvHdl();
         SetOutputDriveHandler(pDH);
 
         bOutputCounter = outputCounter() && isatty(fileno(stderr));
         outputCounterPrefix = bOutputCounter ? "\n" : "";
         outputCounterPostfix = outputStep() ? "\n" : "\r";
 
         /* Si fa dare l'std::ostream al file di output per il log */
         std::ostream& Out = pDM->GetOutFile();
 
         if (eAbortAfter == AFTER_INPUT) {
                 /* Esce */
                 pDM->Output(0, dTime, 0., true);
                 Out << "End of Input; no simulation or assembly is required."
                         << std::endl;
                 return false;
 
         } else if (eAbortAfter == AFTER_ASSEMBLY) {
                 /* Fa l'output dell'assemblaggio iniziale e poi esce */
                 pDM->Output(0, dTime, 0., true);
                 Out << "End of Initial Assembly; no simulation is required."
                         << std::endl;
                 return false;
         }
 
 #ifdef USE_SCHUR
         /* Qui crea le partizioni: principale fra i processi, se parallelo  */
         if (bParallel) {
                 pSDM->CreatePartition();
         }
 #endif // USE_SCHUR
 
         pRegularSteps->SetDriveHandler(pDH);
         if (iDummyStepsNumber) {
                 pDummySteps->SetDriveHandler(pDH);
         }
 
         /* Costruisce i vettori della soluzione ai vari passi */
         DEBUGLCOUT(MYDEBUG_MEM, "creating solution vectors" << std::endl);
 
 #ifdef USE_SCHUR
         if (bParallel) {
                 iNumDofs = pSDM->HowManyDofs(SchurDataManager::TOTAL);
                 pDofs = pSDM->pGetDofsList();
 
                 iNumLocDofs = pSDM->HowManyDofs(SchurDataManager::LOCAL);
                 pLocDofs = pSDM->GetDofsList(SchurDataManager::LOCAL);
 
                 iNumIntDofs = pSDM->HowManyDofs(SchurDataManager::INTERNAL);
                 pIntDofs = pSDM->GetDofsList(SchurDataManager::INTERNAL);
 
         } else
 #endif // USE_SCHUR
         {
                 iNumDofs = pDM->iGetNumDofs();
         }
 
         /* relink those known drive callers that might need
          * the data manager, but were verated ahead of it */
         pTSC->SetDriveHandler(pDM->pGetDrvHdl());
 
         ASSERT(iNumDofs > 0);
 
         integer iRSteps = pRegularSteps->GetIntegratorNumPreviousStates();
         integer iFSteps = 0;
         if (iDummyStepsNumber) {
                 iFSteps = pDummySteps->GetIntegratorNumPreviousStates();
         }
         iNumPreviousVectors = (iRSteps < iFSteps) ? iFSteps : iRSteps;
 
         integer iRUnkStates = pRegularSteps->GetIntegratorNumUnknownStates();
         integer iFUnkStates = 0;
         if (iDummyStepsNumber) {
                 iFUnkStates = pDummySteps->GetIntegratorNumUnknownStates();
         }
         iUnkStates = (iRUnkStates < iFUnkStates) ? iFUnkStates : iRUnkStates;
 
         /* allocate workspace for previous time steps */
         SAFENEWARR(pdWorkSpace, doublereal,
                 2*(iNumPreviousVectors)*iNumDofs);
         /* allocate MyVectorHandlers for previous time steps: use workspace */
         for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
                 SAFENEWWITHCONSTRUCTOR(pX,
                         MyVectorHandler,
                         MyVectorHandler(iNumDofs, pdWorkSpace+ivec*iNumDofs));
                 qX.push_back(pX);
                 SAFENEWWITHCONSTRUCTOR(pXPrime,
                         MyVectorHandler,
                         MyVectorHandler(iNumDofs,
                                 pdWorkSpace+((iNumPreviousVectors)+ivec)*iNumDofs));
                 qXPrime.push_back(pXPrime);
                 pX = 0;
                 pXPrime = 0;
         }
         /* allocate MyVectorHandlers for unknown time step(s): own memory */
         SAFENEWWITHCONSTRUCTOR(pX,
                 MyVectorHandler,
                 MyVectorHandler(iUnkStates*iNumDofs));
         SAFENEWWITHCONSTRUCTOR(pXPrime,
                 MyVectorHandler,
                 MyVectorHandler(iUnkStates*iNumDofs));
 
 
         /* Resetta i vettori */
         for (int ivec = 0; ivec < iNumPreviousVectors; ivec++) {
                 qX[ivec]->Reset();
                 qXPrime[ivec]->Reset();
         }
         pX->Reset();
         pXPrime->Reset();
 
         /*
          * Immediately link DataManager to current solution
          *
          * this should work as long as the last unknown time step is put
          * at the beginning of pX, pXPrime
          */
         pDM->LinkToSolution(*pX, *pXPrime);
 
         /* a questo punto si costruisce il nonlinear solver */
         pNLS = AllocateNonlinearSolver();
 
         Scale.ResizeReset(iNumDofs);
         if (bScale) {
                 /* collects scale factors from data manager */
                 pDM->SetScale(Scale);
         }
 
         /*
          * prepare tests for nonlinear solver;
          *
          * test on residual may allow pre-scaling;
          * test on solution (difference between two iterations) does not
          */
         NonlinearSolverTest *pResTest = 0;
         if (bScale) {
                 NonlinearSolverTestScale *pResTestScale = 0;
 
                 switch (ResTest) {
                 case NonlinearSolverTest::NORM:
                         SAFENEW(pResTestScale, NonlinearSolverTestScaleNorm);
                         break;
 
                 case NonlinearSolverTest::MINMAX:
                         SAFENEW(pResTestScale, NonlinearSolverTestScaleMinMax);
                         break;
 
                 default:
                         ASSERT(0);
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
 
                 /* registers scale factors at nonlinear solver */
                 pResTestScale->SetScale(&Scale);
 
                 pResTest = pResTestScale;
 
 
         } else {
                 switch (ResTest) {
                 case NonlinearSolverTest::NONE:
                         SAFENEW(pResTest, NonlinearSolverTestNone);
                         break;
 
                 case NonlinearSolverTest::NORM:
                         SAFENEW(pResTest, NonlinearSolverTestNorm);
                         break;
 
                 case NonlinearSolverTest::MINMAX:
                         SAFENEW(pResTest, NonlinearSolverTestMinMax);
                         break;
 
                 default:
                         ASSERT(0);
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
         }
 
         NonlinearSolverTest *pSolTest = 0;
         switch (SolTest) {
         case NonlinearSolverTest::NONE:
                 SAFENEW(pSolTest, NonlinearSolverTestNone);
                 break;
 
         case NonlinearSolverTest::NORM:
                 SAFENEW(pSolTest, NonlinearSolverTestNorm);
                 break;
 
         case NonlinearSolverTest::MINMAX:
                 SAFENEW(pSolTest, NonlinearSolverTestMinMax);
                 break;
 
         default:
                 ASSERT(0);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         /* registers tests in nonlinear solver */
         pNLS->SetTest(pResTest, pSolTest);
 
         /*
          * Dell'assemblaggio iniziale dei vincoli se ne occupa il DataManager
          * in quanto e' lui il responsabile dei dati della simulazione,
          * e quindi anche della loro coerenza. Inoltre e' lui a sapere
          * quali equazioni sono di vincolo o meno.
          */
 
         pDM->SetValue(*pX, *pXPrime);
 
 
         /*
          * Prepare output
          */
         pDM->OutputPrepare();
 
         /*
          * If eigenanalysis is requested, prepare output for it
          */
         if (EigAn.bAnalysis) {
                 pDM->OutputEigPrepare(EigAn.Analyses.size(), iNumDofs);
                 if (EigAn.iFNameWidth == EigenAnalysis::EIGAN_WIDTH_COMPUTE) {
                         EigAn.iFNameWidth = int(std::log10(double(EigAn.Analyses.size()))) + 1;
                 }
         }
 
         /*
          * Dialoga con il DataManager per dargli il tempo iniziale
          * e per farsi inizializzare i vettori di soluzione e derivata
          */
         /* FIXME: the time is already set by DataManager, but FileDrivers
          * have not been ServePending'd
          */
         dTime = dInitialTime;
         pDM->SetTime(dTime, dInitialTimeStep, 0);
 
         EigAn.currAnalysis = std::find_if(EigAn.Analyses.begin(), EigAn.Analyses.end(),
                 bind2nd(std::greater<doublereal>(), dTime));
         if (EigAn.currAnalysis != EigAn.Analyses.end() && EigAn.currAnalysis != EigAn.Analyses.begin()) {
                 --EigAn.currAnalysis;
         }
 
         // if eigenanalysis is requested and currAnalysis points
         // past the end of the array, the analysis was requested
         // at Time < initial time; perform *before* derivatives
         if (EigAn.bAnalysis
                 && ((EigAn.currAnalysis == EigAn.Analyses.end()
                                 && EigAn.Analyses.back() < dTime)
                         || (EigAn.currAnalysis != EigAn.Analyses.end()
                                 && *EigAn.currAnalysis < dTime)))
         {
                 Eig();
                 if (EigAn.currAnalysis != EigAn.Analyses.end()) {
                         ++EigAn.currAnalysis;
                 }
         }
 
         /* calcolo delle derivate */
         DEBUGLCOUT(MYDEBUG_DERIVATIVES, "derivatives solution step"
                         << std::endl);
 
         mbdyn_signal_init(0);
 
         /* settaggio degli output Types */
         unsigned OF = OutputFlags;
         if ( DEBUG_LEVEL_MATCH(MYDEBUG_RESIDUAL) ) {
                 OF |= OUTPUT_RES;
         }
         if ( DEBUG_LEVEL_MATCH(MYDEBUG_JAC) ) {
                 OF |= OUTPUT_JAC;
         }
         if ( DEBUG_LEVEL_MATCH(MYDEBUG_SOL) ) {
                 OF |= OUTPUT_SOL;
         }
         pNLS->SetOutputFlags(OF);
         if (pOutputMeter) {
                 pOutputMeter->SetDrvHdl(pDM->pGetDrvHdl());
                 pNLS->SetOutputMeter(pOutputMeter->pCopy());
         }
 
         pDerivativeSteps->SetDataManager(pDM);
         pDerivativeSteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
         if (iDummyStepsNumber) {
                 pFirstDummyStep->SetDataManager(pDM);
                 pFirstDummyStep->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
                 pDummySteps->SetDataManager(pDM);
                 pDummySteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
         }
         pFirstRegularStep->SetDataManager(pDM);
         pFirstRegularStep->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
         pRegularSteps->SetDataManager(pDM);
         pRegularSteps->OutputTypes(DEBUG_LEVEL_MATCH(MYDEBUG_PRED));
 
 #ifdef USE_EXTERNAL
         pNLS->SetExternal(External::EMPTY);
 #endif /* USE_EXTERNAL */
         /* Setup SolutionManager(s) */
         SetupSolmans(pDerivativeSteps->GetIntegratorNumUnknownStates());
 
         /* Derivative steps */
         pCurrStepIntegrator = pDerivativeSteps;
         try {
                 if (outputStep()) {
                         if (outputCounter()) {
                                 silent_cout(std::endl);
                         }
                         silent_cout("Derivatives t=" << dTime << " coef=" << dDerivativesCoef << std::endl);
                 }
                 dTest = pDerivativeSteps->Advance(this,
                         dInitialTimeStep, 1., StepIntegrator::NEWSTEP,
                         qX, qXPrime, pX, pXPrime,
                         iStIter, dTest, dSolTest);
         }
         catch (NonlinearSolver::NoConvergence) {
                 silent_cerr("Initial derivatives calculation " << iStIter
                         << " does not converge; aborting..." << std::endl
                         << "(hint: try playing with the \"derivatives coefficient\" value)" << std::endl);
                 pDM->Output(0, dTime, 0., true);
                 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ErrSimulationDiverged) {
                 /*
                  * Mettere qui eventuali azioni speciali
                  * da intraprendere in caso di errore ...
                  */
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (LinearSolver::ErrFactor& err) {
                 /*
                  * Mettere qui eventuali azioni speciali
                  * da intraprendere in caso di errore ...
                  */
                 silent_cerr("Initial derivatives failed because no pivot element "
                         "could be found for column " << err.iCol
                         << " (" << pDM->GetDofDescription(err.iCol) << "); "
                         "aborting..." << std::endl);
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ConvergenceOnSolution) {
                 bSolConv = true;
         }
         catch (EndOfSimulation& eos) {
                 silent_cerr("Simulation ended during the derivatives steps:\n" << eos.what() << "\n");
                 return false;
         }
 
         SAFEDELETE(pDerivativeSteps);
         pDerivativeSteps = 0;
 
 #if 0
         /* don't sum up the derivatives error */
         dTotErr  += dTest;
 #endif
         iTotIter += iStIter;
 
         if (outputMsg()) {
                 Out << "# Derivatives solution step at time " << dInitialTime
                         << " performed in " << iStIter
                         << " iterations with " << dTest
                         << " error" << std::endl;
         }
 
         DEBUGCOUT("Derivatives solution step has been performed successfully"
                 " in " << iStIter << " iterations" << std::endl);
 
 #ifdef USE_EXTERNAL
         /* comunica che gli ultimi dati inviati sono la condizione iniziale */
         if (iDummyStepsNumber == 0) {
                 External::SendInitial();
         }
 #endif /* USE_EXTERNAL */
 
         if (eAbortAfter == AFTER_DERIVATIVES) {
                 /*
                  * Fa l'output della soluzione delle derivate iniziali ed esce
                  */
 
 
                 pDM->Output(0, dTime, 0., true);
                 Out << "End of derivatives; no simulation is required."
                         << std::endl;
                 return false;
 
         } else if (mbdyn_stop_at_end_of_time_step()) {
                 /*
                  * Fa l'output della soluzione delle derivate iniziali ed esce
                  */
                 pDM->Output(0, dTime, 0., true);
                 Out << "Interrupted during derivatives computation." << std::endl;
                 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
         }
 
         // if eigenanalysis is requested and currAnalysis points
         // past the end of the array, the analysis was requested
         // at Time == initial time; perform *after* derivatives
         if (EigAn.bAnalysis) {
                 ASSERT(EigAn.Analyses.size() > 0);
 
                 if ((EigAn.currAnalysis == EigAn.Analyses.end()
                                 && EigAn.Analyses.back() == dTime)
                         || (EigAn.currAnalysis != EigAn.Analyses.end()
                                 && *EigAn.currAnalysis == dTime))
                 {
                         Eig();
                         if (EigAn.currAnalysis != EigAn.Analyses.end()) {
                                 ++EigAn.currAnalysis;
                         }
                 }
         }
 
         /* Dati comuni a passi fittizi e normali */
         lStep = 1;
 
         if (iDummyStepsNumber > 0) {
                 /* passi fittizi */
 
                 /*
                  * inizio integrazione: primo passo a predizione lineare
                  * con sottopassi di correzione delle accelerazioni
                  * e delle reazioni vincolari
                  */
                 pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
                 Flip();
 
                 dRefTimeStep = dInitialTimeStep*dDummyStepsRatio;
                 dCurrTimeStep = dRefTimeStep;
                 /* FIXME: do we need to serve pending drives in dummy steps? */
                 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 0);
 
                 DEBUGLCOUT(MYDEBUG_FSTEPS, "Current time step: "
                         << dCurrTimeStep << std::endl);
 
                 if (outputStep()) {
                     silent_cout("Dummy Step(" << lStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
                 }
                 
                 ASSERT(pFirstDummyStep != 0);
 
                 /* Setup SolutionManager(s) */
                 SetupSolmans(pFirstDummyStep->GetIntegratorNumUnknownStates());
                 /* pFirstDummyStep */
                 pCurrStepIntegrator = pFirstDummyStep;
                 try {
                         dTest = pFirstDummyStep->Advance(this,
                                 dRefTimeStep, 1.,
                                 StepIntegrator::NEWSTEP,
                                 qX, qXPrime, pX, pXPrime,
                                 iStIter, dTest, dSolTest);
                 }
                 catch (NonlinearSolver::NoConvergence) {
                         silent_cerr("First dummy step does not converge; "
                                 "TimeStep=" << dCurrTimeStep
                                 << " cannot be reduced further; "
                                 "aborting..." << std::endl);
                         pDM->Output(0, dTime, dCurrTimeStep, true);
                         throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
                 }
                 catch (NonlinearSolver::ErrSimulationDiverged) {
                         /*
                          * Mettere qui eventuali azioni speciali
                          * da intraprendere in caso di errore ...
                          */
                         throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
                 }
                 catch (LinearSolver::ErrFactor& err) {
                         /*
                          * Mettere qui eventuali azioni speciali
                          * da intraprendere in caso di errore ...
                          */
                         silent_cerr("First dummy step failed because no pivot element "
                                 "could be found for column " << err.iCol
                                 << " (" << pDM->GetDofDescription(err.iCol) << "); "
                                 "aborting..." << std::endl);
                         throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
                 }
                 catch (NonlinearSolver::ConvergenceOnSolution) {
                         bSolConv = true;
                 }
                 catch (EndOfSimulation& eos) {
                         silent_cerr("Simulation ended during the first dummy step:\n"
                                 << eos.what() << "\n");
                         return false;
                 }
 
                 SAFEDELETE(pFirstDummyStep);
                 pFirstDummyStep = 0;
 
                 dRefTimeStep = dCurrTimeStep;
                 dTime += dRefTimeStep;
 
 #if 0
                 /* don't sum up the derivatives error */
                 dTotErr += dTest;
 #endif
                 iTotIter += iStIter;
 
                 if (mbdyn_stop_at_end_of_time_step()) {
                         /*
                          * Fa l'output della soluzione delle derivate iniziali
                          * ed esce
                          */
 #ifdef DEBUG_FICTITIOUS
                         pDM->Output(0, dTime, dCurrTimeStep, true);
 #endif /* DEBUG_FICTITIOUS */
                         Out << "Interrupted during first dummy step." << std::endl;
                         throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
                 }
 
 #ifdef DEBUG_FICTITIOUS
                 pDM->Output(0, dTime, dCurrTimeStep, true);
 #endif /* DEBUG_FICTITIOUS */
 
                 /* Passi fittizi successivi */
                 if (iDummyStepsNumber > 1) {
                         /* Setup SolutionManager(s) */
                         SetupSolmans(pDummySteps->GetIntegratorNumUnknownStates());
                 }
 
                 for (int iSubStep = 2;
                         iSubStep <= iDummyStepsNumber;
                         iSubStep++)
                 {
                         pDM->BeforePredict(*pX, *pXPrime,
                                 *qX[0], *qXPrime[0]);
                         Flip();
 
                         DEBUGLCOUT(MYDEBUG_FSTEPS, "Dummy step "
                                 << iSubStep
                                 << "; current time step: " << dCurrTimeStep
                                 << std::endl);
 
                         if (outputStep()) {
                             silent_cout("Dummy Step(" << iSubStep << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
                         }
                         
                         pCurrStepIntegrator = pDummySteps;
                         ASSERT(pDummySteps!= 0);
                         try {
                                 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 0);
                                 dTest = pDummySteps->Advance(this,
                                         dRefTimeStep,
                                         dCurrTimeStep/dRefTimeStep,
                                         StepIntegrator::NEWSTEP,
                                         qX, qXPrime, pX, pXPrime,
                                         iStIter, dTest, dSolTest);
                         }
                         catch (NonlinearSolver::NoConvergence) {
                                 silent_cerr("Dummy step " << iSubStep
                                         << " does not converge; "
                                         "TimeStep=" << dCurrTimeStep
                                         << " cannot be reduced further; "
                                         "aborting..." << std::endl);
                                 pDM->Output(0, dTime, dCurrTimeStep, true);
                                 throw ErrMaxIterations(MBDYN_EXCEPT_ARGS);
                         }
 
                         catch (NonlinearSolver::ErrSimulationDiverged) {
                                 /*
                                  * Mettere qui eventuali azioni speciali
                                  * da intraprendere in caso di errore ...
                                  */
                                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
                         }
                         catch (LinearSolver::ErrFactor& err) {
                                 /*
                                  * Mettere qui eventuali azioni speciali
                                  * da intraprendere in caso di errore ...
                                  */
                                 silent_cerr("Dummy step " << iSubStep
                                         << " failed because no pivot element "
                                         "could be found for column " << err.iCol
                                         << " (" << pDM->GetDofDescription(err.iCol) << "); "
                                         "aborting..." << std::endl);
                                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
                         }
                         catch (NonlinearSolver::ConvergenceOnSolution) {
                                 bSolConv = true;
                         }
                         catch (EndOfSimulation& eos) {
                                 silent_cerr("Simulation ended during the dummy steps:\n"
                                         << eos.what() << "\n");
                                 return false;
                         }
 
 #if 0
                         /* don't sum up the derivatives error */
                         dTotErr += dTest;
 #endif
                         iTotIter += iStIter;
 
 #ifdef DEBUG
                         if (DEBUG_LEVEL(MYDEBUG_FSTEPS)) {
                                 Out << "Step " << lStep
                                         << " time " << dTime + dCurrTimeStep
                                         << " step " << dCurrTimeStep
                                         << " iterations " << iStIter
                                         << " error " << dTest << std::endl;
                         }
 #endif /* DEBUG */
 
                         DEBUGLCOUT(MYDEBUG_FSTEPS, "Substep " << iSubStep
                                 << " of step " << lStep
                                 << " has been successfully completed "
                                 "in " << iStIter << " iterations"
                                 << std::endl);
 
                         if (mbdyn_stop_at_end_of_time_step()) {
                                 /* */
 #ifdef DEBUG_FICTITIOUS
                                 pDM->Output(0, dTime, dCurrTimeStep);
 #endif /* DEBUG_FICTITIOUS */
                                 Out << "Interrupted during dummy steps."
                                         << std::endl;
                                 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
                         }
 
                         dTime += dRefTimeStep;
                 }
                 if (outputMsg()) {
                         Out << "# Initial solution after dummy steps "
                                 "at time " << dTime
                                 << " performed in " << iStIter
                                 << " iterations with " << dTest
                                 << " error" << std::endl;
                 }
 
                 DEBUGLCOUT(MYDEBUG_FSTEPS,
                         "Dummy steps have been successfully completed "
                         "in " << iStIter << " iterations" << std::endl);
 #ifdef USE_EXTERNAL
         /* comunica che gli ultimi dati inviati sono la condizione iniziale */
                 External::SendInitial();
 #endif /* USE_EXTERNAL */
         } /* Fine dei passi fittizi */
 
         /* Output delle "condizioni iniziali" */
         bOut = pDM->Output(0, dTime, dCurrTimeStep);
 
         if (outputMsg()) {
                 Out
                         << "# Key for lines starting with \"Step\":"
                                 << std::endl
                         << "# Step Time TStep NIter ResErr SolErr SolConv Out"
                                 << std::endl
                         << "Step " << 0
                         << " " << dTime + dCurrTimeStep
                         << " " << dCurrTimeStep
                         << " " << iStIter
                         << " " << dTest
                         << " " << dSolTest
                         << " " << bSolConv
                         << " " << bOut
                         << std::endl;
         }
 
 
         if (eAbortAfter == AFTER_DUMMY_STEPS) {
                 Out << "End of dummy steps; no simulation is required."
                         << std::endl;
                 return false;
 
         } else if (mbdyn_stop_at_end_of_time_step()) {
                 /* Fa l'output della soluzione ed esce */
                 Out << "Interrupted during dummy steps." << std::endl;
                 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
         }
 
         eStatus = SOLVER_STATUS_PREPARED;
 
         return true;
 }

void Solver::PrintResidual	(	const VectorHandler &	Res,
		integer	iIterCnt
	)		const

virtual

Definition at line 5194 of file solver.cc.

References DataManager::PrintResidual().

Referenced by InverseDynamicsStepSolver::Advance(), LineSearchSolver::Residual(), BiCGStab::Solve(), Gmres::Solve(), NewtonRaphsonSolver::Solve(), and LineSearchSolver::Solve().

 {
         pDM->PrintResidual(Res, iIterCnt);
 }

Here is the call graph for this function:

void Solver::PrintSolution	(	const VectorHandler &	Sol,
		integer	iIterCnt
	)		const

virtual

Definition at line 5200 of file solver.cc.

References DataManager::PrintSolution().

Referenced by InverseDynamicsStepSolver::Advance(), BiCGStab::Solve(), NewtonRaphsonSolver::Solve(), Gmres::Solve(), and LineSearchSolver::Solve().

 {
         pDM->PrintSolution(Sol, iIterCnt);
 }

Here is the call graph for this function:

void Solver::ReadData ( MBDynParser & HP )

protected

Definition at line 1853 of file solver.cc.

Referenced by Prepare().

 {
         DEBUGCOUTFNAME("MultiStepIntegrator::ReadData");
 
         /* parole chiave */
         static const char*const sKeyWords[] = {
                 "begin",
                 "initial" "value",
                 "multistep",            /* deprecated */
                 "end",
 
                 "initial" "time",
                 "final" "time",
                 "time" "step",
                 "min" "time" "step",
                 "max" "time" "step",
                 "tolerance",
                 "max" "residual",
                 "max" "iterations",
                 "modify" "residual" "test",
                 "enforce" "constraint" "equations",
                 /* DEPRECATED */
                 "fictitious" "steps" "number",
                 "fictitious" "steps" "ratio",
                 "fictitious" "steps" "tolerance",
                 "fictitious" "steps" "max" "iterations",
                 /* END OF DEPRECATED */
 
                 "dummy" "steps" "number",
                 "dummy" "steps" "ratio",
                 "dummy" "steps" "tolerance",
                 "dummy" "steps" "max" "iterations",
 
                 "abort" "after",
                         "input",
                         "assembly",
                         "derivatives",
 
                         /* DEPRECATED */ "fictitious" "steps" /* END OF DEPRECATED */ ,
                         "dummy" "steps",
 
                 "output",
                         "none",
                         "iterations",
                         "residual",
                         "solution",
                         /* DEPRECATED */ "jacobian" /* END OF DEPRECATED */ ,
                         "jacobian" "matrix",
                         "bailout",
                         "messages",
                         "counter",
                         "matrix" "condition" "number",
                         "solver" "condition" "number",
                         "cpu" "time",
                 "output" "meter",
 
                 "method",
                 /* DEPRECATED */ "fictitious" "steps" "method" /* END OF DEPRECATED */ ,
                 "dummy" "steps" "method",
 
                 "Crank" "Nicolson",
                 /* DEPRECATED */ "Crank" "Nicholson" /* END OF DEPRECATED */ ,
                         /* DEPRECATED */ "nostro" /* END OF DEPRECATED */ ,
                         "ms",
                         "hope",
                         "bdf",
                         "thirdorder",
                         "implicit" "euler",
 
                 "derivatives" "coefficient",
                 "derivatives" "tolerance",
                 "derivatives" "max" "iterations",
 
                 /* DEPRECATED */
                 "true",
                 "modified",
                 /* END OF DEPRECATED */
 
                 "strategy",
                         "factor",
                         "no" "change",
                         "change",
 
                 "pod",
                 "eigen" "analysis",
 
                 /* DEPRECATED */
                 "solver",
                 "interface" "solver",
                 /* END OF DEPRECATED */
                 "linear" "solver",
                 "interface" "linear" "solver",
 
                 /* DEPRECATED */
                 "preconditioner",
                 /* END OF DEPRECATED */
 
                 "nonlinear" "solver",
                         "default",
                         "newton" "raphson",
                         "line" "search",
                         "matrix" "free",
                                 "bicgstab",
                                 "gmres",
                                         /* DEPRECATED */ "full" "jacobian" /* END OF DEPRECATED */ ,
                                         "full" "jacobian" "matrix",
 
                 /* RTAI stuff */
                 "real" "time",
 
                 /* multithread stuff */
                 "threads",
 
                 NULL
         };
 
         /* enum delle parole chiave */
         enum KeyWords {
                 UNKNOWN = -1,
                 BEGIN = 0,
                 INITIAL_VALUE,
                 MULTISTEP,
                 END,
 
                 INITIALTIME,
                 FINALTIME,
                 TIMESTEP,
                 MINTIMESTEP,
                 MAXTIMESTEP,
                 TOLERANCE,
                 MAXRESIDUAL,
                 MAXITERATIONS,
                 MODIFY_RES_TEST,
                 ENFORCE_CONSTRAINT_EQUATIONS,
                 FICTITIOUSSTEPSNUMBER,
                 FICTITIOUSSTEPSRATIO,
                 FICTITIOUSSTEPSTOLERANCE,
                 FICTITIOUSSTEPSMAXITERATIONS,
 
                 DUMMYSTEPSNUMBER,
                 DUMMYSTEPSRATIO,
                 DUMMYSTEPSTOLERANCE,
                 DUMMYSTEPSMAXITERATIONS,
 
                 ABORTAFTER,
                 INPUT,
                 ASSEMBLY,
                 DERIVATIVES,
                 FICTITIOUSSTEPS,
                 DUMMYSTEPS,
 
                 OUTPUT,
                         NONE,
                         ITERATIONS,
                         RESIDUAL,
                         SOLUTION,
                         JACOBIAN,
                         JACOBIANMATRIX,
                         BAILOUT,
                         MESSAGES,
                         COUNTER,
                         MATRIX_COND_NUM,
                         SOLVER_COND_NUM,
                         CPU_TIME,
                 OUTPUTMETER,
 
                 METHOD,
                 FICTITIOUSSTEPSMETHOD,
                 DUMMYSTEPSMETHOD,
                 CRANKNICOLSON,
                 CRANKNICHOLSON,
                 NOSTRO,
                 MS,
                 HOPE,
                 BDF,
                 THIRDORDER,
                 IMPLICITEULER,
 
                 DERIVATIVESCOEFFICIENT,
                 DERIVATIVESTOLERANCE,
                 DERIVATIVESMAXITERATIONS,
 
                 /* DEPRECATED */
                 NR_TRUE,
                 MODIFIED,
                 /* END OF DEPRECATED */
 
                 STRATEGY,
                 STRATEGYFACTOR,
                 STRATEGYNOCHANGE,
                 STRATEGYCHANGE,
 
                 POD,
                 EIGENANALYSIS,
 
                 /* DEPRECATED */
                 SOLVER,
                 INTERFACESOLVER,
                 /* END OF DEPRECATED */
                 LINEARSOLVER,
                 INTERFACELINEARSOLVER,
 
                 /* DEPRECATED */
                 PRECONDITIONER,
                 /* END OF DEPRECATED */
 
                 NONLINEARSOLVER,
                         DEFAULT,
                         NEWTONRAPHSON,
                         LINESEARCH,
                         MATRIXFREE,
                                 BICGSTAB,
                                 GMRES,
                                         FULLJACOBIAN,
                                         FULLJACOBIANMATRIX,
 
                 /* RTAI stuff */
                 REALTIME,
 
                 THREADS,
 
                 LASTKEYWORD
         };
 
         /* tabella delle parole chiave */
         KeyTable K(HP, sKeyWords);
 
         /* legge i dati della simulazione */
         if (KeyWords(HP.GetDescription()) != BEGIN) {
                 silent_cerr("Error: <begin> expected at line "
                         << HP.GetLineData() << "; aborting..." << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         switch (KeyWords(HP.GetWord())) {
         case MULTISTEP:
                 pedantic_cout("warning: \"begin: multistep\" is deprecated; "
                         "use \"begin: initial value;\" instead." << std::endl);
         case INITIAL_VALUE:
                 break;
 
         default:
                 silent_cerr("Error: \"begin: initial value;\" expected at line "
                         << HP.GetLineData() << "; aborting..." << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         bool bMethod(false);
         bool bDummyStepsMethod(false);
 
         /* dati letti qui ma da passare alle classi
          *      StepIntegration e NonlinearSolver
          */
 
         doublereal dTol = ::dDefaultTol;
         doublereal dSolutionTol = 0.;
         iMaxIterations = ::iDefaultMaxIterations;
         bool bModResTest = false;
         bool bSetScaleAlgebraic = false;
 
         /* Dati dei passi fittizi di trimmaggio iniziale */
         doublereal dDummyStepsTolerance = ::dDefaultDummyStepsTolerance;
         integer iDummyStepsMaxIterations = ::iDefaultMaxIterations;
 
         /* Dati del passo iniziale di calcolo delle derivate */
 
         doublereal dDerivativesTol = ::dDefaultTol;
         integer iDerivativesMaxIterations = ::iDefaultMaxIterations;
         integer iDerivativesCoefMaxIter = 0;
         doublereal dDerivativesCoefFactor = 10.;
 
 #ifdef USE_MULTITHREAD
         bool bSolverThreads(false);
         unsigned nSolverThreads = 0;
 #endif /* USE_MULTITHREAD */
 
 
         /* Ciclo infinito */
         while (true) {
                 KeyWords CurrKeyWord = KeyWords(HP.GetDescription());
 
                 switch (CurrKeyWord) {
                 case INITIALTIME:
                         dInitialTime = HP.GetReal();
                         DEBUGLCOUT(MYDEBUG_INPUT, "Initial time is "
                                 << dInitialTime << std::endl);
                         break;
 
                 case FINALTIME:
                         if (HP.IsKeyWord("forever")) {
                                 dFinalTime = std::numeric_limits<doublereal>::max();
                         } else {
                                 dFinalTime = HP.GetReal();
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT, "Final time is "
                                 << dFinalTime << std::endl);
 
                         if(dFinalTime <= dInitialTime) {
                                 silent_cerr("warning: final time " << dFinalTime
                                         << " is less than initial time "
                                         << dInitialTime << ';' << std::endl
                                         << "this will cause the simulation"
                                         " to abort" << std::endl);
                         }
                         break;
 
                 case TIMESTEP:
                         dInitialTimeStep = HP.GetReal();
                         DEBUGLCOUT(MYDEBUG_INPUT, "Initial time step is "
                                 << dInitialTimeStep << std::endl);
 
                         if (dInitialTimeStep == 0.) {
                                 silent_cerr("warning, null initial time step"
                                         " is not allowed" << std::endl);
                         } else if (dInitialTimeStep < 0.) {
                                 dInitialTimeStep = -dInitialTimeStep;
                                 silent_cerr("warning, negative initial time step"
                                         " is not allowed;" << std::endl
                                         << "its modulus " << dInitialTimeStep
                                         << " will be considered" << std::endl);
                         }
                         break;
 
                 case MINTIMESTEP:
                         try {
                                 dMinTimeStep = HP.GetReal(std::numeric_limits<doublereal>::min(), HighParser::range_gt<doublereal>(0.));
 
                         } catch (HighParser::ErrValueOutOfRange<doublereal> e) {
                                 silent_cerr("invalid min time step " << e.Get() << " (must be positive) [" << e.what() << "] at line " << HP.GetLineData() << std::endl);
                                 throw e;
                         }
                         break;
 
                 case MAXTIMESTEP:
                         if (HP.IsKeyWord("unlimited")) {
                                 DriveCaller *pDC = 0;
                                 SAFENEWWITHCONSTRUCTOR(pDC, ConstDriveCaller, ConstDriveCaller(0.));
                                 MaxTimeStep.Set(pDC);
 
                         } else {
                                 MaxTimeStep.Set(HP.GetDriveCaller());
                         }
 
 #ifdef DEBUG
                         if (typeid(*MaxTimeStep.pGetDriveCaller()) == typeid(PostponedDriveCaller)) {
                                 DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is postponed" << std::endl);
 
                         } else {
                                 DEBUGLCOUT(MYDEBUG_INPUT, "Max time step is " << MaxTimeStep.dGet() << std::endl);
                         }
 #endif // DEBUG
 
                         eTimeStepLimit = (HP.IsKeyWord("hard" "limit")
                                                                 && HP.GetYesNoOrBool())
                                                         ? TS_HARD_LIMIT
                                                         : TS_SOFT_LIMIT;
 
                         if (dGetInitialMaxTimeStep() == 0.) {
                                 silent_cout("no max time step limit will be"
                                         " considered" << std::endl);
 
                         } else if (dGetInitialMaxTimeStep() < 0.) {
                                 silent_cerr("negative max time step"
                                         " is not allowed" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         break;
 
                 case FICTITIOUSSTEPSNUMBER:
                 case DUMMYSTEPSNUMBER:
                         iDummyStepsNumber = HP.GetInt();
                         if (iDummyStepsNumber < 0) {
                                 iDummyStepsNumber =
                                         ::iDefaultDummyStepsNumber;
                                 silent_cerr("negative dummy steps number"
                                         " is illegal" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 
                         } else if (iDummyStepsNumber == 1) {
                                 silent_cerr("warning, a single dummy step"
                                         " may be useless" << std::endl);
                         }
 
                         DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps number: "
                                 << iDummyStepsNumber << std::endl);
                         break;
 
                 case FICTITIOUSSTEPSRATIO:
                 case DUMMYSTEPSRATIO:
                         dDummyStepsRatio = HP.GetReal();
                         if (dDummyStepsRatio < 0.) {
                                 silent_cerr("negative dummy steps ratio"
                                         " is illegal" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
 
                         if (dDummyStepsRatio > 1.) {
                                 silent_cerr("warning, dummy steps ratio"
                                         " is larger than one." << std::endl
                                         << "Something like 1.e-3 should"
                                         " be safer ..." << std::endl);
                         }
 
                         DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps ratio: "
                                 << dDummyStepsRatio << std::endl);
                         break;
 
                 case FICTITIOUSSTEPSTOLERANCE:
                 case DUMMYSTEPSTOLERANCE:
                         dDummyStepsTolerance = HP.GetReal();
                         if (dDummyStepsTolerance <= 0.) {
                                 silent_cerr("negative dummy steps"
                                         " tolerance is illegal" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT,
                                 "Dummy steps tolerance: "
                                 << dDummyStepsTolerance << std::endl);
                         break;
 
                 case ABORTAFTER: {
                         KeyWords WhenToAbort(KeyWords(HP.GetWord()));
                         switch (WhenToAbort) {
                         case INPUT:
                                 eAbortAfter = AFTER_INPUT;
                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Simulation will abort after"
                                         " data input" << std::endl);
                                 break;
 
                         case ASSEMBLY:
                                 eAbortAfter = AFTER_ASSEMBLY;
                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Simulation will abort after"
                                         " initial assembly" << std::endl);
                                 break;
 
                         case DERIVATIVES:
                                 eAbortAfter = AFTER_DERIVATIVES;
                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Simulation will abort after"
                                         " derivatives solution" << std::endl);
                                 break;
 
                         case FICTITIOUSSTEPS:
                         case DUMMYSTEPS:
                                 eAbortAfter = AFTER_DUMMY_STEPS;
                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Simulation will abort after"
                                         " dummy steps solution" << std::endl);
                                 break;
 
                         default:
                                 silent_cerr("Don't know when to abort,"
                                         " so I'm going to abort now" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                 } break;
 
                 case OUTPUT: {
                         unsigned OF = OUTPUT_DEFAULT;
                         bool setOutput = false;
 
                         while (HP.IsArg()) {
                                 KeyWords OutputFlag(KeyWords(HP.GetWord()));
                                 switch (OutputFlag) {
                                 case NONE:
                                         OF = OUTPUT_NONE;
                                         setOutput = true;
                                         break;
 
                                 case ITERATIONS:
                                         OF |= OUTPUT_ITERS;
                                         break;
 
                                 case RESIDUAL:
                                         OF |= OUTPUT_RES;
                                         break;
 
                                 case SOLUTION:
                                         OF |= OUTPUT_SOL;
                                         break;
 
                                 case JACOBIAN:
                                 case JACOBIANMATRIX:
                                         OF |= OUTPUT_JAC;
                                         break;
 
                                 case BAILOUT:
                                         OF |= OUTPUT_BAILOUT;
                                         break;
 
                                 case MESSAGES:
                                         OF |= OUTPUT_MSG;
                                         break;
 
                                 case COUNTER:
                                         OF |= OUTPUT_COUNTER;
                                         break;
 
                                 case MATRIX_COND_NUM:
                                         DelOutputFlags(OUTPUT_MAT_COND_NUM);
                                         if (HP.IsKeyWord("norm")) {
                                                 if (HP.IsKeyWord("inf")) {
                                                         OF |= OUTPUT_MAT_COND_NUM_INF;
                                                 } else {
                                                         const double P = HP.GetReal();
                                                         if (P != 1.) {
                                                                 silent_cerr("Only inf or 1 norm are supported for condition numbers at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                         OF |= OUTPUT_MAT_COND_NUM_1;
                                                 }
                                         } else {
                                                 OF |= OUTPUT_MAT_COND_NUM_1;
                                         }
                                         break;
 
                                 case SOLVER_COND_NUM:
                                         OF |= OUTPUT_SOLVER_COND_NUM;
                                         if (HP.IsKeyWord("stat")) {
                                                 if (HP.GetYesNoOrBool()) {
                                                         OF |= OUTPUT_SOLVER_COND_STAT;
                                                 } else {
                                                         DelOutputFlags(OUTPUT_SOLVER_COND_STAT);
                                                 }
                                         }
                                         break;
                                 case CPU_TIME:
                                         OF |= OUTPUT_CPU_TIME;
                                         break;
 
                                 default:
                                         silent_cerr("Warning, unknown output flag "
                                                 "at line " << HP.GetLineData()
                                                 << "; ignored" << std::endl);
                                         break;
                                 }
                         }
 
                         if (setOutput) {
                                 SetOutputFlags(OF);
 
                         } else {
                                 AddOutputFlags(OF);
                         }
                 } break;
 
                 case OUTPUTMETER:
                         SetOutputMeter(HP.GetDriveCaller(true));
                         break;
 
                 case METHOD: {
                         if (bMethod) {
                                 silent_cerr("error: multiple definition"
                                         " of integration method at line "
                                         << HP.GetLineData() << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         bMethod = true;
 
                         KeyWords KMethod = KeyWords(HP.GetWord());
                         switch (KMethod) {
                         case CRANKNICHOLSON:
                                 silent_cout("warning: \"crank nicolson\" is the correct spelling "
                                         "at line " << HP.GetLineData() << std::endl);
                         case CRANKNICOLSON:
                                 RegularType = INT_CRANKNICOLSON;
                                 break;
 
                         case BDF:
                                 /* default (order 2) */
                                 RegularType = INT_MS2;
 
                                 if (HP.IsKeyWord("order")) {
                                         int iOrder = HP.GetInt();
 
                                         switch (iOrder) {
                                         case 1:
                                                 RegularType = INT_IMPLICITEULER;
                                                 break;
 
                                         case 2:
                                                 break;
 
                                         default:
                                                 silent_cerr("unhandled BDF order " << iOrder << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                 }
 
                                 if (RegularType == INT_MS2) {
                                         SAFENEW(pRhoRegular, NullDriveCaller);
                                         SAFENEW(pRhoAlgebraicRegular, NullDriveCaller);
                                 }
                                 break;
 
                         case NOSTRO:
                                 silent_cerr("integration method \"nostro\" "
                                         "is deprecated; use \"ms\" "
                                         "instead at line "
                                         << HP.GetLineData()
                                         << std::endl);
                         case MS:
                         case HOPE:
                                 pRhoRegular = HP.GetDriveCaller(true);
 
                                 pRhoAlgebraicRegular = 0;
                                 if (HP.IsArg()) {
                                         pRhoAlgebraicRegular = HP.GetDriveCaller(true);
                                 } else {
                                         pRhoAlgebraicRegular = pRhoRegular->pCopy();
                                 }
 
                                 switch (KMethod) {
                                 case NOSTRO:
                                 case MS:
                                         RegularType = INT_MS2;
                                         break;
 
                                 case HOPE:
                                         RegularType = INT_HOPE;
                                         break;
 
                                 default:
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                                 break;
 
                         case THIRDORDER:
                                 if (HP.IsKeyWord("ad" "hoc")) {
                                         /* do nothing */ ;
                                 } else {
                                         pRhoRegular = HP.GetDriveCaller(true);
                                 }
                                 RegularType = INT_THIRDORDER;
                                 break;
 
                         case IMPLICITEULER:
                                 RegularType = INT_IMPLICITEULER;
                                 break;
 
                         default:
                                 silent_cerr("Unknown integration method at line "
                                         << HP.GetLineData() << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         break;
                 }
 
                 case FICTITIOUSSTEPSMETHOD:
                 case DUMMYSTEPSMETHOD: {
                         if (bDummyStepsMethod) {
                                 silent_cerr("error: multiple definition "
                                         "of dummy steps integration method "
                                         "at line " << HP.GetLineData()
                                         << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         bDummyStepsMethod = true;
 
                         KeyWords KMethod = KeyWords(HP.GetWord());
                         switch (KMethod) {
                         case CRANKNICHOLSON:
                                 silent_cout("warning: \"crank nicolson\" is the correct spelling" << std::endl);
                         case CRANKNICOLSON:
                                 DummyType = INT_CRANKNICOLSON;
                                 break;
 
                         case BDF:
                                 /* default (order 2) */
                                 DummyType = INT_MS2;
 
                                 if (HP.IsKeyWord("order")) {
                                         int iOrder = HP.GetInt();
 
                                         switch (iOrder) {
                                         case 1:
                                                 DummyType = INT_IMPLICITEULER;
                                                 break;
 
                                         case 2:
                                                 break;
 
                                         default:
                                                 silent_cerr("unhandled BDF order " << iOrder << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                 }
 
                                 if (DummyType == INT_MS2) {
                                         SAFENEW(pRhoDummy, NullDriveCaller);
                                         SAFENEW(pRhoAlgebraicDummy, NullDriveCaller);
                                 }
                                 break;
 
                         case NOSTRO:
                         case MS:
                         case HOPE:
                                 pRhoDummy = HP.GetDriveCaller(true);
 
                                 if (HP.IsArg()) {
                                         pRhoAlgebraicDummy = HP.GetDriveCaller(true);
                                 } else {
                                         pRhoAlgebraicDummy = pRhoDummy->pCopy();
                                 }
 
                                 switch (KMethod) {
                                 case NOSTRO:
                                 case MS:
                                         DummyType = INT_MS2;
                                         break;
 
                                 case HOPE:
                                         DummyType = INT_HOPE;
                                         break;
 
                                 default:
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                                 break;
 
                         case THIRDORDER:
                                 if (HP.IsKeyWord("ad" "hoc")) {
                                         /* do nothing */ ;
                                 } else {
                                         pRhoDummy = HP.GetDriveCaller(true);
                                 }
                                 DummyType = INT_THIRDORDER;
                                 break;
 
                         case IMPLICITEULER:
                                 DummyType = INT_IMPLICITEULER;
                                 break;
 
                         default:
                                 silent_cerr("Unknown integration method at line " << HP.GetLineData() << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         break;
                 }
 
                 case TOLERANCE: {
                         /*
                          * residual tolerance; can be the keyword "null",
                          * which means that the convergence test
                          * will be computed on the solution, or a number
                          */
                         if (HP.IsKeyWord("null")) {
                                 dTol = 0.;
 
                         } else {
                                 dTol = HP.GetReal();
                                 if (dTol < 0.) {
                                         dTol = ::dDefaultTol;
                                         silent_cerr("warning, residual tolerance "
                                                 "< 0. is illegal; "
                                                 "using default value " << dTol
                                                 << std::endl);
                                 }
                         }
 
                         /* safe default */
                         if (dTol == 0.) {
                                 ResTest = NonlinearSolverTest::NONE;
                         }
 
                         if (HP.IsArg()) {
                                 if (HP.IsKeyWord("test")) {
                                         if (HP.IsKeyWord("norm")) {
                                                 ResTest = NonlinearSolverTest::NORM;
                                         } else if (HP.IsKeyWord("minmax")) {
                                                 ResTest = NonlinearSolverTest::MINMAX;
                                         } else if (HP.IsKeyWord("none")) {
                                                 ResTest = NonlinearSolverTest::NONE;
                                         } else {
                                                 silent_cerr("unknown test "
                                                         "method at line "
                                                         << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         if (HP.IsKeyWord("scale")) {
                                                 if (ResTest == NonlinearSolverTest::NONE) {
                                                         silent_cerr("it's a nonsense "
                                                                 "to scale a disabled test; "
                                                                 "\"scale\" ignored"
                                                                 << std::endl);
                                                         bScale = false;
                                                 } else {
                                                         bScale = true;
                                                 }
                                         }
                                 }
                         }
 
                         if (HP.IsArg()) {
                                 if (!HP.IsKeyWord("null")) {
                                         dSolutionTol = HP.GetReal();
                                 }
 
                                 /* safe default */
                                 if (dSolutionTol != 0.) {
                                         SolTest = NonlinearSolverTest::NORM;
                                 }
 
                                 if (HP.IsArg()) {
                                         if (HP.IsKeyWord("test")) {
                                                 if (HP.IsKeyWord("norm")) {
                                                         SolTest = NonlinearSolverTest::NORM;
                                                 } else if (HP.IsKeyWord("minmax")) {
                                                         SolTest = NonlinearSolverTest::MINMAX;
                                                 } else if (HP.IsKeyWord("none")) {
                                                         SolTest = NonlinearSolverTest::NONE;
                                                 } else {
                                                         silent_cerr("unknown test "
                                                                 "method at line "
                                                                 << HP.GetLineData()
                                                                 << std::endl);
                                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                 }
                                         }
                                 }
 
                         } else if (dTol == 0.) {
                                 silent_cerr("need solution tolerance "
                                         "with null residual tolerance"
                                         << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
 
                         if (dSolutionTol < 0.) {
                                 dSolutionTol = 0.;
                                 silent_cerr("warning, solution tolerance "
                                         "< 0. is illegal; "
                                         "solution test is disabled"
                                         << std::endl);
                         }
 
                         if (dTol == 0. && dSolutionTol == 0.) {
                                 silent_cerr("both residual and solution "
                                         "tolerances are zero" << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
 
                         DEBUGLCOUT(MYDEBUG_INPUT, "tolerance = " << dTol
                                         << ", " << dSolutionTol << std::endl);
                         break;
                 }
 
                 case MAXRESIDUAL: {
                         if (HP.IsKeyWord("differential" "equations")) {
                                 dMaxResidualDiff = HP.GetReal();
 
                                 if (dMaxResidualDiff <= 0.) {
                                         silent_cerr("error: max residual for differential equations "
                                                                 "must be greater than zero at line "
                                                         << HP.GetLineData() << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                         }
 
                         if (HP.IsKeyWord("all" "equations") || HP.IsArg()) {
                                 dMaxResidual = HP.GetReal();
 
                                 if (dMaxResidual <= 0.) {
                                         silent_cerr("error: max residual must be greater than zero at line "
                                                         << HP.GetLineData() << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                         }
                         break;
                 }
 
                 case DERIVATIVESTOLERANCE: {
                         dDerivativesTol = HP.GetReal();
                         if (dDerivativesTol <= 0.) {
                                 dDerivativesTol = ::dDefaultTol;
                                 silent_cerr("warning, derivatives "
                                         "tolerance <= 0.0 is illegal; "
                                         "using default value "
                                         << dDerivativesTol
                                         << std::endl);
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Derivatives tolerance = "
                                         << dDerivativesTol
                                         << std::endl);
                         break;
                 }
 
                 case MAXITERATIONS: {
                         iMaxIterations = HP.GetInt();
                         if (iMaxIterations < 1) {
                                 iMaxIterations = ::iDefaultMaxIterations;
                                 silent_cerr("warning, max iterations "
                                         "< 1 is illegal; using default value "
                                         << iMaxIterations
                                         << std::endl);
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Max iterations = "
                                         << iMaxIterations << std::endl);
                         if (HP.IsKeyWord("at" "most")) {
                                 iMaxIterations = -iMaxIterations;
                         }
                         break;
                 }
 
                 case MODIFY_RES_TEST:
                         if (bParallel) {
                                 silent_cerr("\"modify residual test\" "
                                         "not supported by schur data manager "
                                         "at line " << HP.GetLineData()
                                         << "; ignored" << std::endl);
                         } else {
                                 bModResTest = true;
                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                         "Modify residual test" << std::endl);
                         }
                         break;
 
                 case ENFORCE_CONSTRAINT_EQUATIONS:
                         if (HP.IsKeyWord("constraint" "violations")) {
                                 eScaleFlags = SCALE_ALGEBRAIC_EQUATIONS_YES;
 
                                 bSetScaleAlgebraic = !HP.IsKeyWord("scale" "factor");
 
                                 if (!bSetScaleAlgebraic) {
                                         dScaleAlgebraic = HP.GetReal();
                                 }
                         } else if (HP.IsKeyWord("constraint" "violation" "rates")) {
                                 eScaleFlags = SCALE_ALGEBRAIC_EQUATIONS_NO;
                         } else {
                                 silent_cerr("Keyword \"constraint violations\" or "
                                                     "\"constraint violation rates\" expected at line "
                                                         << HP.GetLineData() << std::endl);
 
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         break;
 
                 case DERIVATIVESMAXITERATIONS: {
                         iDerivativesMaxIterations = HP.GetInt();
                         if (iDerivativesMaxIterations < 1) {
                                 iDerivativesMaxIterations = ::iDefaultMaxIterations;
                                 silent_cerr("warning, derivatives "
                                         "max iterations < 1 is illegal; "
                                         "using default value "
                                         << iDerivativesMaxIterations
                                         << std::endl);
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives "
                                         "max iterations = "
                                         << iDerivativesMaxIterations
                                         << std::endl);
                         break;
                 }
 
                 case FICTITIOUSSTEPSMAXITERATIONS:
                 case DUMMYSTEPSMAXITERATIONS: {
                         iDummyStepsMaxIterations = HP.GetInt();
                         if (iDummyStepsMaxIterations < 1) {
                                 iDummyStepsMaxIterations = ::iDefaultMaxIterations;
                                 silent_cerr("warning, dummy steps "
                                         "max iterations < 1 is illegal; "
                                         "using default value "
                                         << iDummyStepsMaxIterations
                                         << std::endl);
                         }
                         DEBUGLCOUT(MYDEBUG_INPUT, "Dummy steps "
                                         "max iterations = "
                                         << iDummyStepsMaxIterations
                                         << std::endl);
                         break;
                 }
 
                 case DERIVATIVESCOEFFICIENT: {
                         const bool bAutoDerivCoef = HP.IsKeyWord("auto");
 
                         if (!bAutoDerivCoef) {
                                 dDerivativesCoef = HP.GetReal();
                                 if (dDerivativesCoef <= 0.) {
                                         dDerivativesCoef = ::dDefaultDerivativesCoefficient;
                                         silent_cerr("warning, derivatives "
                                                 "coefficient <= 0. is illegal; "
                                                 "using default value "
                                                 << dDerivativesCoef
                                                 << std::endl);
                                 }
                                 DEBUGLCOUT(MYDEBUG_INPUT, "Derivatives coefficient = "
                                         << dDerivativesCoef << std::endl);
                         }
 
                         if (bAutoDerivCoef || HP.IsKeyWord("auto")) {
                                 if (HP.IsKeyWord("max" "iterations")) {
                                         iDerivativesCoefMaxIter = HP.GetInt();
                                         if (iDerivativesCoefMaxIter < 0) {
                                                 silent_cerr("number of iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else {
                                         iDerivativesCoefMaxIter = 6;
                                 }
 
                                 if (HP.IsKeyWord("factor")) {
                                         dDerivativesCoefFactor = HP.GetReal();
 
                                         if (dDerivativesCoefFactor <= 1.) {
                                                 silent_cerr("factor must be greater than one at line " << HP.GetLineData() << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                 }
                         }
 
                         break;
                 }
 
                 case NEWTONRAPHSON: {
                         pedantic_cout("Newton Raphson is deprecated; use "
                                         "\"nonlinear solver: newton raphson "
                                         "[ , modified, <steps> ]\" instead"
                                         << std::endl);
                         KeyWords NewRaph(KeyWords(HP.GetWord()));
                         switch(NewRaph) {
                         case MODIFIED:
                                 bTrueNewtonRaphson = 0;
                                 if (HP.IsArg()) {
                                         iIterationsBeforeAssembly = HP.GetInt();
                                 } else {
                                         iIterationsBeforeAssembly = ::iDefaultIterationsBeforeAssembly;
                                 }
                                 DEBUGLCOUT(MYDEBUG_INPUT, "Modified "
                                                 "Newton-Raphson will be used; "
                                                 "matrix will be assembled "
                                                 "at most after "
                                                 << iIterationsBeforeAssembly
                                                 << " iterations" << std::endl);
                                 break;
 
                         default:
                                 silent_cerr("warning: unknown case; "
                                         "using default" << std::endl);
 
                         /* no break: fall-thru to next case */
                         case NR_TRUE:
                                 bTrueNewtonRaphson = 1;
                                 iIterationsBeforeAssembly = 0;
                                 break;
                         }
                         break;
                 }
 
                 case END:
                         switch (KeyWords(HP.GetWord())) {
                         case MULTISTEP:
                                 pedantic_cout("\"end: multistep;\" is deprecated; "
                                         "use \"end: initial value;\" instead." << std::endl);
                         case INITIAL_VALUE:
                                 break;
 
                         default:
                                 silent_cerr("\"end: initial value;\" expected "
                                         "at line " << HP.GetLineData()
                                         << "; aborting..." << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         goto EndOfCycle;
 
                 case STRATEGY:
                         pTSC = ReadTimeStepData(this, HP);
                         break;
 
                 case POD:
                         silent_cerr("line " << HP.GetLineData()
                                 << ": POD analysis not supported (ignored)"
                                 << std::endl);
                         for (; HP.IsArg();) {
                                 (void)HP.GetReal();
                         }
                         break;
 
                 case EIGENANALYSIS:
                         // initialize output precision: 0 means use default precision
                         EigAn.iMatrixPrecision = 0;
                         EigAn.iResultsPrecision = 0;
 
                         // read eigenanalysis time (to be changed)
                         if (HP.IsKeyWord("list")) {
                                 int iNumTimes = HP.GetInt();
                                 if (iNumTimes <= 0) {
                                         silent_cerr("invalid number of eigenanalysis times "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
 
                                 EigAn.Analyses.resize(iNumTimes);
                                 for (std::vector<doublereal>::iterator i = EigAn.Analyses.begin();
                                         i != EigAn.Analyses.end(); ++i)
                                 {
                                         *i = HP.GetReal();
                                         if (i > EigAn.Analyses.begin() && *i <= *(i-1)) {
                                                 silent_cerr("eigenanalysis times must be in strict ascending order "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                 }
 
                         } else {
                                 EigAn.Analyses.resize(1);
                                 EigAn.Analyses[0] = HP.GetReal();
                         }
 
                         ASSERT(EigAn.Analyses.size() > 0);
                         // initialize EigAn
                         EigAn.currAnalysis = EigAn.Analyses.begin();
                         EigAn.bAnalysis = true;
 
                         // permute is the default; use "balance, no" to disable
                         EigAn.uFlags = EigenAnalysis::EIG_PERMUTE;
 
                         while (HP.IsArg()) {
                                 if (HP.IsKeyWord("parameter")) {
                                         EigAn.dParam = HP.GetReal();
 
                                 } else if (HP.IsKeyWord("output" "matrices")) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_MATRICES;
 
                                 } else if (HP.IsKeyWord("output" "full" "matrices")) {
 #ifndef USE_EIG
                                         silent_cerr("\"output full matrices\" needs eigenanalysis support" << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif // !USE_EIG
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
 
                                 } else if (HP.IsKeyWord("output" "sparse" "matrices")) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
 
                                 } else if (HP.IsKeyWord("output" "eigenvectors")) {
 #ifndef USE_EIG
                                         silent_cerr("\"output eigenvectors\" needs eigenanalysis support" << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif // !USE_EIG
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_EIGENVECTORS;
 
                                 } else if (HP.IsKeyWord("output" "geometry")) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_GEOMETRY;
 
                                 } else if (HP.IsKeyWord("matrix" "output" "precision")) {
                                         EigAn.iMatrixPrecision = HP.GetInt();
                                         if (EigAn.iMatrixPrecision <= 0) {
                                                 silent_cerr("negative or null \"matrix output precision\" "
                                                         "parameter at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else if (HP.IsKeyWord("results" "output" "precision")) {
                                         EigAn.iResultsPrecision = HP.GetInt();
                                         if (EigAn.iResultsPrecision <= 0) {
                                                 silent_cerr("negative or null \"results output precision\" "
                                                         "parameter at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else if (HP.IsKeyWord("upper" "frequency" "limit")) {
                                         EigAn.dUpperFreq = HP.GetReal();
                                         if (EigAn.dUpperFreq < 0.) {
                                                 silent_cerr("invalid \"upper frequency limit\" "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else if (HP.IsKeyWord("lower" "frequency" "limit")) {
                                         EigAn.dLowerFreq = HP.GetReal();
                                         if (EigAn.dLowerFreq < 0.) {
                                                 silent_cerr("invalid \"lower frequency limit\" "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else if (HP.IsKeyWord("use" "lapack")) {
                                         if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
                                                 silent_cerr("eigenanalysis routine already selected "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 #ifdef USE_LAPACK
                                         EigAn.uFlags |= EigenAnalysis::EIG_USE_LAPACK;
 #else // !USE_LAPACK
                                         silent_cerr("\"use lapack\" "
                                                 "needs to configure --with-lapack "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif // !USE_LAPACK
 
                                 } else if (HP.IsKeyWord("use" "arpack")) {
                                         if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
                                                 silent_cerr("eigenanalysis routine already selected "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 #ifdef USE_ARPACK
                                         EigAn.uFlags |= EigenAnalysis::EIG_USE_ARPACK;
 
                                         EigAn.arpack.iNEV = HP.GetInt();
                                         if (EigAn.arpack.iNEV <= 0) {
                                                 silent_cerr("invalid number of eigenvalues "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         EigAn.arpack.iNCV = HP.GetInt();
                                         if (EigAn.arpack.iNCV <= 0
                                                 || EigAn.arpack.iNCV <= EigAn.arpack.iNEV + 2)
                                         {
                                                 silent_cerr("invalid number of Arnoldi vectors "
                                                         "(must be > NEV+2) "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         if (EigAn.arpack.iNCV <= 2*EigAn.arpack.iNEV) {
                                                 silent_cerr("warning, possibly incorrect number of Arnoldi vectors "
                                                         "(should be > 2*NEV) "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                         }
 
                                         EigAn.arpack.dTOL = HP.GetReal();
                                         if (EigAn.arpack.dTOL < 0.) {
                                                 silent_cerr("tolerance must be non-negative "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 EigAn.arpack.dTOL = 0.;
                                         }
 #else // !USE_ARPACK
                                         silent_cerr("\"use arpack\" "
                                                 "needs to configure --with-arpack "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif // !USE_ARPACK
 
                                 } else if (HP.IsKeyWord("use" "jdqz")) {
                                         if (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
                                                 silent_cerr("eigenanalysis routine already selected "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 #ifdef USE_JDQZ
                                         EigAn.uFlags |= EigenAnalysis::EIG_USE_JDQZ;
 
                                         EigAn.jdqz.kmax = HP.GetInt();
                                         if (EigAn.jdqz.kmax <= 0) {
                                                 silent_cerr("invalid number of eigenvalues "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         EigAn.jdqz.jmax = HP.GetInt();
                                         if (EigAn.jdqz.jmax < 20
                                                 || EigAn.jdqz.jmax < 2*EigAn.jdqz.kmax)
                                         {
                                                 silent_cerr("invalid size of the search space "
                                                         "(must be >= 20 && >= 2*kmax) "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                         EigAn.jdqz.jmin = 2*EigAn.jdqz.kmax;
 
                                         EigAn.jdqz.eps = HP.GetReal();
                                         if (EigAn.jdqz.eps <= 0.) {
                                                 silent_cerr("tolerance must be non-negative "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 EigAn.jdqz.eps = std::numeric_limits<doublereal>::epsilon();
                                         }
 #else // !USE_JDQZ
                                         silent_cerr("\"use jdqz\" "
                                                 "needs to configure --with-jdqz "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
 #endif // !USE_JDQZ
 
                                 } else if (HP.IsKeyWord("balance")) {
                                         if (HP.IsKeyWord("no")) {
                                                 EigAn.uFlags &= ~EigenAnalysis::EIG_BALANCE;
 
                                         } else if (HP.IsKeyWord("permute")) {
                                                 EigAn.uFlags |= EigenAnalysis::EIG_PERMUTE;
 
                                         } else if (HP.IsKeyWord("scale")) {
                                                 EigAn.uFlags |= EigenAnalysis::EIG_SCALE;
 
                                         } else if (HP.IsKeyWord("all")) {
                                                 EigAn.uFlags |= EigenAnalysis::EIG_BALANCE;
 
                                         } else {
                                                 silent_cerr("unknown balance option "
                                                         "at line " << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
 
                                 } else if (HP.IsKeyWord("suffix" "width")) {
                                         if (HP.IsKeyWord("compute")) {
                                                 EigAn.iFNameWidth = EigenAnalysis::EIGAN_WIDTH_COMPUTE;
 
                                         } else {
                                                 EigAn.iFNameWidth = HP.GetInt(0, HighParser::range_ge<integer>(0));
                                         }
 
                                 } else if (HP.IsKeyWord("suffix" "format")) {
                                         EigAn.iFNameFormat = HP.GetStringWithDelims();
                                         // check?
 
                                 } else {
                                         silent_cerr("unknown option "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                         }
 
                         // lower must be less than upper
                         if (EigAn.dLowerFreq > EigAn.dUpperFreq) {
                                 silent_cerr("upper frequency limit " << EigAn.dUpperFreq
                                         << " less than lower frequency limit " << EigAn.dLowerFreq
                                         << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
 
                         // if only upper is defined, make lower equal to -upper
                         if (EigAn.dLowerFreq == -1.) {
                                 EigAn.dLowerFreq = -EigAn.dUpperFreq;
                         }
 
                         switch (EigAn.uFlags & EigenAnalysis::EIG_USE_MASK) {
                         case EigenAnalysis::EIG_USE_LAPACK:
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES) {
                                         silent_cerr("sparse matrices output "
                                                 "incompatible with lapack "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
                                 }
                                 break;
 
                         case EigenAnalysis::EIG_USE_ARPACK:
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
                                         silent_cerr("full matrices output "
                                                 "incompatible with arpack "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
                                 }
                                 break;
 
                         case EigenAnalysis::EIG_USE_JDQZ:
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_FULL_MATRICES) {
                                         silent_cerr("full matrices output "
                                                 "incompatible with jdqz "
                                                 "at line " << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
                                 }
                                 break;
 
                         default:
                                 break;
                         }
 
 #ifdef USE_EIG
                         // if an eigenanalysis routine is selected
                         // or sparse matrix output is not requested,
                         // force direct eigensolution
                         if ((EigAn.uFlags & EigenAnalysis::EIG_USE_MASK)
                                 || !(EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES))
                         {
                                 EigAn.uFlags |= EigenAnalysis::EIG_SOLVE;
                         }
 
                         // if no eigenanalysis routine is selected,
                         // force the use of LAPACK's
                         if ((EigAn.uFlags & EigenAnalysis::EIG_SOLVE)
                                 && !(EigAn.uFlags & EigenAnalysis::EIG_USE_MASK))
                         {
                                 EigAn.uFlags |= EigenAnalysis::EIG_USE_LAPACK;
                                 if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
                                         EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_FULL_MATRICES;
                                 }
                         }
 #else // !USE_EIG
                         if (EigAn.uFlags & EigenAnalysis::EIG_OUTPUT_MATRICES) {
                                 EigAn.uFlags |= EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES;
                         }
                         silent_cerr("warning: \"eigenanalysis\" not supported; ignored" << std::endl);
 #endif // !USE_EIG
                         break;
 
                 case SOLVER:
                         silent_cerr("\"solver\" keyword at line "
                                         << HP.GetLineData()
                                         << " is deprecated; "
                                         "use \"linear solver\" instead"
                                         << std::endl);
                 case LINEARSOLVER:
                         ReadLinSol(CurrLinearSolver, HP);
                         break;
 
                 case INTERFACESOLVER:
                         silent_cerr("\"interface solver\" keyword at line "
                                         << HP.GetLineData()
                                         << " is deprecated; "
                                         "use \"interface linear solver\" "
                                         "instead" << std::endl);
                 case INTERFACELINEARSOLVER:
                         ReadLinSol(CurrIntSolver, HP, true);
 
 #ifndef USE_MPI
                         silent_cerr("Interface solver only allowed "
                                 "when compiled with MPI support" << std::endl);
 #endif /* ! USE_MPI */
                         break;
 
                 case NONLINEARSOLVER:
                         switch (KeyWords(HP.GetWord())) {
                         case DEFAULT:
                                 NonlinearSolverType = NonlinearSolver::DEFAULT;
                                 break;
 
                         case NEWTONRAPHSON:
                                 NonlinearSolverType = NonlinearSolver::NEWTONRAPHSON;
                                 break;
 
                         case LINESEARCH:
                                 NonlinearSolverType = NonlinearSolver::LINESEARCH;
                                 break;
 
                         case MATRIXFREE:
                                 NonlinearSolverType = NonlinearSolver::MATRIXFREE;
                                 break;
 
                         default:
                                 silent_cerr("unknown nonlinear solver "
                                         "at line " << HP.GetLineData()
                                         << std::endl);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
 
                         switch (NonlinearSolverType) {
                         case NonlinearSolver::NEWTONRAPHSON:
                         case NonlinearSolver::LINESEARCH:
                                 bTrueNewtonRaphson = true;
                                 bKeepJac = false;
                                 iIterationsBeforeAssembly = 0;
 
                                 if (NonlinearSolverType == NonlinearSolver::NEWTONRAPHSON && HP.IsKeyWord("true")) {
                                         break;
                                 }
 
                                 if (HP.IsKeyWord("modified")) {
                                         bTrueNewtonRaphson = false;
                                         iIterationsBeforeAssembly = HP.GetInt();
 
                                         if (HP.IsKeyWord("keep" "jacobian")) {
                                                 pedantic_cout("Use of deprecated \"keep jacobian\" "
                                                         "at line " << HP.GetLineData() << std::endl);
                                                 bKeepJac = true;
 
                                         } else if (HP.IsKeyWord("keep" "jacobian" "matrix")) {
                                                 bKeepJac = true;
                                         }
 
                                         DEBUGLCOUT(MYDEBUG_INPUT, "modified "
                                                         "Newton-Raphson "
                                                         "will be used; "
                                                         "matrix will be "
                                                         "assembled at most "
                                                         "after "
                                                         << iIterationsBeforeAssembly
                                                         << " iterations"
                                                         << std::endl);
                                         if (HP.IsKeyWord("honor" "element" "requests")) {
                                                 bHonorJacRequest = true;
                                                 DEBUGLCOUT(MYDEBUG_INPUT,
                                                                 "honor elements' "
                                                                 "request to update "
                                                                 "the preconditioner"
                                                                 << std::endl);
                                         }
                                 }
 
                                 if (NonlinearSolver::LINESEARCH == NonlinearSolverType) {
                                         while (HP.IsArg()) {
                                                 if (HP.IsKeyWord("tolerance" "x")) {
                                                         LineSearch.dTolX = HP.GetReal();
                                                         if (LineSearch.dTolX < 0.) {
                                                                 silent_cerr("tolerance x must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("tolerance" "min")) {
                                                         LineSearch.dTolMin = HP.GetReal();
                                                         if (LineSearch.dTolMin < 0.) {
                                                                 silent_cerr("tolerance min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("max" "iterations")) {
                                                         LineSearch.iMaxIterations = HP.GetInt();
                                                         if (LineSearch.iMaxIterations < 0) {
                                                                 silent_cerr("max iterations must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("alpha")) {
                                                         LineSearch.dAlpha = HP.GetReal();
                                                         if (LineSearch.dAlpha < 0.) {
                                                                 silent_cerr("alpha must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("lambda" "min")) {
                                                         LineSearch.dLambdaMin = HP.GetReal();
                                                         if (LineSearch.dLambdaMin < 0.) {
                                                                 silent_cerr("lambda min must be greater than or equal to zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                         if (HP.IsKeyWord("relative")) {
                                                                 if (HP.GetYesNoOrBool()) {
                                                                         LineSearch.uFlags |= LineSearchParameters::RELATIVE_LAMBDA_MIN;
                                                                 } else {
                                                                         LineSearch.uFlags &= ~LineSearchParameters::RELATIVE_LAMBDA_MIN;
                                                                 }
                                                         }
                                                 } else if (HP.IsKeyWord("lambda" "factor" "min")) {
                                                         LineSearch.dLambdaFactMin = HP.GetReal();
                                                         if (LineSearch.dLambdaFactMin <= 0. || LineSearch.dLambdaFactMin >= 1.) {
                                                                 silent_cerr("lambda factor min must be in between zero and one at line" << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("max" "step")) {
                                                         LineSearch.dMaxStep = HP.GetReal();
                                                         if (LineSearch.dMaxStep <= 0.) {
                                                                 silent_cerr("max step must be greater than zero at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("zero" "gradient")) {
                                                         if (HP.IsKeyWord("continue")) {
                                                                 if (HP.GetYesNoOrBool()) {
                                                                         LineSearch.uFlags |= LineSearchParameters::ZERO_GRADIENT_CONTINUE;
                                                                 } else {
                                                                         LineSearch.uFlags &= ~LineSearchParameters::ZERO_GRADIENT_CONTINUE;
                                                                 }
                                                         }
                                                         else {
                                                                 silent_cerr("Keyword \"continue\" expected at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("divergence" "check")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::DIVERGENCE_CHECK;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::DIVERGENCE_CHECK;
                                                         }
                                                         if (HP.IsKeyWord("factor")) {
                                                                 LineSearch.dDivergenceCheck = HP.GetReal();
                                                                 if (LineSearch.dDivergenceCheck <= 0.) {
                                                                         silent_cerr("divergence check factor must be greater than zero at line " << HP.GetLineData() << std::endl);
                                                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                                 }
                                                         }
                                                 } else if (HP.IsKeyWord("algorithm")) {
                                                         LineSearch.uFlags &= ~LineSearchParameters::ALGORITHM;
                                                         if (HP.IsKeyWord("cubic")) {
                                                                 LineSearch.uFlags |= LineSearchParameters::ALGORITHM_CUBIC;
                                                         } else if (HP.IsKeyWord("factor")) {
                                                                 LineSearch.uFlags |= LineSearchParameters::ALGORITHM_FACTOR;
                                                         } else {
                                                                 silent_cerr("Keyword \"cubic\" or \"factor\" expected at line " << HP.GetLineData() << std::endl);
                                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                         }
                                                 } else if (HP.IsKeyWord("scale" "newton" "step")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::SCALE_NEWTON_STEP;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::SCALE_NEWTON_STEP;
                                                         }
                                                         if (HP.IsKeyWord("min" "scale")) {
                                                                 LineSearch.dMinStepScale = HP.GetReal();
                                                                 if (LineSearch.dMinStepScale < 0. || LineSearch.dMinStepScale > 1.) {
                                                                         silent_cerr("min scale must be in range [0-1] at line " << HP.GetLineData() << std::endl);
                                                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                                 }
                                                         }
                                                 } else if (HP.IsKeyWord("print" "convergence" "info")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::PRINT_CONVERGENCE_INFO;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::PRINT_CONVERGENCE_INFO;
                                                         }
                                                 } else if (HP.IsKeyWord("verbose")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::VERBOSE_MODE;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::VERBOSE_MODE;
                                                         }
                                                 } else if (HP.IsKeyWord("abort" "at" "lambda" "min")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::ABORT_AT_LAMBDA_MIN;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::ABORT_AT_LAMBDA_MIN;
                                                         }
                                                 } else if (HP.IsKeyWord("non" "negative" "slope" "continue")) {
                                                         if (HP.GetYesNoOrBool()) {
                                                                 LineSearch.uFlags |= LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE;
                                                         } else {
                                                                 LineSearch.uFlags &= ~LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE;
                                                         }
                                                 } else {
                                                         silent_cerr("Keyword \"tolerance x\", "
                                                                 "\"tolerance min\", \"max iterations\", \"alpha\", "
                                                                 "\"lambda min\" \"lambda factor min\", \"max step\" "
                                                                 "\"divergence check\", \"algorithm\", \"scale newton step\" "
                                                                 "\"print convergence info\", \"verbose\", "
                                                                 "\"abort at lambda min\" "
                                                                 "or \"zero gradient\" expected at line " << HP.GetLineData() << std::endl);
                                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                                 }
                                         }
                                 }
                                 break;
 
                         case NonlinearSolver::MATRIXFREE:
                                 switch (KeyWords(HP.GetWord())) {
                                 case DEFAULT:
                                         MFSolverType = MatrixFreeSolver::DEFAULT;
                                         break;
 
 
                                 case BICGSTAB:
                                         MFSolverType = MatrixFreeSolver::BICGSTAB;
                                         break;
 
                                 case GMRES:
                                         MFSolverType = MatrixFreeSolver::GMRES;
                                         break;
 
                                 default:
                                         silent_cerr("unknown iterative "
                                                 "solver at line "
                                                 << HP.GetLineData()
                                                 << std::endl);
                                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                 }
 
                                 if (HP.IsKeyWord("tolerance")) {
                                         dIterTol = HP.GetReal();
                                         DEBUGLCOUT(MYDEBUG_INPUT,"inner "
                                                         "iterative solver "
                                                         "tolerance: "
                                                         << dIterTol
                                                         << std::endl);
                                 }
 
                                 if (HP.IsKeyWord("steps")) {
                                         iIterativeMaxSteps = HP.GetInt();
                                         DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
                                                         "number of inner "
                                                         "steps for iterative "
                                                         "solver: "
                                                         << iIterativeMaxSteps
                                                         << std::endl);
                                 }
 
                                 if (HP.IsKeyWord("tau")) {
                                         dIterertiveTau = HP.GetReal();
                                         DEBUGLCOUT(MYDEBUG_INPUT,
                                                         "tau scaling "
                                                         "coefficient "
                                                         "for iterative "
                                                         "solver: "
                                                         << dIterertiveTau
                                                         << std::endl);
                                 }
 
                                 if (HP.IsKeyWord("eta")) {
                                         dIterertiveEtaMax = HP.GetReal();
                                         DEBUGLCOUT(MYDEBUG_INPUT, "maximum "
                                                         "eta coefficient "
                                                         "for iterative "
                                                         "solver: "
                                                         << dIterertiveEtaMax
                                                         << std::endl);
                                 }
 
                                 if (HP.IsKeyWord("preconditioner")) {
                                         KeyWords KPrecond = KeyWords(HP.GetWord());
                                         switch (KPrecond) {
                                         case FULLJACOBIAN:
                                         case FULLJACOBIANMATRIX:
                                                 PcType = Preconditioner::FULLJACOBIANMATRIX;
                                                 if (HP.IsKeyWord("steps")) {
                                                         iPrecondSteps = HP.GetInt();
                                                         DEBUGLCOUT(MYDEBUG_INPUT,
                                                                         "number of steps "
                                                                         "before recomputing "
                                                                         "the preconditioner: "
                                                                         << iPrecondSteps
                                                                         << std::endl);
                                                 }
                                                 if (HP.IsKeyWord("honor" "element" "requests")) {
                                                         bHonorJacRequest = true;
                                                         DEBUGLCOUT(MYDEBUG_INPUT,
                                                                         "honor elements' "
                                                                         "request to update "
                                                                         "the preconditioner"
                                                                         << std::endl);
                                                 }
                                                 break;
 
                                                 /* add other preconditioners
                                                  * here */
 
                                         default:
                                                 silent_cerr("unknown "
                                                         "preconditioner "
                                                         "at line "
                                                         << HP.GetLineData()
                                                         << std::endl);
                                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                                         }
                                         break;
                                 }
                                 break;
 
                         default:
                                 ASSERT(0);
                                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         }
                         break;
 
                 case REALTIME:
                         pRTSolver = ReadRTSolver(this, HP);
                         break;
 
                 case THREADS:
                         if (HP.IsKeyWord("auto")) {
 #ifdef USE_MULTITHREAD
                                 int n = get_nprocs();
                                 /* sanity checks ... */
                                 if (n <= 0) {
                                         silent_cerr("got " << n << " CPUs "
                                                         "at line "
                                                         << HP.GetLineData()
                                                         << std::endl);
                                         nThreads = 1;
                                 } else {
                                         nThreads = n;
                                 }
 #else /* ! USE_MULTITHREAD */
                                 silent_cerr("configure with "
                                                 "--enable-multithread "
                                                 "for multithreaded assembly"
                                                 << std::endl);
 #endif /* ! USE_MULTITHREAD */
 
                         } else if (HP.IsKeyWord("disable")) {
 #ifdef USE_MULTITHREAD
                                 nThreads = 1;
 #endif /* USE_MULTITHREAD */
 
                         } else {
 #ifdef USE_MULTITHREAD
                                 bool bAssembly = false;
                                 bool bSolver = false;
                                 bool bAll = true;
                                 unsigned nt;
 #endif // USE_MULTITHREAD
 
                                 if (HP.IsKeyWord("assembly")) {
 #ifdef USE_MULTITHREAD
                                         bAll = false;
                                         bAssembly = true;
 #endif // USE_MULTITHREAD
 
                                 } else if (HP.IsKeyWord("solver")) {
 #ifdef USE_MULTITHREAD
                                         bAll = false;
                                         bSolver = true;
 #endif // USE_MULTITHREAD
                                 }
 
 #ifdef USE_MULTITHREAD
                                 nt =
 #endif // USE_MULTITHREAD
                                 HP.GetInt();
 
 #ifdef USE_MULTITHREAD
                                 if (bAll || bAssembly) {
                                         nThreads = nt;
                                 }
 
                                 if (bAll || bSolver) {
                                         bSolverThreads = true;
                                         nSolverThreads = nt;
                                 }
 #else /* ! USE_MULTITHREAD */
                                 silent_cerr("configure with "
                                                 "--enable-multithread "
                                                 "for multithreaded assembly"
                                                 << std::endl);
 #endif /* ! USE_MULTITHREAD */
                         }
                         break;
 
 
                 default:
                         silent_cerr("unknown description at line "
                                 << HP.GetLineData() << "; aborting..."
                                 << std::endl);
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 }
         }
 
 EndOfCycle: /* esce dal ciclo di lettura */
 
         if (pTSC == 0) {
                 pedantic_cout("No time step control strategy defined; defaulting to NoChange time step control" );
                 pTSC = new NoChange(this);
         }
 
         if (dFinalTime < dInitialTime) {
                 eAbortAfter = AFTER_ASSEMBLY;
         }
 
         if (dFinalTime == dInitialTime) {
                 eAbortAfter = AFTER_DERIVATIVES;
         }
 
         /* Metodo di integrazione di default */
         if (!bMethod) {
                 ASSERT(RegularType == INT_UNKNOWN);
 
                 /* FIXME: maybe we should use a better value
                  * like 0.6; however, BDF should be conservative */
                 SAFENEW(pRhoRegular, NullDriveCaller);
 
                 /* DriveCaller per Rho asintotico per variabili algebriche */
                 pRhoAlgebraicRegular = pRhoRegular->pCopy();
 
                 RegularType = INT_MS2;
         }
 
         /* Metodo di integrazione di default */
         if (iDummyStepsNumber && !bDummyStepsMethod) {
                 ASSERT(DummyType == INT_UNKNOWN);
 
                 SAFENEW(pRhoDummy, NullDriveCaller);
 
                 /* DriveCaller per Rho asintotico per variabili algebriche */
                 pRhoAlgebraicDummy = pRhoDummy->pCopy();
 
                 DummyType = INT_MS2;
         }
 
         /* costruzione dello step solver derivative */
         SAFENEWWITHCONSTRUCTOR(pDerivativeSteps,
                         DerivativeSolver,
                         DerivativeSolver(dDerivativesTol,
                                 0.,
                                 dInitialTimeStep*dDerivativesCoef,
                                 iDerivativesMaxIterations,
                                 bModResTest,
                                 iDerivativesCoefMaxIter,
                                 dDerivativesCoefFactor));
 
         /* First step prediction must always be Crank-Nicolson for accuracy */
         if (iDummyStepsNumber) {
                 SAFENEWWITHCONSTRUCTOR(pFirstDummyStep,
                                 CrankNicolsonIntegrator,
                                 CrankNicolsonIntegrator(dDummyStepsTolerance,
                                         0.,
                                         iDummyStepsMaxIterations,
                                         bModResTest));
 
                 /* costruzione dello step solver dummy */
                 switch (DummyType) {
                 case INT_CRANKNICOLSON:
                         SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                         CrankNicolsonIntegrator,
                                         CrankNicolsonIntegrator(dDummyStepsTolerance,
                                                 0.,
                                                 iDummyStepsMaxIterations,
                                                 bModResTest));
                         break;
 
                 case INT_MS2:
                         SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                         MultistepSolver,
                                         MultistepSolver(dDummyStepsTolerance,
                                                 0.,
                                                 iDummyStepsMaxIterations,
                                                 pRhoDummy,
                                                 pRhoAlgebraicDummy,
                                                 bModResTest));
                         break;
 
                 case INT_HOPE:
                         SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                         HopeSolver,
                                         HopeSolver(dDummyStepsTolerance,
                                                 dSolutionTol,
                                                 iDummyStepsMaxIterations,
                                                 pRhoDummy,
                                                 pRhoAlgebraicDummy,
                                                 bModResTest));
                         break;
 
                 case INT_THIRDORDER:
                         if (pRhoDummy == 0) {
                                 SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                                 AdHocThirdOrderIntegrator,
                                                 AdHocThirdOrderIntegrator(dDummyStepsTolerance,
                                                         dSolutionTol,
                                                         iDummyStepsMaxIterations,
                                                         bModResTest));
                         } else {
                                 SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                                 TunableThirdOrderIntegrator,
                                                 TunableThirdOrderIntegrator(dDummyStepsTolerance,
                                                         dSolutionTol,
                                                         iDummyStepsMaxIterations,
                                                         pRhoDummy,
                                                         bModResTest));
                         }
                         break;
 
                 case INT_IMPLICITEULER:
                         SAFENEWWITHCONSTRUCTOR(pDummySteps,
                                         ImplicitEulerIntegrator,
                                         ImplicitEulerIntegrator(dDummyStepsTolerance,
                                                 dSolutionTol, iDummyStepsMaxIterations,
                                                 bModResTest));
                         break;
 
                 default:
                         silent_cerr("unknown dummy steps integration method"
                                 << std::endl);
                         throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                         break;
                 }
         }
 
         SAFENEWWITHCONSTRUCTOR(pFirstRegularStep,
                         CrankNicolsonIntegrator,
                         CrankNicolsonIntegrator(dTol,
                                 dSolutionTol,
                                 iMaxIterations,
                                 bModResTest));
 
         /* costruzione dello step solver per i passi normali */
         switch (RegularType) {
         case INT_CRANKNICOLSON:
                 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                         CrankNicolsonIntegrator,
                         CrankNicolsonIntegrator(dTol,
                                 dSolutionTol,
                                 iMaxIterations,
                                 bModResTest));
                 break;
 
         case INT_MS2:
                 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                                 MultistepSolver,
                                 MultistepSolver(dTol,
                                         dSolutionTol,
                                         iMaxIterations,
                                         pRhoRegular,
                                         pRhoAlgebraicRegular,
                                         bModResTest));
                 break;
 
         case INT_HOPE:
                 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                                 HopeSolver,
                                 HopeSolver(dTol,
                                         dSolutionTol,
                                         iMaxIterations,
                                         pRhoRegular,
                                         pRhoAlgebraicRegular,
                                         bModResTest));
                 break;
 
         case INT_THIRDORDER:
                 if (pRhoRegular == 0) {
                         SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                                         AdHocThirdOrderIntegrator,
                                         AdHocThirdOrderIntegrator(dTol,
                                                 dSolutionTol,
                                                 iMaxIterations,
                                                 bModResTest));
                 } else {
                         SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                                         TunableThirdOrderIntegrator,
                                         TunableThirdOrderIntegrator(dTol,
                                                 dSolutionTol,
                                                 iMaxIterations,
                                                 pRhoRegular,
                                                 bModResTest));
                 }
                 break;
 
         case INT_IMPLICITEULER:
                 SAFENEWWITHCONSTRUCTOR(pRegularSteps,
                                 ImplicitEulerIntegrator,
                                 ImplicitEulerIntegrator(dTol,
                                         dSolutionTol,
                                         iMaxIterations,
                                         bModResTest));
                 break;
 
         default:
                 silent_cerr("Unknown integration method" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
                 break;
         }
 
         if (bSetScaleAlgebraic) {
                 dScaleAlgebraic = 1. / dInitialTimeStep;
         }
 
 #ifdef USE_MULTITHREAD
         if (bSolverThreads) {
                 if (CurrLinearSolver.SetNumThreads(nSolverThreads)) {
                         silent_cerr("linear solver "
                                         << CurrLinearSolver.GetSolverName()
                                         << " does not support "
                                         "threaded solution" << std::endl);
                 }
         }
 #endif /* USE_MULTITHREAD */
 }

Here is the call graph for this function:

std::ostream & Solver::Restart	(	std::ostream &	out,
		DataManager::eRestart	type
	)		const

Definition at line 1717 of file solver.cc.

References ASSERT, DataManager::ATEND, NonlinearSolverOptions::bHonorJacRequest, bKeepJac, bScale, bTrueNewtonRaphson, CurrLinearSolver, dDerivativesCoef, dDummyStepsRatio, dFinalTime, DataManager::dGetTime(), dInitialTimeStep, StepIntegrator::GetIntegratorDSolTol(), StepIntegrator::GetIntegratorDTol(), StepIntegrator::GetIntegratorMaxIters(), iDummyStepsNumber, iIterationsBeforeAssembly, INT_CRANKNICOLSON, INT_HOPE, INT_IMPLICITEULER, INT_MS2, INT_THIRDORDER, DataManager::ITERATIONS, NonlinearSolver::MATRIXFREE, MBDYN_EXCEPT_ARGS, NonlinearSolverTest::MINMAX, NonlinearSolver::NEWTONRAPHSON, NO_OP, NonlinearSolverTest::NONE, NonlinearSolverType, NonlinearSolverTest::NORM, pDerivativeSteps, pDM, pDummySteps, pRegularSteps, pRhoAlgebraicRegular, pRhoRegular, RegularType, DriveCaller::Restart(), RestartLinSol(), ResTest, SolTest, DataManager::TIME, and DataManager::TIMES.

Referenced by DataManager::MakeRestart().

 {
 
         out << "begin: initial value;" << std::endl;
         switch(type) {
         case DataManager::ATEND:
                 out << "  #  initial time: " << pDM->dGetTime() << ";"
                         << std::endl
                         << "  #  final time: " << dFinalTime << ";"
                         << std::endl
                         << "  #  time step: " << dInitialTimeStep << ";"
                         << std::endl;
                 break;
         case DataManager::ITERATIONS:
         case DataManager::TIME:
         case DataManager::TIMES:
                 out << "  initial time: " << pDM->dGetTime()<< ";" << std::endl
                         << "  final time: " << dFinalTime << ";" << std::endl
                         << "  time step: " << dInitialTimeStep << ";"
                         << std::endl;
                 break;
         default:
                 ASSERT(0);
         }
 
         out << "  method: ";
         switch(RegularType) {
         case INT_CRANKNICOLSON:
                 out << "Crank Nicolson; " << std::endl;
                 break;
         case INT_MS2:
                 out << "ms, ";
                 pRhoRegular->Restart(out) << ", ";
                 pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
                 break;
         case INT_HOPE:
                 out << "hope, ";
                 pRhoRegular->Restart(out) << ", ";
                 pRhoAlgebraicRegular->Restart(out) << ";" << std::endl;
                 break;
 
         case INT_THIRDORDER:
                 out << "thirdorder, ";
                 if (!pRhoRegular)
                         out << "ad hoc;" << std::endl;
                         else
                         pRhoRegular->Restart(out) << ";" << std::endl;
                 break;
         case INT_IMPLICITEULER:
                 out << "implicit euler;" << std::endl;
                 break;
         default:
                 ASSERT(0);
         }
 
         integer iMI = pRegularSteps->GetIntegratorMaxIters();
         out << "  max iterations: " << std::abs(iMI);
         if (iMI < 0) out << ", at most";
         out
                 << ";" << std::endl
                 << "  tolerance: " << pRegularSteps->GetIntegratorDTol();
         switch(ResTest) {
         case NonlinearSolverTest::NORM:
                 out << ", test, norm" ;
                 break;
         case NonlinearSolverTest::MINMAX:
                 out << ", test, minmax" ;
                 break;
         case NonlinearSolverTest::NONE:
                 NO_OP;
         default:
                 silent_cerr("unhandled nonlinear solver test type" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         if (bScale) {
                 out << ", scale"
                         << ", " << pRegularSteps->GetIntegratorDSolTol();
         }
 
         switch (SolTest) {
         case NonlinearSolverTest::NORM:
                 out << ", test, norm" ;
                 break;
         case NonlinearSolverTest::MINMAX:
                 out << ", test, minmax" ;
                 break;
         case NonlinearSolverTest::NONE:
                 NO_OP;
         default:
                 ASSERT(0);
         }
         out
                 << ";" << std::endl;
 
         if ( pDerivativeSteps != 0 ) {
                 out
                         << "  derivatives max iterations: " << pDerivativeSteps->GetIntegratorMaxIters() << ";" << std::endl
                         << "  derivatives tolerance: " << pDerivativeSteps->GetIntegratorDTol() << ";" << std::endl
                         << "  derivatives coefficient: " << dDerivativesCoef << ";" << std::endl;
         }
 
         if (iDummyStepsNumber) {
                 out
                         << "  dummy steps max iterations: " << pDummySteps->GetIntegratorMaxIters() << ";" << std::endl
                         << "  dummy steps tolerance: " << pDummySteps->GetIntegratorDTol() << ";" << std::endl;
         }
         out
                 << "  dummy steps number: " << iDummyStepsNumber << ";" << std::endl
                 << "  dummy steps ratio: " << dDummyStepsRatio << ";" << std::endl;
         switch (NonlinearSolverType) {
         case NonlinearSolver::MATRIXFREE:
                 out << "  #  nonlinear solver: matrix free;" << std::endl;
                 break;
         case NonlinearSolver::NEWTONRAPHSON:
         default :
                 out << "  nonlinear solver: newton raphson";
                 if (!bTrueNewtonRaphson) {
                         out << ", modified, " << iIterationsBeforeAssembly;
                         if (bKeepJac) {
                                 out << ", keep jacobian matrix";
                         }
                         if (bHonorJacRequest) {
                                 out << ", honor element requests";
                         }
                 }
                 out << ";" << std::endl;
         }
         out << "  solver: ";
         RestartLinSol(out, CurrLinearSolver);
         out << "end: initial value;" << std::endl << std::endl;
         return out;
 }

Here is the call graph for this function:

void Solver::Run ( void )

Definition at line 1633 of file solver.cc.

References Advance(), DEBUGCOUTFNAME, NO_OP, Prepare(), and Start().

Referenced by RunMBDyn().

 {
         DEBUGCOUTFNAME("Solver::Run");
 
         if (Prepare()) {
                 if (Start()) {
                         while (Advance()) {
                                 NO_OP;
                         }
                 }
         }
 }

Here is the call graph for this function:

void Solver::SetupSolmans	(	integer	iStates,
		bool	bCanBeParallel = `false`
	)

protected

Definition at line 5135 of file solver.cc.

References AllocateSchurSolman(), AllocateSolman(), bParallel, CurrLinearSolver, DEBUGLCOUT, LinSol::iGetWorkSpaceSize(), iNumDofs, iNumLocDofs, iUnkStates, MBDYN_EXCEPT_ARGS, MYDEBUG_MEM, pLocalSM, and SAFEDELETE.

Referenced by InverseSolver::Prepare(), Prepare(), and Start().

 {
         DEBUGLCOUT(MYDEBUG_MEM, "creating SolutionManager\n\tsize = "
                 << iNumDofs*iUnkStates <<
                 "\n\tnumdofs = " << iNumDofs
                 << "\n\tnumstates = " << iStates << std::endl);
 
         /* delete previous solmans */
         if (pSM != 0) {
                 SAFEDELETE(pSM);
                 pSM = 0;
         }
         if (pLocalSM != 0) {
                 SAFEDELETE(pLocalSM);
                 pLocalSM = 0;
         }
 
         integer iWorkSpaceSize = CurrLinearSolver.iGetWorkSpaceSize();
         integer iLWS = iWorkSpaceSize;
         integer iNLD = iNumDofs*iStates;
         if (bCanBeParallel && bParallel) {
                 /* FIXME BEPPE! */
                 iLWS = iWorkSpaceSize*iNumLocDofs/(iNumDofs*iNumDofs);
                 /* FIXME: GIUSTO QUESTO? */
                 iNLD = iNumLocDofs*iStates;
         }
 
         SolutionManager *pCurrSM = AllocateSolman(iNLD, iLWS);
 
         /*
          * This is the LOCAL solver if instantiating a parallel
          * integrator; otherwise it is the MAIN solver
          */
         if (bCanBeParallel && bParallel) {
                 pLocalSM = pCurrSM;
 
                 /* Crea il solutore di Schur globale */
                 pSM = AllocateSchurSolman(iStates);
 
         } else {
                 pSM = pCurrSM;
         }
         /*
          * FIXME: at present there MUST be a pSM
          * (even for matrix-free nonlinear solvers)
          */
         if (pSM == 0) {
                 silent_cerr("No linear solver defined" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 }

Here is the call graph for this function:

bool Solver::Start ( void )

virtual

Reimplemented in InverseSolver.

Definition at line 1202 of file solver.cc.

References StepIntegrator::Advance(), Solver::EigenAnalysis::Analyses, ASSERT, Solver::EigenAnalysis::bAnalysis, DataManager::BeforePredict(), bOut, bSolConv, Solver::EigenAnalysis::currAnalysis, CurrStep, dCurrTimeStep, DEBUGCOUT, DEBUGCOUTFNAME, TimeStepControl::dGetNewStepTime(), dInitialTimeStep, dMinTimeStep, dRefTimeStep, dSolTest, dTest, dTime, dTotErr, Eig(), EigAn, eStatus, Flip(), DataManager::GetDofDescription(), StepIntegrator::GetIntegratorMaxIters(), StepIntegrator::GetIntegratorNumUnknownStates(), DataManager::GetOutFile(), LinearSolver::ErrFactor::iCol, iMaxIterations, TimeStepControl::Init(), RTSolverBase::Init(), iStIter, iTotIter, lStep, MaxTimeStep, MBDYN_EXCEPT_ARGS, mbdyn_stop_at_end_of_time_step(), StepIntegrator::NEWSTEP, DataManager::Output(), SolverDiagnostics::outputCounter(), SolverDiagnostics::outputMsg(), SolverDiagnostics::outputStep(), pCurrStepIntegrator, pDM, pFirstRegularStep, pNLS, pRegularSteps, pRTSolver, pTSC, pX, pXPrime, qX, qXPrime, StepIntegrator::REPEATSTEP, SAFEDELETE, DataManager::SetTime(), SetupSolmans(), SOLVER_STATUS_PREPARED, SOLVER_STATUS_STARTED, and MBDynErrBase::what().

Referenced by mb_sol_start(), and Run().

 {
         DEBUGCOUTFNAME("Solver::Start");
 
         // consistency check
         if (eStatus != SOLVER_STATUS_PREPARED) {
                 silent_cerr("Start() must be called after Prepare()" << std::endl);
                 throw ErrGeneric(MBDYN_EXCEPT_ARGS);
         }
 
         /* primo passo regolare */
 
 #ifdef USE_EXTERNAL
         /* il prossimo passo e' un regular */
         pNLS->SetExternal(External::REGULAR);
 #endif /* USE_EXTERNAL */
 
         lStep = 1; /* Resetto di nuovo lStep */
 
         DEBUGCOUT("Step " << lStep << " has been successfully completed "
                         "in " << iStIter << " iterations" << std::endl);
 
 
         DEBUGCOUT("Current time step: " << dCurrTimeStep << std::endl);
 
         pDM->BeforePredict(*pX, *pXPrime, *qX[0], *qXPrime[0]);
 
         Flip();
         dRefTimeStep = dInitialTimeStep;
         dCurrTimeStep = dRefTimeStep;
 
         CurrStep = StepIntegrator::NEWSTEP;
         pTSC->Init(iMaxIterations, dMinTimeStep, MaxTimeStep, dInitialTimeStep);
 
         /* Setup SolutionManager(s) */
         ASSERT(pFirstRegularStep!= 0);
         SetupSolmans(pFirstRegularStep->GetIntegratorNumUnknownStates(), true);
         pCurrStepIntegrator = pFirstRegularStep;
 
 IfFirstStepIsToBeRepeated:
         try {
                 pDM->SetTime(dTime + dCurrTimeStep, dCurrTimeStep, 1);
                 if (outputStep()) {
                         if (outputCounter()) {
                                 silent_cout(std::endl);
                         }
                         silent_cout("Step(" << 1 << ':' << 0 << ") t=" << dTime + dCurrTimeStep << " dt=" << dCurrTimeStep << std::endl);
                 }
                 dTest = pFirstRegularStep->Advance(this, dRefTimeStep,
                                 dCurrTimeStep/dRefTimeStep, CurrStep,
                                 qX, qXPrime, pX, pXPrime,
                                 iStIter, dTest, dSolTest);
         }
         catch (NonlinearSolver::NoConvergence) {
                 if (dCurrTimeStep > dMinTimeStep) {
                         /* Riduce il passo */
                         CurrStep = StepIntegrator::REPEATSTEP;
                         doublereal dOldCurrTimeStep = dCurrTimeStep;
                         dCurrTimeStep = pTSC->dGetNewStepTime(CurrStep, iStIter);
                         if (dCurrTimeStep < dOldCurrTimeStep) {
                                 DEBUGCOUT("Changing time step"
                                         " from " << dOldCurrTimeStep
                                         << " to " << dCurrTimeStep
                                         << " during first step after "
                                         << iStIter << " iterations"
                                         << std::endl);
                                 goto IfFirstStepIsToBeRepeated;
                         }
                 }
 
                 silent_cerr("Max iterations number "
                         << std::abs(pFirstRegularStep->GetIntegratorMaxIters())
                         << " has been reached during "
                         "first step, Time=" << dTime << "; "
                         << "TimeStep=" << dCurrTimeStep
                         << " cannot be reduced further; "
                         "aborting..." << std::endl);
                 pDM->Output(0, dTime, dCurrTimeStep, true);
 
                 throw Solver::ErrMaxIterations(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ErrSimulationDiverged) {
                 /*
                  * Mettere qui eventuali azioni speciali
                  * da intraprendere in caso di errore ...
                  */
 
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (LinearSolver::ErrFactor& err) {
                 /*
                  * Mettere qui eventuali azioni speciali
                  * da intraprendere in caso di errore ...
                  */
                 silent_cerr("First step failed because no pivot element "
                         "could be found for column " << err.iCol
                         << " (" << pDM->GetDofDescription(err.iCol) << "); "
                         "aborting..." << std::endl);
                 throw SimulationDiverged(MBDYN_EXCEPT_ARGS);
         }
         catch (NonlinearSolver::ConvergenceOnSolution) {
                 bSolConv = true;
         }
         catch (EndOfSimulation& eos) {
                 silent_cerr("Simulation ended during the first regular step:\n"
                         << eos.what() << "\n");
                 return false;
         }
 
         SAFEDELETE(pFirstRegularStep);
         pFirstRegularStep = 0;
 
         bOut = pDM->Output(lStep, dTime + dCurrTimeStep, dCurrTimeStep);
 
         /* Si fa dare l'std::ostream al file di output per il log */
         std::ostream& Out = pDM->GetOutFile();
 
         if (mbdyn_stop_at_end_of_time_step()) {
                 /* Fa l'output della soluzione al primo passo ed esce */
                 Out << "Interrupted during first step." << std::endl;
                 throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
         }
 
         if (outputMsg()) {
                 Out
                         << "Step " << lStep
                         << " " << dTime + dCurrTimeStep
                         << " " << dCurrTimeStep
                         << " " << iStIter
                         << " " << dTest
                         << " " << dSolTest
                         << " " << bSolConv
                         << " " << bOut
                         << std::endl;
         }
 
         bSolConv = false;
 
         dRefTimeStep = dCurrTimeStep;
         dTime += dRefTimeStep;
 
         dTotErr += dTest;
         iTotIter += iStIter;
 
         if (EigAn.bAnalysis
                 && EigAn.currAnalysis != EigAn.Analyses.end()
                 && *EigAn.currAnalysis <= dTime)
         {
                 std::vector<doublereal>::iterator i = std::find_if(EigAn.Analyses.begin(),
                         EigAn.Analyses.end(), bind2nd(std::greater<doublereal>(), dTime));
                 if (i != EigAn.Analyses.end()) {
                         EigAn.currAnalysis = --i;
                 }
                 Eig();
                 ++EigAn.currAnalysis;
         }
 
         if (pRTSolver) {
                 pRTSolver->Init();
         }
 
         /* Altri passi regolari */
         ASSERT(pRegularSteps != 0);
 
         /* Setup SolutionManager(s) */
         SetupSolmans(pRegularSteps->GetIntegratorNumUnknownStates(), true);
         pCurrStepIntegrator = pRegularSteps;
 
         eStatus = SOLVER_STATUS_STARTED;
 
         return true;
 }