10 #ifndef EIGEN_ITERSCALING_H
11 #define EIGEN_ITERSCALING_H
46 template<
typename _MatrixType>
50 typedef _MatrixType MatrixType;
51 typedef typename MatrixType::Scalar Scalar;
52 typedef typename MatrixType::Index Index;
55 IterScaling() { init(); }
57 IterScaling(
const MatrixType& matrix)
72 void compute (
const MatrixType& mat)
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)
104 Dr(i) = std::sqrt(Dr(i));
105 Dc(i) = std::sqrt(Dc(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();
133 }
while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
134 m_isInitialized =
true;
141 void computeRef (MatrixType& mat)
148 VectorXd& LeftScaling()
155 VectorXd& RightScaling()
162 void setTolerance(
double tol)
173 m_isInitialized =
false;
177 mutable ComputationInfo m_info;
178 bool m_isInitialized;