00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_SCALING_H
00026 #define EIGEN_SCALING_H
00027
00059 using std::abs;
00060 using namespace Eigen;
00061 template<typename _MatrixType>
00062 class Scaling
00063 {
00064 public:
00065 typedef _MatrixType MatrixType;
00066 typedef typename MatrixType::Scalar Scalar;
00067 typedef typename MatrixType::Index Index;
00068
00069 public:
00070 Scaling() { init(); }
00071
00072 Scaling(const MatrixType& matrix)
00073 {
00074 init();
00075 compute(matrix);
00076 }
00077
00078 ~Scaling() { }
00079
00087 void compute (const MatrixType& mat)
00088 {
00089 int m = mat.rows();
00090 int n = mat.cols();
00091 assert((m>0 && m == n) && "Please give a non - empty matrix");
00092 m_left.resize(m);
00093 m_right.resize(n);
00094 m_left.setOnes();
00095 m_right.setOnes();
00096 m_matrix = mat;
00097 VectorXd Dr, Dc, DrRes, DcRes;
00098 Dr.resize(m); Dc.resize(n);
00099 DrRes.resize(m); DcRes.resize(n);
00100 double EpsRow = 1.0, EpsCol = 1.0;
00101 int its = 0;
00102 do
00103 {
00104
00105 Dr.setZero(); Dc.setZero();
00106 for (int k=0; k<m_matrix.outerSize(); ++k)
00107 {
00108 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00109 {
00110 if ( Dr(it.row()) < abs(it.value()) )
00111 Dr(it.row()) = abs(it.value());
00112
00113 if ( Dc(it.col()) < abs(it.value()) )
00114 Dc(it.col()) = abs(it.value());
00115 }
00116 }
00117 for (int i = 0; i < m; ++i)
00118 {
00119 Dr(i) = std::sqrt(Dr(i));
00120 Dc(i) = std::sqrt(Dc(i));
00121 }
00122
00123 for (int i = 0; i < m; ++i)
00124 {
00125 m_left(i) /= Dr(i);
00126 m_right(i) /= Dc(i);
00127 }
00128
00129 DrRes.setZero(); DcRes.setZero();
00130 for (int k=0; k<m_matrix.outerSize(); ++k)
00131 {
00132 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00133 {
00134 it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
00135
00136 if ( DrRes(it.row()) < abs(it.value()) )
00137 DrRes(it.row()) = abs(it.value());
00138
00139 if ( DcRes(it.col()) < abs(it.value()) )
00140 DcRes(it.col()) = abs(it.value());
00141 }
00142 }
00143 DrRes.array() = (1-DrRes.array()).abs();
00144 EpsRow = DrRes.maxCoeff();
00145 DcRes.array() = (1-DcRes.array()).abs();
00146 EpsCol = DcRes.maxCoeff();
00147 its++;
00148 }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
00149 m_isInitialized = true;
00150 }
00156 void computeRef (MatrixType& mat)
00157 {
00158 compute (mat);
00159 mat = m_matrix;
00160 }
00163 VectorXd& LeftScaling()
00164 {
00165 return m_left;
00166 }
00167
00170 VectorXd& RightScaling()
00171 {
00172 return m_right;
00173 }
00174
00177 void setTolerance(double tol)
00178 {
00179 m_tol = tol;
00180 }
00181
00182 protected:
00183
00184 void init()
00185 {
00186 m_tol = 1e-10;
00187 m_maxits = 5;
00188 m_isInitialized = false;
00189 }
00190
00191 MatrixType m_matrix;
00192 mutable ComputationInfo m_info;
00193 bool m_isInitialized;
00194 VectorXd m_left;
00195 VectorXd m_right;
00196 double m_tol;
00197 int m_maxits;
00198 };
00199
00200 #endif