This commit is contained in:
Mark Borgerding
2010-02-22 21:44:30 -05:00
66 changed files with 676 additions and 557 deletions

View File

@@ -563,6 +563,8 @@ template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> >
AddCost = 1,
MulCost = 1
};
inline static Real epsilon() { return std::numeric_limits<Real>::epsilon(); }
inline static Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
};
}

View File

@@ -313,7 +313,7 @@ template<typename Derived> struct MatrixExponentialReturnValue
inline void evalTo(ResultType& result) const
{
const typename ei_eval<Derived>::type srcEvaluated = m_src.eval();
MatrixExponential<typename Derived::PlainMatrixType> me(srcEvaluated);
MatrixExponential<typename Derived::PlainObject> me(srcEvaluated);
me.compute(result);
}
@@ -327,7 +327,7 @@ template<typename Derived> struct MatrixExponentialReturnValue
template<typename Derived>
struct ei_traits<MatrixExponentialReturnValue<Derived> >
{
typedef typename Derived::PlainMatrixType ReturnMatrixType;
typedef typename Derived::PlainObject ReturnType;
};
/** \ingroup MatrixFunctions_Module

View File

@@ -178,9 +178,9 @@ class MatrixFunction<MatrixType, 1>
*
* This is morally a \c static \c const \c Scalar, but only
* integers can be static constant class members in C++. The
* separation constant is set to 0.01, a value taken from the
* separation constant is set to 0.1, a value taken from the
* paper by Davies and Higham. */
static const RealScalar separation() { return static_cast<RealScalar>(0.01); }
static const RealScalar separation() { return static_cast<RealScalar>(0.1); }
};
/** \brief Constructor.
@@ -492,14 +492,12 @@ typename MatrixFunction<MatrixType,1>::DynMatrixType MatrixFunction<MatrixType,1
template<typename Derived> class MatrixFunctionReturnValue
: public ReturnByValue<MatrixFunctionReturnValue<Derived> >
{
private:
public:
typedef typename ei_traits<Derived>::Scalar Scalar;
typedef typename ei_stem_function<Scalar>::type StemFunction;
public:
/** \brief Constructor.
/** \brief Constructor.
*
* \param[in] A %Matrix (expression) forming the argument of the
* matrix function.
@@ -516,7 +514,7 @@ template<typename Derived> class MatrixFunctionReturnValue
inline void evalTo(ResultType& result) const
{
const typename ei_eval<Derived>::type Aevaluated = m_A.eval();
MatrixFunction<typename Derived::PlainMatrixType> mf(Aevaluated, m_f);
MatrixFunction<typename Derived::PlainObject> mf(Aevaluated, m_f);
mf.compute(result);
}
@@ -531,7 +529,7 @@ template<typename Derived> class MatrixFunctionReturnValue
template<typename Derived>
struct ei_traits<MatrixFunctionReturnValue<Derived> >
{
typedef typename Derived::PlainMatrixType ReturnMatrixType;
typedef typename Derived::PlainObject ReturnType;
};

View File

@@ -57,7 +57,7 @@ class HybridNonLinearSolver
{
public:
HybridNonLinearSolver(FunctorType &_functor)
: functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; }
: functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; useExternalScaling=false;}
struct Parameters {
Parameters()
@@ -84,36 +84,18 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
HybridNonLinearSolverSpace::Status solveInit(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveOneStep(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solve(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveInit(FVectorType &x);
HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x);
HybridNonLinearSolverSpace::Status solve(FVectorType &x);
HybridNonLinearSolverSpace::Status hybrd1(
FVectorType &x,
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiff(
FVectorType &x,
const int mode=1
);
HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x);
HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x);
HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x);
void resetParameters(void) { parameters = Parameters(); }
Parameters parameters;
@@ -124,6 +106,7 @@ public:
int njev;
int iter;
Scalar fnorm;
bool useExternalScaling;
private:
FunctorType &functor;
int n;
@@ -160,18 +143,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrj1(
parameters.maxfev = 100*(n+1);
parameters.xtol = tol;
diag.setConstant(n, 1.);
return solve(
x,
2
);
useExternalScaling = true;
return solve(x);
}
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveInit(FVectorType &x)
{
n = x.size();
@@ -179,9 +157,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
fvec.resize(n);
qtf.resize(n);
fjac.resize(n, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
/* Function Body */
nfev = 0;
@@ -190,7 +168,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
/* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. )
return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return HybridNonLinearSolverSpace::ImproperInputParameters;
@@ -214,10 +192,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType &x)
{
int j;
std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -231,10 +206,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
@@ -260,7 +235,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
qtf = fjac.transpose() * fvec;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
while (true) {
@@ -372,14 +347,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solve(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solve(FVectorType &x)
{
HybridNonLinearSolverSpace::Status status = solveInit(x, mode);
HybridNonLinearSolverSpace::Status status = solveInit(x);
while (status==HybridNonLinearSolverSpace::Running)
status = solveOneStep(x, mode);
status = solveOneStep(x);
return status;
}
@@ -403,18 +375,13 @@ HybridNonLinearSolver<FunctorType,Scalar>::hybrd1(
parameters.xtol = tol;
diag.setConstant(n, 1.);
return solveNumericalDiff(
x,
2
);
useExternalScaling = true;
return solveNumericalDiff(x);
}
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(FVectorType &x)
{
n = x.size();
@@ -425,10 +392,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
qtf.resize(n);
fjac.resize(n, n);
fvec.resize(n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
/* Function Body */
nfev = 0;
@@ -437,7 +403,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
/* check the input parameters for errors. */
if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. )
return HybridNonLinearSolverSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return HybridNonLinearSolverSpace::ImproperInputParameters;
@@ -461,10 +427,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType &x)
{
int j;
std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -480,10 +443,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
wa2 = fjac.colwise().blueNorm();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.) ? 1. : wa2[j];
@@ -509,7 +472,7 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
qtf = fjac.transpose() * fvec;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
while (true) {
@@ -621,14 +584,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
template<typename FunctorType, typename Scalar>
HybridNonLinearSolverSpace::Status
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(
FVectorType &x,
const int mode
)
HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(FVectorType &x)
{
HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x, mode);
HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x);
while (status==HybridNonLinearSolverSpace::Running)
status = solveNumericalDiffOneStep(x, mode);
status = solveNumericalDiffOneStep(x);
return status;
}

View File

@@ -61,7 +61,7 @@ class LevenbergMarquardt
{
public:
LevenbergMarquardt(FunctorType &_functor)
: functor(_functor) { nfev = njev = iter = 0; fnorm=gnorm = 0.; }
: functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; }
struct Parameters {
Parameters()
@@ -87,18 +87,9 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
LevenbergMarquardtSpace::Status minimize(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeInit(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOneStep(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimize(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x);
static LevenbergMarquardtSpace::Status lmdif1(
FunctorType &functor,
@@ -112,18 +103,9 @@ public:
const Scalar tol = ei_sqrt(NumTraits<Scalar>::epsilon())
);
LevenbergMarquardtSpace::Status minimizeOptimumStorage(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(
FVectorType &x,
const int mode=1
);
LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType &x);
LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType &x);
void resetParameters(void) { parameters = Parameters(); }
@@ -135,6 +117,7 @@ public:
int njev;
int iter;
Scalar fnorm, gnorm;
bool useExternalScaling;
Scalar lm_param(void) { return par; }
private:
@@ -175,24 +158,18 @@ LevenbergMarquardt<FunctorType,Scalar>::lmder1(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimize(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimize(FVectorType &x)
{
LevenbergMarquardtSpace::Status status = minimizeInit(x, mode);
LevenbergMarquardtSpace::Status status = minimizeInit(x);
do {
status = minimizeOneStep(x, mode);
status = minimizeOneStep(x);
} while (status==LevenbergMarquardtSpace::Running);
return status;
}
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType &x)
{
n = x.size();
m = functor.values();
@@ -201,9 +178,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
wa4.resize(m);
fvec.resize(m);
fjac.resize(m, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
qtf.resize(n);
/* Function Body */
@@ -214,7 +191,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
@@ -235,10 +212,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType &x)
{
int j;
@@ -257,10 +231,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
fjac = qrfac.matrixQR();
permutation = qrfac.colsPermutation();
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.)? 1. : wa2[j];
@@ -290,7 +264,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
do {
@@ -406,10 +380,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType &x)
{
n = x.size();
m = functor.values();
@@ -423,9 +394,9 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
// The purpose it to only use a nxn matrix, instead of mxn here, so
// that we can handle cases where m>>n :
fjac.resize(n, n);
if (mode != 2)
if (!useExternalScaling)
diag.resize(n);
assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'");
assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'");
qtf.resize(n);
/* Function Body */
@@ -436,7 +407,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
if (mode == 2)
if (useExternalScaling)
for (int j = 0; j < n; ++j)
if (diag[j] <= 0.)
return LevenbergMarquardtSpace::ImproperInputParameters;
@@ -458,10 +429,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType &x)
{
int i, j;
bool sing;
@@ -514,10 +482,10 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
}
}
/* on the first iteration and if mode is 1, scale according */
/* on the first iteration and if external scaling is not used, scale according */
/* to the norms of the columns of the initial jacobian. */
if (iter == 1) {
if (mode != 2)
if (!useExternalScaling)
for (j = 0; j < n; ++j)
diag[j] = (wa2[j]==0.)? 1. : wa2[j];
@@ -541,7 +509,7 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
return LevenbergMarquardtSpace::CosinusTooSmall;
/* rescale if necessary. */
if (mode != 2)
if (!useExternalScaling)
diag = diag.cwiseMax(wa2);
do {
@@ -635,14 +603,11 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
template<typename FunctorType, typename Scalar>
LevenbergMarquardtSpace::Status
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(
FVectorType &x,
const int mode
)
LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(FVectorType &x)
{
LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x, mode);
LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x);
do {
status = minimizeOptimumStorageOneStep(x, mode);
status = minimizeOptimumStorageOneStep(x);
} while (status==LevenbergMarquardtSpace::Running);
return status;
}

View File

@@ -36,7 +36,7 @@
* \param Derived
*
*/
template<typename Derived> class SkylineMatrixBase : public AnyMatrixBase<Derived> {
template<typename Derived> class SkylineMatrixBase : public EigenBase<Derived> {
public:
typedef typename ei_traits<Derived>::Scalar Scalar;