Solver Class Reference

#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.

References Solver::EigenAnalysis::Analyses, DataManager::AssJac(), buf, Solver::EigenAnalysis::currAnalysis, DEBUGCOUT, DEBUGCOUTFNAME, Solver::EigenAnalysis::dParam, dTime, Solver::EigenAnalysis::EIG_OUTPUT, Solver::EigenAnalysis::EIG_OUTPUT_FULL_MATRICES, Solver::EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES, Solver::EigenAnalysis::EIG_USE_ARPACK, Solver::EigenAnalysis::EIG_USE_JDQZ, Solver::EigenAnalysis::EIG_USE_LAPACK, Solver::EigenAnalysis::EIG_USE_MASK, EigAn, Solver::EigenAnalysis::iFNameFormat, Solver::EigenAnalysis::iFNameWidth, Solver::EigenAnalysis::iMatrixPrecision, iNumDofs, Solver::EigenAnalysis::iResultsPrecision, DataManager::OutputEigClose(), DataManager::OutputEigFullMatrices(), DataManager::OutputEigNaiveMatrices(), DataManager::OutputEigOpen(), DataManager::OutputEigParams(), DataManager::OutputEigSparseMatrices(), SolutionManager::pMatHdl(), MatrixHandler::Reset(), SAFEDELETE, SAFENEWWITHCONSTRUCTOR, and Solver::EigenAnalysis::uFlags.

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.

References LinSol::AddSolverFlags(), StepIntegrator::Advance(), AFTER_ASSEMBLY, AFTER_DERIVATIVES, AFTER_DUMMY_STEPS, AFTER_INPUT, AllocateNonlinearSolver(), Solver::EigenAnalysis::Analyses, ASSERT, Solver::EigenAnalysis::bAnalysis, DataManager::BeforePredict(), bOut, bOutputCounter, bParallel, bScale, bSolConv, IncludeParser::Close(), SchurDataManager::CreatePartition(), Solver::EigenAnalysis::currAnalysis, CurrLinearSolver, dCurrTimeStep, dDerivativesCoef, dDummyStepsRatio, DEBUG_LEVEL, DEBUG_LEVEL_MATCH, DEBUGCOUT, DEBUGCOUTFNAME, DEBUGLCOUT, dInitialTime, dInitialTimeStep, dRefTimeStep, dSolTest, dTest, dTime, dTotErr, eAbortAfter, Eig(), EigAn, Solver::EigenAnalysis::EIGAN_WIDTH_COMPUTE, EMPTY, environ, eStatus, Flip(), DataManager::GetDofDescription(), SchurDataManager::GetDofsList(), StepIntegrator::GetIntegratorNumPreviousStates(), StepIntegrator::GetIntegratorNumUnknownStates(), DataManager::GetLogFile(), HighParser::GetMathParser(), DataManager::GetOutFile(), LinSol::GetSolverFlags(), LinSol::GetSolverName(), MathParser::GetSymbolTable(), SchurDataManager::HowManyDofs(), HP, LinearSolver::ErrFactor::iCol, iDummyStepsNumber, Solver::EigenAnalysis::iFNameWidth, DataManager::iGetNumDofs(), SchurDataManager::INTERNAL, iNumDofs, iNumIntDofs, iNumLocDofs, iNumPreviousVectors, iStIter, iTotIter, iUnkStates, DataManager::LinkToSolution(), SchurDataManager::LOCAL, grad::log(), lStep, MBDYN_EXCEPT_ARGS, mbdyn_signal_init(), mbdyn_stop_at_end_of_time_step(), NonlinearSolverTest::MINMAX, MYDEBUG_DERIVATIVES, MYDEBUG_FSTEPS, MYDEBUG_JAC, MYDEBUG_MEM, MYDEBUG_PRED, MYDEBUG_RESIDUAL, MYDEBUG_SOL, StepIntegrator::NEWSTEP, NonlinearSolverTest::NONE, NonlinearSolverTest::NORM, DataManager::Output(), SolverDiagnostics::OUTPUT_JAC, SolverDiagnostics::OUTPUT_RES, SolverDiagnostics::OUTPUT_SOL, SolverDiagnostics::outputCounter(), outputCounterPostfix, outputCounterPrefix, DataManager::OutputEigPrepare(), SolverDiagnostics::OutputFlags, SolverDiagnostics::outputMsg(), DataManager::OutputPrepare(), SolverDiagnostics::outputStep(), StepIntegrator::OutputTypes(), DriveCaller::pCopy(), pCurrStepIntegrator, pDerivativeSteps, pDM, pDofs, pDummySteps, pdWorkSpace, pFirstDummyStep, pFirstRegularStep, SchurDataManager::pGetDofsList(), DataManager::pGetDrvHdl(), pIntDofs, pLocDofs, pNLS, SolverDiagnostics::pOutputMeter, pRegularSteps, pRTSolver, pSDM, pTSC, pX, pXPrime, qX, qXPrime, ReadData(), MyVectorHandler::Reset(), VectorHandler::ResizeReset(), ResTest, SAFEDELETE, SAFENEW, SAFENEWARR, SAFENEWWITHCONSTRUCTOR, Scale, StepIntegrator::SetDataManager(), TimeStepControl::SetDriveHandler(), StepIntegrator::SetDriveHandler(), DriveCaller::SetDrvHdl(), SolverDiagnostics::SetOutputDriveHandler(), SolverDiagnostics::SetOutputFlags(), SolverDiagnostics::SetOutputMeter(), NonlinearSolverTestScale::SetScale(), DataManager::SetScale(), LinSol::SetSolver(), NonlinearSolver::SetTest(), DataManager::SetTime(), RTSolverBase::Setup(), SetupSolmans(), DataManager::SetValue(), sInputFileName, SolTest, LinSol::SOLVER_FLAGS_ALLOWS_MT_ASS, SOLVER_STATUS_PREPARED, SOLVER_STATUS_UNINITIALIZED, sOutputFileName, SchurDataManager::TOTAL, LinSol::UMFPACK_SOLVER, MBDynErrBase::what(), and LinSol::Y12_SOLVER.

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.

References LineSearchParameters::ABORT_AT_LAMBDA_MIN, SolverDiagnostics::AddOutputFlags(), AFTER_ASSEMBLY, AFTER_DERIVATIVES, AFTER_DUMMY_STEPS, AFTER_INPUT, LineSearchParameters::ALGORITHM_CUBIC, LineSearchParameters::ALGORITHM_FACTOR, Solver::EigenAnalysis::Analyses, Solver::EigenAnalysis::arpack, ASSERT, Solver::EigenAnalysis::bAnalysis, NonlinearSolverOptions::bHonorJacRequest, MatrixFreeSolver::BICGSTAB, bKeepJac, bParallel, bScale, bTrueNewtonRaphson, Solver::EigenAnalysis::currAnalysis, CurrIntSolver, CurrLinearSolver, LineSearchParameters::dAlpha, dDefaultDerivativesCoefficient, dDefaultDummyStepsTolerance, dDefaultTol, dDerivativesCoef, LineSearchParameters::dDivergenceCheck, dDummyStepsRatio, DEBUGCOUTFNAME, DEBUGLCOUT, MatrixFreeSolver::DEFAULT, NonlinearSolver::DEFAULT, SolverDiagnostics::DelOutputFlags(), dFinalTime, DriveOwner::dGet(), dGetInitialMaxTimeStep(), dInitialTime, dInitialTimeStep, dIterertiveEtaMax, dIterertiveTau, dIterTol, LineSearchParameters::DIVERGENCE_CHECK, LineSearchParameters::dLambdaFactMin, LineSearchParameters::dLambdaMin, Solver::EigenAnalysis::dLowerFreq, dMaxResidual, dMaxResidualDiff, LineSearchParameters::dMaxStep, LineSearchParameters::dMinStepScale, dMinTimeStep, Solver::EigenAnalysis::dParam, NonlinearSolverOptions::dScaleAlgebraic, Solver::EigenAnalysis::ARPACK::dTOL, LineSearchParameters::dTolMin, LineSearchParameters::dTolX, DummyType, Solver::EigenAnalysis::dUpperFreq, eAbortAfter, Solver::EigenAnalysis::EIG_BALANCE, Solver::EigenAnalysis::EIG_OUTPUT_EIGENVECTORS, Solver::EigenAnalysis::EIG_OUTPUT_FULL_MATRICES, Solver::EigenAnalysis::EIG_OUTPUT_GEOMETRY, Solver::EigenAnalysis::EIG_OUTPUT_MATRICES, Solver::EigenAnalysis::EIG_OUTPUT_SPARSE_MATRICES, Solver::EigenAnalysis::EIG_PERMUTE, Solver::EigenAnalysis::EIG_SCALE, Solver::EigenAnalysis::EIG_SOLVE, Solver::EigenAnalysis::EIG_USE_ARPACK, Solver::EigenAnalysis::EIG_USE_JDQZ, Solver::EigenAnalysis::EIG_USE_LAPACK, Solver::EigenAnalysis::EIG_USE_MASK, EigAn, Solver::EigenAnalysis::EIGAN_WIDTH_COMPUTE, END, Solver::EigenAnalysis::JDQZ::eps, NonlinearSolverOptions::eScaleFlags, eTimeStepLimit, Preconditioner::FULLJACOBIANMATRIX, HighParser::ErrValueOutOfRange< T >::Get(), get_nprocs(), HighParser::GetDescription(), MBDynParser::GetDriveCaller(), HighParser::GetInt(), IncludeParser::GetLineData(), HighParser::GetReal(), LinSol::GetSolverName(), HighParser::GetStringWithDelims(), HighParser::GetWord(), HighParser::GetYesNoOrBool(), MatrixFreeSolver::GMRES, iDefaultDummyStepsNumber, iDefaultIterationsBeforeAssembly, iDefaultMaxIterations, iDummyStepsNumber, Solver::EigenAnalysis::iFNameFormat, Solver::EigenAnalysis::iFNameWidth, iIterationsBeforeAssembly, iIterativeMaxSteps, Solver::EigenAnalysis::iMatrixPrecision, LineSearchParameters::iMaxIterations, iMaxIterations, Solver::EigenAnalysis::ARPACK::iNCV, Solver::EigenAnalysis::ARPACK::iNEV, INT_CRANKNICOLSON, INT_HOPE, INT_IMPLICITEULER, INT_MS2, INT_THIRDORDER, INT_UNKNOWN, iPrecondSteps, Solver::EigenAnalysis::iResultsPrecision, HighParser::IsArg(), HighParser::IsKeyWord(), Solver::EigenAnalysis::jdqz, Solver::EigenAnalysis::JDQZ::jmax, Solver::EigenAnalysis::JDQZ::jmin, Solver::EigenAnalysis::JDQZ::kmax, LASTKEYWORD, NonlinearSolver::LINESEARCH, LineSearch, NonlinearSolver::MATRIXFREE, MaxTimeStep, MBDYN_EXCEPT_ARGS, MFSolverType, NonlinearSolverTest::MINMAX, MYDEBUG_INPUT, NonlinearSolver::NEWTONRAPHSON, LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE, NonlinearSolverTest::NONE, NonlinearSolverType, NonlinearSolverTest::NORM, OUTPUT, SolverDiagnostics::OUTPUT_BAILOUT, SolverDiagnostics::OUTPUT_COUNTER, SolverDiagnostics::OUTPUT_CPU_TIME, SolverDiagnostics::OUTPUT_DEFAULT, SolverDiagnostics::OUTPUT_ITERS, SolverDiagnostics::OUTPUT_JAC, SolverDiagnostics::OUTPUT_MAT_COND_NUM, SolverDiagnostics::OUTPUT_MAT_COND_NUM_1, SolverDiagnostics::OUTPUT_MAT_COND_NUM_INF, SolverDiagnostics::OUTPUT_MSG, SolverDiagnostics::OUTPUT_NONE, SolverDiagnostics::OUTPUT_RES, SolverDiagnostics::OUTPUT_SOL, SolverDiagnostics::OUTPUT_SOLVER_COND_NUM, SolverDiagnostics::OUTPUT_SOLVER_COND_STAT, DriveCaller::pCopy(), PcType, pDerivativeSteps, pDummySteps, pFirstDummyStep, pFirstRegularStep, DriveOwner::pGetDriveCaller(), pRegularSteps, pRhoAlgebraicDummy, pRhoAlgebraicRegular, pRhoDummy, pRhoRegular, LineSearchParameters::PRINT_CONVERGENCE_INFO, pRTSolver, pTSC, ReadLinSol(), ReadRTSolver(), ReadTimeStepData(), RegularType, LineSearchParameters::RELATIVE_LAMBDA_MIN, ResTest, SAFENEW, SAFENEWWITHCONSTRUCTOR, NonlinearSolverOptions::SCALE_ALGEBRAIC_EQUATIONS_NO, NonlinearSolverOptions::SCALE_ALGEBRAIC_EQUATIONS_YES, LineSearchParameters::SCALE_NEWTON_STEP, DriveOwner::Set(), LinSol::SetNumThreads(), SolverDiagnostics::SetOutputFlags(), SolverDiagnostics::SetOutputMeter(), SolTest, TS_HARD_LIMIT, TS_SOFT_LIMIT, LineSearchParameters::uFlags, Solver::EigenAnalysis::uFlags, LineSearchParameters::VERBOSE_MODE, MBDynErrBase::what(), and LineSearchParameters::ZERO_GRADIENT_CONTINUE.

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()) {