Page MenuHomec4science
Files F63765274

estimate_cond_number.hpp
No OneTemporary
Actions

File Metadata

Created
Wed, May 22, 08:34

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 {
//! \brief Estimate the condition number of a dense square triangular matrix
//!
//! References :
//! - A. K. Cline, C. B. Moler, G. W. Stewart, and J. H. Wilkinson.
//! An Estimate for the Condition Number of a Matrix.
//! SIAM Journal on Numerical Analysis, 16(2):368-375, April 1979.
//! - W. W. Hager. Condition estimates.
//! SIAM Journal on Scientific and statistical Computing,
//! 5(2):311-316, June 1984.
//!
//! @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;
while (cnt < maxit)
{
y = tmatrix.solve(x);
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;
++cnt;
}
y = tmatrix.solve(x);
// 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

Event Timeline

c4science · Help