10 #ifndef EIGEN_ITERSCALING_H 11 #define EIGEN_ITERSCALING_H 46 template<
typename _MatrixType>
51 typedef typename MatrixType::Scalar
Scalar;
52 typedef typename MatrixType::Index
Index;
76 eigen_assert((m>0 && m == n) &&
"Please give a non - empty matrix");
82 VectorXd Dr, Dc, DrRes, DcRes;
83 Dr.resize(m); Dc.resize(n);
84 DrRes.resize(m); DcRes.resize(n);
85 double EpsRow = 1.0, EpsCol = 1.0;
90 Dr.setZero(); Dc.setZero();
91 for (
int k=0; k<
m_matrix.outerSize(); ++k)
93 for (
typename MatrixType::InnerIterator it(
m_matrix, k); it; ++it)
95 if ( Dr(it.row()) <
abs(it.value()) )
96 Dr(it.row()) =
abs(it.value());
98 if ( Dc(it.col()) <
abs(it.value()) )
99 Dc(it.col()) =
abs(it.value());
102 for (
int i = 0; i < m; ++i)
108 for (
int i = 0; i < m; ++i)
114 DrRes.setZero(); DcRes.setZero();
115 for (
int k=0; k<
m_matrix.outerSize(); ++k)
117 for (
typename MatrixType::InnerIterator it(
m_matrix, k); it; ++it)
119 it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
121 if ( DrRes(it.row()) <
abs(it.value()) )
122 DrRes(it.row()) =
abs(it.value());
124 if ( DcRes(it.col()) <
abs(it.value()) )
125 DcRes(it.col()) =
abs(it.value());
128 DrRes.array() = (1-DrRes.array()).
abs();
129 EpsRow = DrRes.maxCoeff();
130 DcRes.array() = (1-DcRes.array()).
abs();
131 EpsCol = DcRes.maxCoeff();
IntermediateState sqrt(const Expression &arg)
iterative scaling algorithm to equilibrate rows and column norms in matrices
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
MatrixType::Scalar Scalar
void compute(const MatrixType &mat)
VectorXd & RightScaling()
void computeRef(MatrixType &mat)
void setTolerance(double tol)
IterScaling(const MatrixType &matrix)