Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F63767024
estimate_cond_number.hpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, May 22, 08:44
Size
2 KB
Mime Type
text/x-c++
Expires
Fri, May 24, 08:44 (2 d)
Engine
blob
Format
Raw Data
Handle
17816545
Attached To
rSPECMICP SpecMiCP / ReactMiCP
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
.
transpose
().
solveInPlace
(
y
);
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
);
// 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
Log In to Comment