Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/micpsolver/estimate_cond_number.hpp
c551311e303ereactmicp/unsaturated

Last Change
Raw File

estimate_cond_number.hpp
View Options

/*-------------------------------------------------------
- Module : micpsolver
- File : estimate_cond_number.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#ifndef SPECMICP_MICPSOLVER_ESTIMATECONDNUMBER_HPP
#define SPECMICP_MICPSOLVER_ESTIMATECONDNUMBER_HPP
#include <Eigen/Dense>
//! \file estimate_cond_number.hpp Estimate the condition number of a matrix
namespace specmicp {
//! \namespace micpsolver Namespace for the MiCP Solver
namespace micpsolver {
//! \brief Estimate the condition number of a dense square triangular matrix
//!
//! References :
//! - \cite Cline1979
//! - \cite Hager1984
//!
//! @param tmatrix a triangular view of a dense matrix
//! @param maxit maximum number of iterations done by the algorithm (2 triangular solve by iteration)
//!
template <class MatrixType, unsigned int mode>
double estimate_condition_number(Eigen::TriangularView<MatrixType, mode> tmatrix, int maxit=10)
{
const int n = tmatrix.cols();
Eigen::VectorXd x = 1/n*Eigen::VectorXd::Ones(n);
Eigen::VectorXd y(n);
int cnt = 0;
y = tmatrix.solve(x);
while (cnt < maxit)
{
for (int i=0; i<n; ++i)
{
if (y(i) >= 0) y(i) = 1;
else y(i) = -1;
}
tmatrix.solveInPlace(y);
y.reverseInPlace(); // transpose
int j;
const double maxi = y.array().abs().maxCoeff(&j);
const double gamma = y.dot(x);
if (maxi <= gamma) break; // This is a local maximum
x= Eigen::VectorXd::Zero(n);
x(j) = 1;
y = tmatrix.solve(x);
++cnt;
}
// norm tmatrix
// ------------
double nnorm;
if (mode == Eigen::Lower)
{
nnorm = std::abs(tmatrix(0, 0));
for (int i=1; i<n; ++i)
{
double normrow = 0;
for (int j=0; j<i+1; ++j)
{
normrow += std::abs(tmatrix(i, j));
}
nnorm = std::max(nnorm, normrow);
}
}
else if (mode == Eigen::Upper)
{
nnorm = std::abs(tmatrix(n-1, n-1));
for (int i=n-2; i>-1; --i)
{
double normrow = 0;
for (int j=i; j<n; ++j)
{
normrow += std::abs(tmatrix(i, j));
}
nnorm = std::max(nnorm, normrow);
}
}
// Return ||A||*||A^-1||
return nnorm*x.lpNorm<1>();
}
} // end namespace micpsolver
} // end namespace specmicp
#endif // SPECMICP_MICPSOLVER_ESTIMATECONDNUMBER_HPP
c4science · Help