LineSearchSolver Class Reference

#include <linesearch.h>

Inheritance diagram for LineSearchSolver:

Collaboration diagram for LineSearchSolver:

Classes
class	MaxIterations
class	ResidualNotDecreased
class	SlopeNotNegative
class	ZeroGradient

Public Member Functions
	LineSearchSolver (DataManager *pDM, const NonlinearSolverOptions &options, const struct LineSearchParameters &param)
	~LineSearchSolver (void)
void	Solve (const NonlinearProblem *pNLP, Solver *pS, const integer iMaxIter, const doublereal &Tol, integer &iIterCnt, doublereal &dErr, const doublereal &SolTol, doublereal &dSolErr)
Public Member Functions inherited from NonlinearSolver
	NonlinearSolver (const NonlinearSolverOptions &options)
virtual void	SetTest (NonlinearSolverTest *pr, NonlinearSolverTest *ps)
virtual	~NonlinearSolver (void)
virtual bool	MakeResTest (Solver *pS, const NonlinearProblem *pNLP, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest, doublereal &dTestDiff)
virtual integer	TotalAssembledJacobian (void)
virtual NonlinearSolverTest *	pGetResTest (void)
virtual NonlinearSolverTest *	pGetSolTest (void)
Public Member Functions inherited from SolverDiagnostics
	SolverDiagnostics (unsigned OF=OUTPUT_DEFAULT, DriveCaller *pOM=0)
virtual	~SolverDiagnostics (void)
void	SetNoOutput (void)
void	SetOutputMeter (DriveCaller *pOM)
void	SetOutputDriveHandler (const DriveHandler *pDH)
void	SetOutputFlags (unsigned OF)
void	AddOutputFlags (unsigned OF)
void	DelOutputFlags (unsigned OF)
MatrixHandler::Norm_t	GetCondMatNorm (void) const
bool	outputMeter (void) const
bool	outputIters (void) const
bool	outputRes (void) const
bool	outputSol (void) const
bool	outputJac (void) const
bool	outputStep (void) const
bool	outputBailout (void) const
bool	outputCounter (void) const
bool	outputMatrixConditionNumber (void) const
bool	outputSolverConditionNumber (void) const
bool	outputSolverConditionStat (void) const
bool	outputCPUTime (void) const
bool	outputMsg (void) const

Private Member Functions
void	LineSearch (doublereal stpmax, doublereal fold, doublereal &f, bool &check, const integer &iIterCnt)
void	Residual (doublereal &f, integer iIterCnt)
void	Jacobian ()
Private Member Functions inherited from LineSearchParameters
	LineSearchParameters ()

Private Attributes
VectorHandler *	pRes
VectorHandler *	pSol
MyVectorHandler	p
MyVectorHandler	g
MyVectorHandler	dXneg
const NonlinearProblem *	pNLP
Solver *	pS
const DataManager *const	pDM
Private Attributes inherited from LineSearchParameters
doublereal	dTolX
doublereal	dTolMin
integer	iMaxIterations
doublereal	dMaxStep
doublereal	dAlpha
doublereal	dLambdaMin
doublereal	dLambdaFactMin
doublereal	dDivergenceCheck
doublereal	dMinStepScale
unsigned	uFlags

Additional Inherited Members
Public Types inherited from NonlinearSolver
enum	Type {   UNKNOWN = -1, NEWTONRAPHSON, MATRIXFREE, LINESEARCH,   DEFAULT = NEWTONRAPHSON, LASTSOLVERTYPE }
Protected Types inherited from NonlinearSolver
enum	CPUTimeType { CPU_RESIDUAL, CPU_JACOBIAN, CPU_LINEAR_SOLVER, CPU_LAST_TYPE }
Protected Types inherited from SolverDiagnostics
enum	{   OUTPUT_NONE = 0x0000, OUTPUT_ITERS = 0x0001, OUTPUT_RES = 0x0002, OUTPUT_SOL = 0x0004,   OUTPUT_JAC = 0x0008, OUTPUT_BAILOUT = 0x0010, OUTPUT_MSG = 0x0020, OUTPUT_COUNTER = 0x0040,   OUTPUT_MAT_COND_NUM_1 = 0x0080, OUTPUT_MAT_COND_NUM_INF = 0x0100, OUTPUT_SOLVER_COND_NUM = 0x0200, OUTPUT_SOLVER_COND_STAT = 0x400,   OUTPUT_CPU_TIME = 0x800, OUTPUT_MAT_COND_NUM = OUTPUT_MAT_COND_NUM_1 | OUTPUT_MAT_COND_NUM_INF, OUTPUT_DEFAULT = OUTPUT_MSG, OUTPUT_STEP = (OUTPUT_ITERS | OUTPUT_RES | OUTPUT_SOL | OUTPUT_JAC | OUTPUT_MAT_COND_NUM | OUTPUT_SOLVER_COND_NUM),   OUTPUT_MASK = 0x07FF }
Protected Types inherited from NonlinearSolverOptions
enum	ScaleFlags { SCALE_ALGEBRAIC_EQUATIONS_NO = 0, SCALE_ALGEBRAIC_EQUATIONS_YES = 1 }
Protected Member Functions inherited from NonlinearSolver
virtual bool	MakeSolTest (Solver *pS, const VectorHandler &Vec, const doublereal &dTol, doublereal &dTest)
doublereal	dGetCondMax () const
doublereal	dGetCondMin () const
doublereal	dGetCondAvg () const
doublereal	dGetTimeCPU (CPUTimeType eType) const
void	AddCond (doublereal dCond)
void	AddTimeCPU (doublereal dTime, CPUTimeType eType)
Protected Member Functions inherited from NonlinearSolverOptions
	NonlinearSolverOptions (bool bHonorJacRequest=false, enum ScaleFlags eScaleFlags=SCALE_ALGEBRAIC_EQUATIONS_NO, doublereal dScaleAlgebraic=1.)
Protected Attributes inherited from NonlinearSolver
integer	Size
integer	TotJac
NonlinearSolverTest *	pResTest
NonlinearSolverTest *	pSolTest
Protected Attributes inherited from SolverDiagnostics
unsigned	OutputFlags
DriveCaller *	pOutputMeter
Protected Attributes inherited from NonlinearSolverOptions
bool	bHonorJacRequest
enum NonlinearSolverOptions::ScaleFlags	eScaleFlags
doublereal	dScaleAlgebraic
Private Types inherited from LineSearchParameters
enum	{   ZERO_GRADIENT_CONTINUE = 0x001, DIVERGENCE_CHECK = 0x002, ALGORITHM_CUBIC = 0x004, ALGORITHM_FACTOR = 0x008,   PRINT_CONVERGENCE_INFO = 0x010, SCALE_NEWTON_STEP = 0x020, RELATIVE_LAMBDA_MIN = 0x040, ABORT_AT_LAMBDA_MIN = 0x080,   VERBOSE_MODE = 0x100, NON_NEGATIVE_SLOPE_CONTINUE = 0x200, ALGORITHM = ALGORITHM_CUBIC | ALGORITHM_FACTOR }

Detailed Description

Definition at line 81 of file linesearch.h.

Constructor & Destructor Documentation

LineSearchSolver::LineSearchSolver	(	DataManager *	pDM,
		const NonlinearSolverOptions &	options,
		const struct LineSearchParameters &	param
	)

Definition at line 102 of file linesearch.cc.

References LineSearchParameters::ALGORITHM_CUBIC, LineSearchParameters::ALGORITHM_FACTOR, NonlinearSolverOptions::bHonorJacRequest, LineSearchParameters::dAlpha, LineSearchParameters::dDivergenceCheck, LineSearchParameters::DIVERGENCE_CHECK, LineSearchParameters::dLambdaFactMin, LineSearchParameters::dLambdaMin, LineSearchParameters::dMaxStep, LineSearchParameters::dMinStepScale, LineSearchParameters::dTolMin, LineSearchParameters::dTolX, LineSearchParameters::iMaxIterations, LineSearchParameters::PRINT_CONVERGENCE_INFO, LineSearchParameters::RELATIVE_LAMBDA_MIN, LineSearchParameters::SCALE_NEWTON_STEP, TRACE_FLAG, TRACE_VAR, LineSearchParameters::uFlags, and LineSearchParameters::ZERO_GRADIENT_CONTINUE.

: NonlinearSolver(options),
LineSearchParameters(param),
pRes(0),
pSol(0),
pNLP(0),
pS(0),
pDM(pDM)
{
        TRACE_VAR(dTolX);
        // TRACE_VAR(dTolF);
        TRACE_VAR(dTolMin);
        TRACE_VAR(iMaxIterations);
        TRACE_VAR(dMaxStep);
        TRACE_VAR(dAlpha);
        TRACE_VAR(dLambdaMin);
        TRACE_VAR(dLambdaFactMin);
        TRACE_VAR(dMinStepScale);
        TRACE_VAR(dDivergenceCheck);
        TRACE_VAR(bHonorJacRequest);
        TRACE_VAR(uFlags);
        TRACE_FLAG(uFlags, ALGORITHM_CUBIC);
        TRACE_FLAG(uFlags, ALGORITHM_FACTOR);
        TRACE_FLAG(uFlags, SCALE_NEWTON_STEP);
        TRACE_FLAG(uFlags, RELATIVE_LAMBDA_MIN);
        TRACE_FLAG(uFlags, ZERO_GRADIENT_CONTINUE);
        TRACE_FLAG(uFlags, DIVERGENCE_CHECK);
        TRACE_FLAG(uFlags, PRINT_CONVERGENCE_INFO);
}

LineSearchSolver::~LineSearchSolver

(

void

)

Definition at line 134 of file linesearch.cc.

References NO_OP.

{
        NO_OP;
}

Member Function Documentation

void LineSearchSolver::Jacobian ( )

private

Definition at line 163 of file linesearch.cc.

References MatrixHandler::ConditionNumber(), SolverDiagnostics::GetCondMatNorm(), NonlinearProblem::Jacobian(), SolutionManager::MatrInitialize(), SolutionManager::MatrReset(), MBDYN_EXCEPT_ARGS, SolverDiagnostics::outputJac(), SolverDiagnostics::outputMatrixConditionNumber(), Solver::pGetSolutionManager(), SolutionManager::pMatHdl(), pNLP, pS, and NonlinearSolver::TotJac.

Referenced by Solve().

{
        SolutionManager *const pSM = pS->pGetSolutionManager();

        const integer iMaxIterRebuild = 10;

        pSM->MatrReset();
        integer iIter = 0;

        do {
                try {
                        pNLP->Jacobian(pSM->pMatHdl());
                } catch (const MatrixHandler::ErrRebuildMatrix&) {
                        silent_cout("NewtonRaphsonSolver: "
                                        "rebuilding matrix..."
                                        << std::endl);

                        /* need to rebuild the matrix... */
                        pSM->MatrInitialize();
                        continue;
                }
                break;
        } while (++iIter < iMaxIterRebuild);

        if (iIter >= iMaxIterRebuild) {
                silent_cerr("Maximum number of iterations exceeded when rebuilding the Jacobian matrix" << std::endl);
                throw MatrixHandler::ErrRebuildMatrix(MBDYN_EXCEPT_ARGS);
        }

        TotJac++;

#ifdef USE_MPI
        if (!bParallel || MBDynComm.Get_rank() == 0)
#endif /* USE_MPI */
        {
                if (outputJac()) {
                        silent_cout("Jacobian:" << std::endl
                                        << *(pSM->pMatHdl()));
                }

                if (outputMatrixConditionNumber()) {
                        silent_cout(" cond=" << pSM->pMatHdl()->ConditionNumber(GetCondMatNorm()) << std::endl);
                }
        }
}

Here is the call graph for this function:

void LineSearchSolver::LineSearch ( doublereal  stpmax, doublereal  fold, doublereal &  f, bool &  check, const integer &  iIterCnt  )

private

Definition at line 210 of file linesearch.cc.

References a, LineSearchParameters::ALGORITHM_CUBIC, ASSERT, Solver::CheckTimeStepLimit(), LineSearchParameters::dAlpha, DataManager::dGetTime(), LineSearchParameters::dLambdaFactMin, LineSearchParameters::dLambdaMin, LineSearchParameters::dMinStepScale, LineSearchParameters::dTolX, dXneg, g, DataManager::GetpXCurr(), DataManager::GetpXPCurr(), LineSearchParameters::iMaxIterations, VectorHandler::InnerProd(), NonlinearSolver::MakeResTest(), MBDYN_EXCEPT_ARGS, mbdyn_stop_at_end_of_iteration(), LineSearchParameters::NON_NEGATIVE_SLOPE_CONTINUE, VectorHandler::Norm(), SolverDiagnostics::outputIters(), p, pDM, pNLP, pRes, LineSearchParameters::PRINT_CONVERGENCE_INFO, pS, pSol, LineSearchParameters::RELATIVE_LAMBDA_MIN, VectorHandler::Reset(), Residual(), VectorHandler::ScalarMul(), LineSearchParameters::SCALE_NEWTON_STEP, NonlinearSolver::Size, grad::sqrt(), TRACE, TRACE_VAR, LineSearchParameters::uFlags, NonlinearProblem::Update(), and LineSearchParameters::VERBOSE_MODE.

Referenced by Solve().

{
        check = false;

        if (uFlags & SCALE_NEWTON_STEP) {
                const doublereal dNormSol = pSol->Norm();

                if (dNormSol > stpmax) {
                        const doublereal dScale = std::max(dMinStepScale, stpmax / dNormSol);
                        if (uFlags & VERBOSE_MODE) {
                                silent_cerr("line search warning: "
                                        "Newton increment is reduced by factor " << dScale
                        << " at time step " << pDM->dGetTime()
                                        << " at iteration " << iIterCnt
                                        << " The time step is probably too large" << std::endl);
                        }
                        ASSERT(stpmax >= 0);
                        ASSERT(dScale <= 1.);
                        ASSERT(dScale > 0.);
                        *pSol *= dScale;
                }
        }

        p = *pSol; // save the Newton increment

        doublereal slope = g.InnerProd(p);

        TRACE_VAR(slope);

        doublereal dLambdaMinEff = -1.;

        if (slope >= 0.) {
                if (uFlags & VERBOSE_MODE) {
                        silent_cerr("line search warning: slope=" << slope << " >= 0"
                        " at time step " << pDM->dGetTime()
                                                << " at iteration " << iIterCnt << std::endl
                                                << "This could be a roundoff problem" << std::endl);
                }

                if (uFlags & NON_NEGATIVE_SLOPE_CONTINUE) {
                        // It seems to be a numerical problem.
                        // Line search may not work in this situation.
                        // Resort to the ordinary Newton Raphson algorithm
                        slope = 0.;
                        dLambdaMinEff = 1.;
                } else {
                        throw SlopeNotNegative(MBDYN_EXCEPT_ARGS);
                }
        }

        if (dLambdaMinEff < 0) {
                // dLambdaMinEff has to be detected
                if (uFlags & RELATIVE_LAMBDA_MIN) {
                        doublereal test = 0.;

                        for (integer i = 1; i <= Size; ++i) {
                                const doublereal temp = std::abs(p(i)) / std::max(std::max(std::abs((*pDM->GetpXCurr())(i)), std::abs((*pDM->GetpXPCurr())(i))), 1.);
                                if (temp > test) {
                                        test = temp;
                                }
                        }

                        dLambdaMinEff = std::max(dLambdaMin, std::min(dTolX / test, 1.));
                } else {
                        dLambdaMinEff = dLambdaMin;
                }
        }

        doublereal lambda = 1.;
        doublereal tmplam;
        doublereal lambda2 = 0; //f2, lambda2 initialized to silence g++ possible uninitialized warning
        doublereal f2 = 0.;

#ifdef DEBUG
        lambda2 = f2 = NAN;
#endif

        integer iIter = 0;

        TRACE_VAR(dLambdaMinEff);

        do {
                // FIXME: dLambdaMax not defined
                // ASSERT(lambda <= dLambdaMax);

                if (iIter > 0) {
                        TRACE("Start new step from Xold, XPold with lambda = " << lambda << " ..." << std::endl);
                        pNLP->Update(&dXneg); // restore the previous state

                        pSol->Reset();
                        pSol->ScalarMul(p, lambda); // scale the Newton increment by lambda
                }

                for (integer i = 1; i <= Size; ++i)
                        dXneg(i) = -(*pSol)(i); // save the increment; so we can restore the previous state

                TRACE("Update the nonlinear problem ... pSol->Norm()=" << pSol->Norm() << std::endl);

                pNLP->Update(pSol);
                Residual(f, iIter);

                ++iIter;

                TRACE("New value for f:" << f << std::endl);

                doublereal dErr = 0., dErrDiff = 0.;
                MakeResTest(pS, pNLP, *pRes, 0., dErr, dErrDiff);

                if (outputIters() && (uFlags & PRINT_CONVERGENCE_INFO)) {
#ifdef USE_MPI
                        if (!bParallel || MBDynComm.Get_rank() == 0)
#endif /* USE_MPI */
                        {
                                silent_cout("\t\tf(" << iIterCnt << ":" << iIter << ")=" << std::setw(12) << f
                                        << "\tErr=" << std::setw(12) << dErr
                                        << "\tlambda=" << std::setw(12) << lambda
                                        <<"\tslope=" << slope << std::endl);
                        }
                }

                pS->CheckTimeStepLimit(dErr, dErrDiff);

                if (f <= fold + dAlpha * lambda * slope) {
                        TRACE("Sufficient decrease in f: backtrack" << std::endl);
                        return;
                } else if (lambda <= dLambdaMinEff) {
                        TRACE("Checking for convergence: lambda=" << lambda << " < lambdaMin=" << dLambdaMinEff << std::endl);
                        check = true; // check for convergence
                        return;
                } else {
                        if (uFlags & ALGORITHM_CUBIC) {
                                TRACE("Calculate new value for alam ..." << std::endl);
                                if (lambda == 1.) {
                                        tmplam = -slope / (2 * (f - fold - slope));
                                        TRACE_VAR(tmplam);
                                } else {
                                        const doublereal rhs1 = f - fold - lambda * slope;
                                        const doublereal rhs2 = f2 - fold - lambda2 * slope;
                                        const doublereal a = (rhs1 / (lambda * lambda) - rhs2 / (lambda2 * lambda2)) / (lambda - lambda2);
                                        const doublereal b = (-lambda2 * rhs1 / (lambda * lambda) + lambda * rhs2 / (lambda2 * lambda2)) / (lambda - lambda2);

                                        if (a == 0.) {
                                                tmplam = -slope / (2. * b);
                                                TRACE_VAR(tmplam);
                                        } else {
                                                const doublereal disc = b * b - 3. * a * slope;

                                                if (disc < 0.) {
                                                        tmplam = 0.5 * lambda;
                                                        TRACE_VAR(tmplam);
                                                } else if (b <= 0.) {
                                                        tmplam = (-b + sqrt(disc)) / (3. * a);
                                                        TRACE_VAR(tmplam);
                                                } else {
                                                        tmplam = -slope / (b + sqrt(disc));
                                                        TRACE_VAR(tmplam);
                                                }

                                                if (tmplam > 0.5 * lambda) {
                                                        tmplam = 0.5 * lambda;
                                                        TRACE_VAR(tmplam);
                                                }
                                        }
                                }
                                lambda2 = lambda;
                                f2 = f;
                                TRACE_VAR(tmplam);
                                lambda = std::max(tmplam, dLambdaFactMin * lambda);
                        } else {
                                lambda *= dLambdaFactMin;
                        }
                }

                if (mbdyn_stop_at_end_of_iteration()) {
                        throw ErrInterrupted(MBDYN_EXCEPT_ARGS);
                }
        } while (iIter < iMaxIterations);

        if (uFlags & VERBOSE_MODE) {
                silent_cerr("line search warning: Maximum number of line search iterations="
                        << iMaxIterations
                        << " exceeded in line search at time step " << pDM->dGetTime() 
            << " at iteration " << iIterCnt << std::endl);
        }

        throw MaxIterations(MBDYN_EXCEPT_ARGS);
}

Here is the call graph for this function:

void LineSearchSolver::Residual ( doublereal &  f, integer  iIterCnt  )

private

Definition at line 140 of file linesearch.cc.

References VectorHandler::Dot(), SolverDiagnostics::outputRes(), pNLP, pRes, Solver::PrintResidual(), pS, VectorHandler::Reset(), and NonlinearProblem::Residual().

Referenced by LineSearch(), and Solve().

{
#ifdef USE_EXTERNAL
        SendExternal();
#endif /* USE_EXTERNAL */
        
        pRes->Reset();

        try {
        pNLP->Residual(pRes);
        }
        catch (const SolutionDataManager::ChangedEquationStructure&) {
                // The Jacobian will be rebuilt anyway
        }

        if (outputRes()) {
                pS->PrintResidual(*pRes, iIterCnt);
        }

        f = 0.5 * pRes->Dot();
}

Here is the call graph for this function:

void LineSearchSolver::Solve ( const NonlinearProblem *  pNLP, Solver *  pS, const integer  iMaxIter, const doublereal &  Tol, integer &  iIterCnt, doublereal &  dErr, const doublereal &  SolTol, doublereal &  dSolErr  )

virtual

Implements NonlinearSolver.

Definition at line 400 of file linesearch.cc.

References LineSearchParameters::ABORT_AT_LAMBDA_MIN, NonlinearSolver::AddCond(), NonlinearSolver::AddTimeCPU(), ASSERT, SolutionManager::bGetConditionNumber(), Solver::CheckTimeStepLimit(), NonlinearSolver::CPU_JACOBIAN, NonlinearSolver::CPU_LINEAR_SOLVER, NonlinearSolver::CPU_RESIDUAL, LineSearchParameters::dDivergenceCheck, NonlinearSolver::dGetCondAvg(), NonlinearSolver::dGetCondMax(), NonlinearSolver::dGetCondMin(), DataManager::dGetTime(), NonlinearSolver::dGetTimeCPU(), LineSearchParameters::DIVERGENCE_CHECK, LineSearchParameters::dMaxStep, VectorHandler::Dot(), LineSearchParameters::dTolMin, dXneg, g, DataManager::GetpXCurr(), DataManager::GetpXPCurr(), MatrixHandler::iGetNumCols(), MatrixHandler::iGetNumRows(), VectorHandler::iGetSize(), MyVectorHandler::iGetSize(), Jacobian(), LineSearch(), NonlinearSolver::MakeResTest(), NonlinearSolver::MakeSolTest(), MatrixHandler::MatTVecDecMul(), mbdyn_clock_time(), MBDYN_EXCEPT_ARGS, SolverDiagnostics::outputBailout(), SolverDiagnostics::outputCPUTime(), SolverDiagnostics::outputIters(), SolverDiagnostics::outputSol(), SolverDiagnostics::outputSolverConditionNumber(), SolverDiagnostics::outputSolverConditionStat(), p, pDM, Solver::pGetSolutionManager(), SolutionManager::pMatHdl(), pNLP, pRes, SolutionManager::pResHdl(), LineSearchParameters::PRINT_CONVERGENCE_INFO, Solver::PrintResidual(), Solver::PrintSolution(), pS, pSol, SolutionManager::pSolHdl(), MyVectorHandler::Reset(), Residual(), MyVectorHandler::Resize(), LineSearchParameters::SCALE_NEWTON_STEP, NonlinearSolver::Size, SolutionManager::Solve(), grad::sqrt(), TRACE, TRACE_VAR, LineSearchParameters::uFlags, LineSearchParameters::VERBOSE_MODE, and LineSearchParameters::ZERO_GRADIENT_CONTINUE.

{
        ASSERT(pS != NULL);
        SolutionManager *const pSM = pS->pGetSolutionManager();
        this->pS = pS;
        pNLP = pNonLinProblem;
        pRes = pSM->pResHdl();
        pSol = pSM->pSolHdl();
        Size = pRes->iGetSize();

        if (g.iGetSize() != Size) {
                TRACE("Resize temporary vectors ..." << std::endl);
                g.Resize(Size);
                p.Resize(Size);
                dXneg.Resize(Size);
        }

        bool check = false;
        iIterCnt = 0;
        dSolErr = 0.;
        dErr = 0.;
        doublereal dErrDiff = 0.;
        doublereal f;
        Residual(f, iIterCnt);

        TRACE("\t\tf(0) = " << f << std::endl);

        bool bTest = MakeResTest(pS, pNLP, *pRes, 1e-2 * Tol, dErr, dErrDiff); // use a more stringent test for the first iteration
        doublereal dPrevErr = std::numeric_limits<doublereal>::max(); // disable error test for the first iteration
        doublereal dErrFactor = 1.;

        if (outputIters()) {
#ifdef USE_MPI
                if (!bParallel || MBDynComm.Get_rank() == 0)
#endif /* USE_MPI */
                {
                        silent_cout("\tIteration(" << iIterCnt << ") " << std::setw(12) << dErr);

                        if (uFlags & PRINT_CONVERGENCE_INFO) {
                                silent_cout(" f=" << std::setw(12) << f);
                        }

                        silent_cout(std::endl);
                }
        }

        pS->CheckTimeStepLimit(dErr, dErrDiff);

        if (bTest) {
                return;
        }

        const doublereal stpmax = (uFlags & SCALE_NEWTON_STEP)
                ? dMaxStep * std::max(sqrt(pDM->GetpXPCurr()->Dot() + pDM->GetpXCurr()->Dot()), static_cast<doublereal>(Size))
                : std::numeric_limits<doublereal>::max();

        TRACE_VAR(stpmax);

        doublereal dStartTimeCPU, dEndTimeCPU, dJacobianCPU, dResidualCPU, dLinSolveCPU;
        
        while (true) {
                if (outputCPUTime()) {
                        dStartTimeCPU = mbdyn_clock_time();
                }

                Jacobian();

                ASSERT(pSM->pMatHdl()->iGetNumCols() == Size);
                ASSERT(pSM->pMatHdl()->iGetNumCols() == pRes->iGetSize());
                ASSERT(pSM->pMatHdl()->iGetNumRows() == pRes->iGetSize());

                g.Reset();
                pSM->pMatHdl()->MatTVecDecMul(g, *pRes); // compute gradient g = \nabla f = fjac^T \, fvec = -Jac^T \, pRes

                if (outputCPUTime()) {
                        dEndTimeCPU = mbdyn_clock_time();
                        dJacobianCPU = dEndTimeCPU - dStartTimeCPU;
                        AddTimeCPU(dJacobianCPU, CPU_JACOBIAN);
                        dStartTimeCPU = dEndTimeCPU;
                }

                const doublereal fold = f;

                pSM->Solve();

                if (outputCPUTime()) {
                        dEndTimeCPU = mbdyn_clock_time();
                        dLinSolveCPU = dEndTimeCPU - dStartTimeCPU;
                        AddTimeCPU(dLinSolveCPU, CPU_LINEAR_SOLVER);
                }

                if (outputSol()) {
                        pS->PrintSolution(*pSol, iIterCnt);
                }

                if (outputCPUTime()) {
                        dStartTimeCPU = mbdyn_clock_time();
                }

                LineSearch(stpmax, fold, f, check, iIterCnt);

                bTest = MakeResTest(pS, pNLP, *pRes, Tol, dErr, dErrDiff);

                if (outputCPUTime()) {
                        dEndTimeCPU = mbdyn_clock_time();
                        dResidualCPU = dEndTimeCPU - dStartTimeCPU;
                        AddTimeCPU(dResidualCPU, CPU_RESIDUAL);
                }

                if (iIterCnt > 0) {
                        dErrFactor *= dErr / dPrevErr;
                }

                iIterCnt++;

#ifdef USE_MPI
                if (!bParallel || MBDynComm.Get_rank() == 0)
#endif /* USE_MPI */
                {
                        if (outputIters() || outputSolverConditionNumber()) {
                                if (outputIters()) {
                                        silent_cout("\tIteration(" << iIterCnt << ") " << std::setw(12) << dErr << " J");
                                }

                                if (outputSolverConditionNumber()) {
                                        silent_cout(" cond=");
                                        doublereal dCond;
                                        if (pSM->bGetConditionNumber(dCond)) {
                                                silent_cout(dCond);
                                                if (outputSolverConditionStat()) {
                                                        AddCond(dCond);
                                                        silent_cout(" " << dGetCondMin() << " " << dGetCondMax() << " " << dGetCondAvg());
                                                }
                                        } else {
                                                silent_cout("NA");
                                        }
                                }

                                if (outputCPUTime()) {
                                        silent_cout(" CPU:" << dResidualCPU << "/" << dGetTimeCPU(CPU_RESIDUAL)
                                                << "+" << dJacobianCPU << "/" << dGetTimeCPU(CPU_JACOBIAN)
                                                << "+" << dLinSolveCPU << "/" << dGetTimeCPU(CPU_LINEAR_SOLVER));
                                }

                                if (uFlags & PRINT_CONVERGENCE_INFO) {
                                        silent_cout(" f=" << std::setw(12) << f);

                                        if (iIterCnt > 1) {
                                                silent_cout(" Err(n-1)=" << std::setw(12) << dPrevErr
                                                                 << " Err(n)/Err(n-1)=" << std::setw(12) << dErr / dPrevErr
                                                                 << " Err(n)/Err(1)=" << std::setw(12) << dErrFactor);
                                        }
                                }

                                if (outputIters()
                                                || outputSolverConditionNumber()
                                                || (uFlags & PRINT_CONVERGENCE_INFO)) {
                                        silent_cout(std::endl);
                                }
                        }
                }

                if (bTest) {
                        return;
                }

                bTest = MakeSolTest(pS, *pSol, SolTol, dSolErr);

                if (outputIters()) {
#ifdef USE_MPI
                        if (!bParallel || MBDynComm.Get_rank() == 0)
#endif /* USE_MPI */
                        {
                                silent_cout("\t\tSolErr "
                                        << dSolErr << std::endl);
                        }
                }

                if (bTest) {
                        if (uFlags & VERBOSE_MODE) {
                                silent_cerr("line search warning: Convergence on solution "
                    "at time step " << pDM->dGetTime()
                                        << "at iteration " << iIterCnt
                                        << "\tErr(n)=" << dErr);

                                if (iIterCnt >= 1) {
                                        silent_cerr("\tErr(n-1)=" << dPrevErr
                                        << "\tErr(n)/Err(n-1)=" << dErr/dPrevErr
                                        << "\tErr(n)/Err(1) " << dErrFactor);
                                }

                                silent_cerr(std::endl);
                        }
                        throw ConvergenceOnSolution(MBDYN_EXCEPT_ARGS);
                }

                if (check) { // lambda <= dLambdaMinEff: check for gradient zero
                        doublereal test = 0.;
                        const doublereal den = std::max(f, 0.5 * Size);

                        for (integer i = 1; i <= Size; ++i) {
                                const doublereal absX = std::max(std::abs((*pDM->GetpXCurr())(i)),
                                                                                                 std::abs((*pDM->GetpXPCurr())(i)));
                                const doublereal temp = std::abs(g(i)) * std::max(absX, 1.) / den;
                                if (temp > test) {
                                        test = temp;
                                }
                        }

                        if (test < dTolMin) { // we are at a local minimum of the function f
                                if (uFlags & VERBOSE_MODE) {
                                        silent_cerr("line search warning: Zero gradient detected at time step "
                        << pDM->dGetTime() << " at iteration "
                                                << iIterCnt << " (spurious convergence) test=" << test
                                                << " < tol=" << dTolMin
                                                << "\tErr(n)=" << dErr
                                                << " > Tol = " << Tol << std::endl);
                                }

                                if (uFlags & ZERO_GRADIENT_CONTINUE) {
                                        return;
                                } else {
                                        throw ZeroGradient(MBDYN_EXCEPT_ARGS);
                                }
                        } else {
                                if (uFlags & VERBOSE_MODE) {
                                        silent_cerr("line search warning: lambda min"
                                                " has been reached at time step " << pDM->dGetTime() 
                        << " at iteration " << iIterCnt
                                                << " but the gradient is not zero" << std::endl);
                                }

                                if (uFlags & ABORT_AT_LAMBDA_MIN) {
                                        throw NoConvergence(MBDYN_EXCEPT_ARGS);
                                }
                        }
                } else if (dErr > dPrevErr) { // f <= fold + dAlpha * lambda * slope
                        // We should not get here since f was decreased also Err should be decreased
                        // If not it could be a numerical problem (e.g. the tolerance could be too small)
                        if (uFlags & VERBOSE_MODE) {
                                silent_cerr("line search warning: f has been reduced during line search but the residual could not be reduced at time step "
                                                << pDM->dGetTime()
                                                << " at iteration " << iIterCnt << std::endl);
                        }

                        // Do not throw an exception here because if we specify for example
                        //              tolerance: <<Tol>>, test, minmax;
                        // it could happen that the the norm of the residual vector is decreased
                        // but the maximum residual is not!

                        // In any case we will check for divergence if DIVERGENCE_CHECK is enabled.
                }

                if (dErrFactor > dDivergenceCheck) {
                        if (uFlags & DIVERGENCE_CHECK) {
                                if (uFlags & VERBOSE_MODE) {
                                        silent_cerr("line search warning: The residual could not be decreased"
                                                " at time step " << pDM->dGetTime() 
                                << " at iteration " << iIterCnt << std::endl);

                                        if (iIterCnt > 1) {
                                                silent_cerr("Err(n-1)=" << dPrevErr
                                                                        << "\tErr(n)=" << dErr
                                                                        << "\tErr(n)/Err(n-1)=" << dErr/dPrevErr
                                                                        << "\tErr(n)/Err(1)=" << dErrFactor << std::endl);
                                        }
                                }
                                throw ResidualNotDecreased(MBDYN_EXCEPT_ARGS);
                        }
                }

                if (iIterCnt >= std::abs(iMaxIter)) {
                        if (iMaxIter < 0 && dErrFactor < 1.) {
                                return;
                        }
                        if (outputBailout()) {
                                pS->PrintResidual(*pRes, iIterCnt);
                        }
                        throw NoConvergence(MBDYN_EXCEPT_ARGS);
                }

                dPrevErr = dErr;
        }
}

Here is the call graph for this function:

Member Data Documentation

MyVectorHandler LineSearchSolver::dXneg

private

Definition at line 87 of file linesearch.h.

Referenced by LineSearch(), and Solve().

MyVectorHandler LineSearchSolver::g

private

Definition at line 86 of file linesearch.h.

Referenced by LineSearch(), and Solve().

MyVectorHandler LineSearchSolver::p

private

Definition at line 85 of file linesearch.h.

Referenced by LineSearch(), and Solve().

const DataManager* const LineSearchSolver::pDM

private

Definition at line 90 of file linesearch.h.

Referenced by LineSearch(), and Solve().

const NonlinearProblem* LineSearchSolver::pNLP

private

Definition at line 88 of file linesearch.h.

Referenced by Jacobian(), LineSearch(), Residual(), and Solve().

VectorHandler* LineSearchSolver::pRes

private

Definition at line 83 of file linesearch.h.

Referenced by LineSearch(), Residual(), and Solve().

Solver* LineSearchSolver::pS

private

Definition at line 89 of file linesearch.h.

Referenced by Jacobian(), LineSearch(), Residual(), and Solve().

VectorHandler* LineSearchSolver::pSol

private

Definition at line 84 of file linesearch.h.

Referenced by LineSearch(), and Solve().

---

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

• linesearch.h
• linesearch.cc