diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-03-21 14:04:53 +0400 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2014-03-21 14:04:53 +0400 |
commit | 3411146984316c97f56543333a46f47aeb7f9d35 (patch) | |
tree | 5de608a3c18ff2a5459fd2191609dd882ad86213 /extern/Eigen3/Eigen/src/SVD/JacobiSVD.h | |
parent | 1781928f9d720fa1dc4811afb69935b35aa89878 (diff) |
Update Eigen to version 3.2.1
Main purpose of this is to have SparseLU solver which
we can use now as a replacement to opennl library.
Diffstat (limited to 'extern/Eigen3/Eigen/src/SVD/JacobiSVD.h')
-rw-r--r-- | extern/Eigen3/Eigen/src/SVD/JacobiSVD.h | 89 |
1 files changed, 52 insertions, 37 deletions
diff --git a/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h b/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h index a7dbf073766..f44995cd39c 100644 --- a/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h +++ b/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h @@ -78,7 +78,8 @@ public: { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { - m_qr = FullPivHouseholderQR<MatrixType>(svd.rows(), svd.cols()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.rows(), svd.cols()); } if (svd.m_computeFullU) m_workspace.resize(svd.rows()); } @@ -96,7 +97,8 @@ public: return false; } private: - FullPivHouseholderQR<MatrixType> m_qr; + typedef FullPivHouseholderQR<MatrixType> QRType; + QRType m_qr; WorkspaceType m_workspace; }; @@ -121,7 +123,8 @@ public: { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { - m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.cols(), svd.rows()); } m_adjoint.resize(svd.cols(), svd.rows()); if (svd.m_computeFullV) m_workspace.resize(svd.cols()); @@ -141,7 +144,8 @@ public: else return false; } private: - FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; + typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; + QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; typename internal::plain_row_type<MatrixType>::type m_workspace; }; @@ -158,7 +162,8 @@ public: { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { - m_qr = ColPivHouseholderQR<MatrixType>(svd.rows(), svd.cols()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.rows(), svd.cols()); } if (svd.m_computeFullU) m_workspace.resize(svd.rows()); else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); @@ -183,7 +188,8 @@ public: } private: - ColPivHouseholderQR<MatrixType> m_qr; + typedef ColPivHouseholderQR<MatrixType> QRType; + QRType m_qr; typename internal::plain_col_type<MatrixType>::type m_workspace; }; @@ -209,7 +215,8 @@ public: { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { - m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.cols(), svd.rows()); } if (svd.m_computeFullV) m_workspace.resize(svd.cols()); else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); @@ -237,7 +244,8 @@ public: } private: - ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; + typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; + QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; typename internal::plain_row_type<MatrixType>::type m_workspace; }; @@ -254,7 +262,8 @@ public: { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { - m_qr = HouseholderQR<MatrixType>(svd.rows(), svd.cols()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.rows(), svd.cols()); } if (svd.m_computeFullU) m_workspace.resize(svd.rows()); else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); @@ -278,7 +287,8 @@ public: return false; } private: - HouseholderQR<MatrixType> m_qr; + typedef HouseholderQR<MatrixType> QRType; + QRType m_qr; typename internal::plain_col_type<MatrixType>::type m_workspace; }; @@ -304,7 +314,8 @@ public: { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { - m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(svd.cols(), svd.rows()); + m_qr.~QRType(); + ::new (&m_qr) QRType(svd.cols(), svd.rows()); } if (svd.m_computeFullV) m_workspace.resize(svd.cols()); else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); @@ -332,7 +343,8 @@ public: } private: - HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; + typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; + QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; typename internal::plain_row_type<MatrixType>::type m_workspace; }; @@ -359,15 +371,19 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> typedef typename SVD::Index Index; static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) { + using std::sqrt; Scalar z; JacobiRotation<Scalar> rot; - RealScalar n = sqrt(abs2(work_matrix.coeff(p,p)) + abs2(work_matrix.coeff(q,p))); + RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); if(n==0) { z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); - z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + if(work_matrix.coeff(q,q)!=Scalar(0)) + z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + else + z = Scalar(0); work_matrix.row(q) *= z; if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); } @@ -398,9 +414,10 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation<RealScalar> *j_left, JacobiRotation<RealScalar> *j_right) { + using std::sqrt; Matrix<RealScalar,2,2> m; - m << real(matrix.coeff(p,p)), real(matrix.coeff(p,q)), - real(matrix.coeff(q,p)), real(matrix.coeff(q,q)); + m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), + numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); JacobiRotation<RealScalar> rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); @@ -412,7 +429,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, else { RealScalar u = d / t; - rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u)); + rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); rot1.s() = rot1.c() * u; } m.applyOnTheLeft(0,1,rot1); @@ -709,12 +726,14 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u } m_diagSize = (std::min)(m_rows, m_cols); m_singularValues.resize(m_diagSize); - m_matrixU.resize(m_rows, m_computeFullU ? m_rows - : m_computeThinU ? m_diagSize - : 0); - m_matrixV.resize(m_cols, m_computeFullV ? m_cols - : m_computeThinV ? m_diagSize - : 0); + if(RowsAtCompileTime==Dynamic) + m_matrixU.resize(m_rows, m_computeFullU ? m_rows + : m_computeThinU ? m_diagSize + : 0); + if(ColsAtCompileTime==Dynamic) + m_matrixV.resize(m_cols, m_computeFullV ? m_cols + : m_computeThinV ? m_diagSize + : 0); m_workMatrix.resize(m_diagSize, m_diagSize); if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); @@ -725,6 +744,7 @@ template<typename MatrixType, int QRPreconditioner> JacobiSVD<MatrixType, QRPreconditioner>& JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) { + using std::abs; allocate(matrix.rows(), matrix.cols(), computationOptions); // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, @@ -762,9 +782,9 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't // keep us iterating forever. Similarly, small denormal numbers are considered zero. using std::max; - RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)), - internal::abs(m_workMatrix.coeff(q,q)))); - if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold) + RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), + abs(m_workMatrix.coeff(q,q)))); + if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) { finished = false; @@ -788,7 +808,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig for(Index i = 0; i < m_diagSize; ++i) { - RealScalar a = internal::abs(m_workMatrix.coeff(i,i)); + RealScalar a = abs(m_workMatrix.coeff(i,i)); m_singularValues.coeffRef(i) = a; if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; } @@ -833,17 +853,12 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> // A = U S V^* // So A^{-1} = V S^{-1} U^* - Index diagSize = (std::min)(dec().rows(), dec().cols()); - typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize); - + Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp; Index nonzeroSingVals = dec().nonzeroSingularValues(); - invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse(); - invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); - - dst = dec().matrixV().leftCols(diagSize) - * invertedSingVals.asDiagonal() - * dec().matrixU().leftCols(diagSize).adjoint() - * rhs(); + + tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs(); + tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp; + dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp; } }; } // end namespace internal |