25 #ifndef EIGEN_SCALING_H
26 #define EIGEN_SCALING_H
60 using namespace Eigen;
61 template<
typename _MatrixType>
65 typedef _MatrixType MatrixType;
66 typedef typename MatrixType::Scalar Scalar;
67 typedef typename MatrixType::Index Index;
72 Scaling(
const MatrixType& matrix)
87 void compute (
const MatrixType& mat)
91 assert((m>0 && m == n) &&
"Please give a non - empty matrix");
97 VectorXd Dr, Dc, DrRes, DcRes;
98 Dr.resize(m); Dc.resize(n);
99 DrRes.resize(m); DcRes.resize(n);
100 double EpsRow = 1.0, EpsCol = 1.0;
105 Dr.setZero(); Dc.setZero();
106 for (
int k=0; k<m_matrix.outerSize(); ++k)
108 for (
typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
110 if ( Dr(it.row()) <
abs(it.value()) )
111 Dr(it.row()) =
abs(it.value());
113 if ( Dc(it.col()) <
abs(it.value()) )
114 Dc(it.col()) =
abs(it.value());
117 for (
int i = 0; i < m; ++i)
119 Dr(i) = std::sqrt(Dr(i));
120 Dc(i) = std::sqrt(Dc(i));
123 for (
int i = 0; i < m; ++i)
129 DrRes.setZero(); DcRes.setZero();
130 for (
int k=0; k<m_matrix.outerSize(); ++k)
132 for (
typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
134 it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
136 if ( DrRes(it.row()) <
abs(it.value()) )
137 DrRes(it.row()) =
abs(it.value());
139 if ( DcRes(it.col()) <
abs(it.value()) )
140 DcRes(it.col()) =
abs(it.value());
143 DrRes.array() = (1-DrRes.array()).
abs();
144 EpsRow = DrRes.maxCoeff();
145 DcRes.array() = (1-DcRes.array()).
abs();
146 EpsCol = DcRes.maxCoeff();
148 }
while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
149 m_isInitialized =
true;
156 void computeRef (MatrixType& mat)
163 VectorXd& LeftScaling()
170 VectorXd& RightScaling()
177 void setTolerance(
double tol)
188 m_isInitialized =
false;
193 bool m_isInitialized;