13 #ifndef EIGEN_HYBRIDNONLINEARSOLVER_H
14 #define EIGEN_HYBRIDNONLINEARSOLVER_H
20 namespace HybridNonLinearSolverSpace {
44 template<
typename FunctorType,
typename Scalar=
double>
122 template<
typename FunctorType,
typename Scalar>
132 if (
n <= 0 || tol < 0.)
136 parameters.maxfev = 100*(
n+1);
137 parameters.xtol = tol;
138 diag.setConstant(
n, 1.);
139 useExternalScaling =
true;
143 template<
typename FunctorType,
typename Scalar>
149 wa1.resize(
n); wa2.resize(
n); wa3.resize(
n); wa4.resize(
n);
153 if (!useExternalScaling)
155 eigen_assert( (!useExternalScaling || diag.size()==
n) &&
"When useExternalScaling is set, the caller must provide a valid 'diag'");
162 if (
n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
164 if (useExternalScaling)
172 if ( functor(
x, fvec) < 0)
174 fnorm = fvec.stableNorm();
186 template<
typename FunctorType,
typename Scalar>
195 std::vector<JacobiRotation<Scalar> > v_givens(
n), w_givens(
n);
200 if ( functor.df(
x, fjac) < 0)
204 wa2 = fjac.colwise().blueNorm();
209 if (!useExternalScaling)
210 for (
j = 0;
j <
n; ++
j)
211 diag[
j] = (wa2[
j]==0.) ? 1. : wa2[
j];
215 xnorm = diag.cwiseProduct(
x).stableNorm();
216 delta = parameters.factor * xnorm;
218 delta = parameters.factor;
234 if (!useExternalScaling)
235 diag = diag.cwiseMax(wa2);
239 internal::dogleg<Scalar>(R, diag, qtf, delta, wa1);
244 pnorm = diag.cwiseProduct(wa1).stableNorm();
251 if ( functor(wa2, wa4) < 0)
254 fnorm1 = wa4.stableNorm();
262 wa3 = R.template triangularView<Upper>()*wa1 + qtf;
263 temp = wa3.stableNorm();
271 ratio = actred / prered;
274 if (ratio < Scalar(.1)) {
277 delta = Scalar(.5) * delta;
281 if (ratio >= Scalar(.5) || ncsuc > 1)
282 delta = (
std::max)(delta, pnorm / Scalar(.5));
283 if (
abs(ratio - 1.) <= Scalar(.1)) {
284 delta = pnorm / Scalar(.5);
289 if (ratio >= Scalar(1
e-4)) {
292 wa2 = diag.cwiseProduct(
x);
294 xnorm = wa2.stableNorm();
301 if (actred >= Scalar(.001))
305 if (actred >= Scalar(.1))
309 if (delta <= parameters.xtol * xnorm || fnorm == 0.)
313 if (nfev >= parameters.maxfev)
328 wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm );
329 wa2 = fjac.transpose() * wa4;
330 if (ratio >= Scalar(1
e-4))
332 wa2 = (wa2-wa3)/pnorm;
335 internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing);
336 internal::r1mpyq<Scalar>(
n,
n, fjac.data(), v_givens, w_givens);
337 internal::r1mpyq<Scalar>(1,
n, qtf.data(), v_givens, w_givens);
344 template<
typename FunctorType,
typename Scalar>
352 status = solveOneStep(
x);
358 template<
typename FunctorType,
typename Scalar>
368 if (
n <= 0 || tol < 0.)
372 parameters.maxfev = 200*(
n+1);
373 parameters.xtol = tol;
375 diag.setConstant(
n, 1.);
376 useExternalScaling =
true;
377 return solveNumericalDiff(
x);
380 template<
typename FunctorType,
typename Scalar>
386 if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals=
n-1;
387 if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals=
n-1;
389 wa1.resize(
n); wa2.resize(
n); wa3.resize(
n); wa4.resize(
n);
393 if (!useExternalScaling)
395 eigen_assert( (!useExternalScaling || diag.size()==
n) &&
"When useExternalScaling is set, the caller must provide a valid 'diag'");
402 if (
n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
404 if (useExternalScaling)
412 if ( functor(
x, fvec) < 0)
414 fnorm = fvec.stableNorm();
426 template<
typename FunctorType,
typename Scalar>
436 std::vector<JacobiRotation<Scalar> > v_givens(
n), w_givens(
n);
439 if (parameters.nb_of_subdiagonals<0) parameters.nb_of_subdiagonals=
n-1;
440 if (parameters.nb_of_superdiagonals<0) parameters.nb_of_superdiagonals=
n-1;
443 if (
internal::fdjac1(functor,
x, fvec, fjac, parameters.nb_of_subdiagonals, parameters.nb_of_superdiagonals, parameters.epsfcn) <0)
445 nfev += (
std::min)(parameters.nb_of_subdiagonals+parameters.nb_of_superdiagonals+ 1,
n);
447 wa2 = fjac.colwise().blueNorm();
452 if (!useExternalScaling)
453 for (
j = 0;
j <
n; ++
j)
454 diag[
j] = (wa2[
j]==0.) ? 1. : wa2[
j];
458 xnorm = diag.cwiseProduct(
x).stableNorm();
459 delta = parameters.factor * xnorm;
461 delta = parameters.factor;
477 if (!useExternalScaling)
478 diag = diag.cwiseMax(wa2);
482 internal::dogleg<Scalar>(R, diag, qtf, delta, wa1);
487 pnorm = diag.cwiseProduct(wa1).stableNorm();
494 if ( functor(wa2, wa4) < 0)
497 fnorm1 = wa4.stableNorm();
505 wa3 = R.template triangularView<Upper>()*wa1 + qtf;
506 temp = wa3.stableNorm();
514 ratio = actred / prered;
517 if (ratio < Scalar(.1)) {
520 delta = Scalar(.5) * delta;
524 if (ratio >= Scalar(.5) || ncsuc > 1)
525 delta = (
std::max)(delta, pnorm / Scalar(.5));
526 if (
abs(ratio - 1.) <= Scalar(.1)) {
527 delta = pnorm / Scalar(.5);
532 if (ratio >= Scalar(1
e-4)) {
535 wa2 = diag.cwiseProduct(
x);
537 xnorm = wa2.stableNorm();
544 if (actred >= Scalar(.001))
548 if (actred >= Scalar(.1))
552 if (delta <= parameters.xtol * xnorm || fnorm == 0.)
556 if (nfev >= parameters.maxfev)
571 wa1 = diag.cwiseProduct( diag.cwiseProduct(wa1)/pnorm );
572 wa2 = fjac.transpose() * wa4;
573 if (ratio >= Scalar(1
e-4))
575 wa2 = (wa2-wa3)/pnorm;
578 internal::r1updt<Scalar>(R, wa1, v_givens, w_givens, wa2, wa3, &sing);
579 internal::r1mpyq<Scalar>(
n,
n, fjac.data(), v_givens, w_givens);
580 internal::r1mpyq<Scalar>(1,
n, qtf.data(), v_givens, w_givens);
587 template<
typename FunctorType,
typename Scalar>
595 status = solveNumericalDiffOneStep(
x);
Array< double, 1, 3 > e(1./3., 0.5, 2.)
HouseholderSequenceType householderQ() const
const MatrixType & matrixQR() const
TransposeReturnType transpose() const
Finds a zero of a system of n nonlinear functions in n variables by a modification of the Powell hybr...
HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x)
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x)
void resetParameters(void)
HybridNonLinearSolverSpace::Status solveInit(FVectorType &x)
HybridNonLinearSolver & operator=(const HybridNonLinearSolver &)
HybridNonLinearSolverSpace::Status solve(FVectorType &x)
Matrix< Scalar, Dynamic, Dynamic > JacobianType
HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x)
Matrix< Scalar, Dynamic, 1 > FVectorType
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x)
HybridNonLinearSolverSpace::Status hybrj1(FVectorType &x, const Scalar tol=numext::sqrt(NumTraits< Scalar >::epsilon()))
HybridNonLinearSolver(FunctorType &_functor)
HybridNonLinearSolverSpace::Status hybrd1(FVectorType &x, const Scalar tol=numext::sqrt(NumTraits< Scalar >::epsilon()))
Matrix< Scalar, Dynamic, Dynamic > UpperTriangularType
@ ImproperInputParameters
@ NotMakingProgressIterations
@ TooManyFunctionEvaluation
@ NotMakingProgressJacobian
DenseIndex fdjac1(const FunctorType &Functor, Matrix< Scalar, Dynamic, 1 > &x, Matrix< Scalar, Dynamic, 1 > &fvec, Matrix< Scalar, Dynamic, Dynamic > &fjac, DenseIndex ml, DenseIndex mu, Scalar epsfcn)
EIGEN_ALWAYS_INLINE double sqrt(const double &x)
: TensorContractionSycl.h, provides various tensor contraction kernel for SYCL backend
Eigen::AutoDiffScalar< EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Eigen::internal::remove_all_t< DerType >, typename Eigen::internal::traits< Eigen::internal::remove_all_t< DerType >>::Scalar, product) > sqrt(const Eigen::AutoDiffScalar< DerType > &x)
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)
CleanedUpDerType< DerType >::type() max(const AutoDiffScalar< DerType > &x, const T &y)
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
adouble abs(const adouble &x)
Index nb_of_superdiagonals