![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 00011 #define EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00026 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> 00027 EIGEN_DONT_INLINE 00028 void least_square_conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, 00029 const Preconditioner& precond, Index& iters, 00030 typename Dest::RealScalar& tol_error) 00031 { 00032 using std::sqrt; 00033 using std::abs; 00034 typedef typename Dest::RealScalar RealScalar; 00035 typedef typename Dest::Scalar Scalar; 00036 typedef Matrix<Scalar,Dynamic,1> VectorType; 00037 00038 RealScalar tol = tol_error; 00039 Index maxIters = iters; 00040 00041 Index m = mat.rows(), n = mat.cols(); 00042 00043 VectorType residual = rhs - mat * x; 00044 VectorType normal_residual = mat.adjoint() * residual; 00045 00046 RealScalar rhsNorm2 = (mat.adjoint()*rhs).squaredNorm(); 00047 if(rhsNorm2 == 0) 00048 { 00049 x.setZero(); 00050 iters = 0; 00051 tol_error = 0; 00052 return; 00053 } 00054 RealScalar threshold = tol*tol*rhsNorm2; 00055 RealScalar residualNorm2 = normal_residual.squaredNorm(); 00056 if (residualNorm2 < threshold) 00057 { 00058 iters = 0; 00059 tol_error = sqrt(residualNorm2 / rhsNorm2); 00060 return; 00061 } 00062 00063 VectorType p(n); 00064 p = precond.solve(normal_residual); // initial search direction 00065 00066 VectorType z(n), tmp(m); 00067 RealScalar absNew = numext::real(normal_residual.dot(p)); // the square of the absolute value of r scaled by invM 00068 Index i = 0; 00069 while(i < maxIters) 00070 { 00071 tmp.noalias() = mat * p; 00072 00073 Scalar alpha = absNew / tmp.squaredNorm(); // the amount we travel on dir 00074 x += alpha * p; // update solution 00075 residual -= alpha * tmp; // update residual 00076 normal_residual = mat.adjoint() * residual; // update residual of the normal equation 00077 00078 residualNorm2 = normal_residual.squaredNorm(); 00079 if(residualNorm2 < threshold) 00080 break; 00081 00082 z = precond.solve(normal_residual); // approximately solve for "A'A z = normal_residual" 00083 00084 RealScalar absOld = absNew; 00085 absNew = numext::real(normal_residual.dot(z)); // update the absolute value of r 00086 RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction 00087 p = z + beta * p; // update search direction 00088 i++; 00089 } 00090 tol_error = sqrt(residualNorm2 / rhsNorm2); 00091 iters = i; 00092 } 00093 00094 } 00095 00096 template< typename _MatrixType, 00097 typename _Preconditioner = LeastSquareDiagonalPreconditioner<typename _MatrixType::Scalar> > 00098 class LeastSquaresConjugateGradient; 00099 00100 namespace internal { 00101 00102 template< typename _MatrixType, typename _Preconditioner> 00103 struct traits<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> > 00104 { 00105 typedef _MatrixType MatrixType; 00106 typedef _Preconditioner Preconditioner; 00107 }; 00108 00109 } 00110 00148 template< typename _MatrixType, typename _Preconditioner> 00149 class LeastSquaresConjugateGradient : public IterativeSolverBase<LeastSquaresConjugateGradient<_MatrixType,_Preconditioner> > 00150 { 00151 typedef IterativeSolverBase<LeastSquaresConjugateGradient> Base; 00152 using Base::matrix; 00153 using Base::m_error; 00154 using Base::m_iterations; 00155 using Base::m_info; 00156 using Base::m_isInitialized; 00157 public: 00158 typedef _MatrixType MatrixType; 00159 typedef typename MatrixType::Scalar Scalar; 00160 typedef typename MatrixType::RealScalar RealScalar; 00161 typedef _Preconditioner Preconditioner; 00162 00163 public: 00164 00166 LeastSquaresConjugateGradient() : Base() {} 00167 00178 template<typename MatrixDerived> 00179 explicit LeastSquaresConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {} 00180 00181 ~LeastSquaresConjugateGradient() {} 00182 00184 template<typename Rhs,typename Dest> 00185 void _solve_with_guess_impl(const Rhs& b, Dest& x) const 00186 { 00187 m_iterations = Base::maxIterations(); 00188 m_error = Base::m_tolerance; 00189 00190 for(Index j=0; j<b.cols(); ++j) 00191 { 00192 m_iterations = Base::maxIterations(); 00193 m_error = Base::m_tolerance; 00194 00195 typename Dest::ColXpr xj(x,j); 00196 internal::least_square_conjugate_gradient(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); 00197 } 00198 00199 m_isInitialized = true; 00200 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; 00201 } 00202 00204 using Base::_solve_impl; 00205 template<typename Rhs,typename Dest> 00206 void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const 00207 { 00208 x.setZero(); 00209 _solve_with_guess_impl(b.derived(),x); 00210 } 00211 00212 }; 00213 00214 } // end namespace Eigen 00215 00216 #endif // EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H