diff options
Diffstat (limited to 'Eigen/src/SVD')
-rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 51 | ||||
-rw-r--r-- | Eigen/src/SVD/SVD.h | 59 |
2 files changed, 76 insertions, 34 deletions
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index b7cab023b..9323c0180 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -93,10 +93,35 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD public: + /** \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via JacobiSVD::compute(const MatrixType&). + */ JacobiSVD() : m_isInitialized(false) {} - JacobiSVD(const MatrixType& matrix) : m_isInitialized(false) + + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa JacobiSVD() + */ + JacobiSVD(int rows, int cols) : m_matrixU(rows, rows), + m_matrixV(cols, cols), + m_singularValues(std::min(rows, cols)), + m_workMatrix(rows, cols), + m_isInitialized(false) {} + + JacobiSVD(const MatrixType& matrix) : m_matrixU(matrix.rows(), matrix.rows()), + m_matrixV(matrix.cols(), matrix.cols()), + m_singularValues(), + m_workMatrix(), + m_isInitialized(false) { + const int minSize = std::min(matrix.rows(), matrix.cols()); + m_singularValues.resize(minSize); + m_workMatrix.resize(minSize, minSize); compute(matrix); } @@ -124,6 +149,7 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD MatrixUType m_matrixU; MatrixVType m_matrixV; SingularValuesType m_singularValues; + WorkMatrixType m_workMatrix; bool m_isInitialized; template<typename _MatrixType, unsigned int _Options, bool _IsComplex> @@ -289,17 +315,16 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true> template<typename MatrixType, unsigned int Options> JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix) { - WorkMatrixType work_matrix; int rows = matrix.rows(); int cols = matrix.cols(); int diagSize = std::min(rows, cols); m_singularValues.resize(diagSize); const RealScalar precision = 2 * NumTraits<Scalar>::epsilon(); - if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this) - && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this)) + if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, m_workMatrix, *this) + && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, m_workMatrix, *this)) { - work_matrix = matrix.block(0,0,diagSize,diagSize); + m_workMatrix = matrix.block(0,0,diagSize,diagSize); if(ComputeU) m_matrixU.setIdentity(rows,rows); if(ComputeV) m_matrixV.setIdentity(cols,cols); } @@ -312,19 +337,19 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma { for(int q = 0; q < p; ++q) { - if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) + if(std::max(ei_abs(m_workMatrix.coeff(p,q)),ei_abs(m_workMatrix.coeff(q,p))) + > std::max(ei_abs(m_workMatrix.coeff(p,p)),ei_abs(m_workMatrix.coeff(q,q)))*precision) { finished = false; - ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q); + ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, q); PlanarRotation<RealScalar> j_left, j_right; - ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); + ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); - work_matrix.applyOnTheLeft(p,q,j_left); + m_workMatrix.applyOnTheLeft(p,q,j_left); if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); - work_matrix.applyOnTheRight(p,q,j_right); + m_workMatrix.applyOnTheRight(p,q,j_right); if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); } } @@ -333,9 +358,9 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma for(int i = 0; i < diagSize; ++i) { - RealScalar a = ei_abs(work_matrix.coeff(i,i)); + RealScalar a = ei_abs(m_workMatrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; - if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; + if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; } for(int i = 0; i < diagSize; i++) diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index ed0e69f91..a9e22dfd4 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -73,11 +73,25 @@ template<typename _MatrixType> class SVD */ SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {} - SVD(const MatrixType& matrix) - : m_matU(matrix.rows(), matrix.rows()), - m_matV(matrix.cols(),matrix.cols()), - m_sigma(matrix.cols()), - m_isInitialized(false) + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem \a size. + * \sa JacobiSVD() + */ + SVD(int rows, int cols) : m_matU(rows, rows), + m_matV(cols,cols), + m_sigma(std::min(rows, cols)), + m_workMatrix(rows, cols), + m_rv1(cols), + m_isInitialized(false) {} + + SVD(const MatrixType& matrix) : m_matU(matrix.rows(), matrix.rows()), + m_matV(matrix.cols(),matrix.cols()), + m_sigma(std::min(matrix.rows(), matrix.cols())), + m_workMatrix(matrix.rows(), matrix.cols()), + m_rv1(matrix.cols()), + m_isInitialized(false) { compute(matrix); } @@ -165,6 +179,8 @@ template<typename _MatrixType> class SVD MatrixVType m_matV; /** \internal */ SingularValuesType m_sigma; + MatrixType m_workMatrix; + RowVector m_rv1; bool m_isInitialized; int m_rows, m_cols; }; @@ -185,11 +201,12 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) m_matU.setZero(); m_sigma.resize(n); m_matV.resize(n,n); + m_workMatrix = matrix; int max_iters = 30; MatrixVType& V = m_matV; - MatrixType A = matrix; + MatrixType& A = m_workMatrix; SingularValuesType& W = m_sigma; bool flag; @@ -198,14 +215,14 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) bool convergence = true; Scalar eps = NumTraits<Scalar>::dummy_precision(); - RowVector rv1(n); + m_rv1.resize(n); g = scale = anorm = 0; // Householder reduction to bidiagonal form. for (i=0; i<n; i++) { l = i+2; - rv1[i] = scale*g; + m_rv1[i] = scale*g; g = s = scale = 0.0; if (i < m) { @@ -246,16 +263,16 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) g = -sign(ei_sqrt(s),f); h = f*g - s; A(i,l-1) = f-g; - rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h; + m_rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h; for (j=l-1; j<m; j++) { s = A.row(i).tail(n-l+1).dot(A.row(j).tail(n-l+1)); - A.row(j).tail(n-l+1) += s*rv1.tail(n-l+1).transpose(); + A.row(j).tail(n-l+1) += s*m_rv1.tail(n-l+1).transpose(); } A.row(i).tail(n-l+1) *= scale; } } - anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) ); + anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(m_rv1[i])) ); } // Accumulation of right-hand transformations. for (i=n-1; i>=0; i--) @@ -277,7 +294,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) V.col(i).tail(n-l).setZero(); } V(i, i) = 1.0; - g = rv1[i]; + g = m_rv1[i]; l = i; } // Accumulation of left-hand transformations. @@ -318,7 +335,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) nm = l-1; // Note that rv1[1] is always zero. //if ((double)(ei_abs(rv1[l])+anorm) == anorm) - if (l==0 || ei_abs(rv1[l]) <= eps*anorm) + if (l==0 || ei_abs(m_rv1[l]) <= eps*anorm) { flag = false; break; @@ -333,8 +350,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) s = 1.0; for (i=l ;i<k+1; i++) { - f = s*rv1[i]; - rv1[i] = c*rv1[i]; + f = s*m_rv1[i]; + m_rv1[i] = c*m_rv1[i]; //if ((double)(ei_abs(f)+anorm) == anorm) if (ei_abs(f) <= eps*anorm) break; @@ -363,8 +380,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) x = W[l]; // Shift from bottom 2-by-2 minor. nm = k-1; y = W[nm]; - g = rv1[nm]; - h = rv1[k]; + g = m_rv1[nm]; + h = m_rv1[k]; f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y); g = pythag(f,1.0); f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x; @@ -373,13 +390,13 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) for (j=l; j<=nm; j++) { i = j+1; - g = rv1[i]; + g = m_rv1[i]; y = W[i]; h = s*g; g = c*g; z = pythag(f,h); - rv1[j] = z; + m_rv1[j] = z; c = f/z; s = h/z; f = x*c + g*s; @@ -401,8 +418,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) x = c*y - s*g; A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s)); } - rv1[l] = 0.0; - rv1[k] = f; + m_rv1[l] = 0.0; + m_rv1[k] = f; W[k] = x; } } |