1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2020 Chris Schoutrop <c.e.m.schoutrop@tue.nl> 5 // Copyright (C) 2020 Jens Wehner <j.wehner@esciencecenter.nl> 6 // Copyright (C) 2020 Jan van Dijk <j.v.dijk@tue.nl> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 13 #ifndef EIGEN_IDRS_H 14 #define EIGEN_IDRS_H 15 16 namespace Eigen 17 { 18 19 namespace internal 20 { 21 /** \internal Low-level Induced Dimension Reduction algoritm 22 \param A The matrix A 23 \param b The right hand side vector b 24 \param x On input and initial solution, on output the computed solution. 25 \param precond A preconditioner being able to efficiently solve for an 26 approximation of Ax=b (regardless of b) 27 \param iter On input the max number of iteration, on output the number of performed iterations. 28 \param relres On input the tolerance error, on output an estimation of the relative error. 29 \param S On input Number of the dimension of the shadow space. 30 \param smoothing switches residual smoothing on. 31 \param angle small omega lead to faster convergence at the expense of numerical stability 32 \param replacement switches on a residual replacement strategy to increase accuracy of residual at the expense of more Mat*vec products 33 \return false in the case of numerical issue, for example a break down of IDRS. 34 */ 35 template<typename Vector, typename RealScalar> omega(const Vector & t,const Vector & s,RealScalar angle)36 typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle) 37 { 38 using numext::abs; 39 typedef typename Vector::Scalar Scalar; 40 const RealScalar ns = s.norm(); 41 const RealScalar nt = t.norm(); 42 const Scalar ts = t.dot(s); 43 const RealScalar rho = abs(ts / (nt * ns)); 44 45 if (rho < angle) { 46 if (ts == Scalar(0)) { 47 return Scalar(0); 48 } 49 // Original relation for om is given by 50 // om = om * angle / rho; 51 // To alleviate potential (near) division by zero this can be rewritten as 52 // om = angle * (ns / nt) * (ts / abs(ts)) = angle * (ns / nt) * sgn(ts) 53 return angle * (ns / nt) * (ts / abs(ts)); 54 } 55 return ts / (nt * nt); 56 } 57 58 template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner> idrs(const MatrixType & A,const Rhs & b,Dest & x,const Preconditioner & precond,Index & iter,typename Dest::RealScalar & relres,Index S,bool smoothing,typename Dest::RealScalar angle,bool replacement)59 bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond, 60 Index& iter, 61 typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle, bool replacement) 62 { 63 typedef typename Dest::RealScalar RealScalar; 64 typedef typename Dest::Scalar Scalar; 65 typedef Matrix<Scalar, Dynamic, 1> VectorType; 66 typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType; 67 const Index N = b.size(); 68 S = S < x.rows() ? S : x.rows(); 69 const RealScalar tol = relres; 70 const Index maxit = iter; 71 72 Index replacements = 0; 73 bool trueres = false; 74 75 FullPivLU<DenseMatrixType> lu_solver; 76 77 DenseMatrixType P; 78 { 79 HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S)); 80 P = (qr.householderQ() * DenseMatrixType::Identity(N, S)); 81 } 82 83 const RealScalar normb = b.norm(); 84 85 if (internal::isApprox(normb, RealScalar(0))) 86 { 87 //Solution is the zero vector 88 x.setZero(); 89 iter = 0; 90 relres = 0; 91 return true; 92 } 93 // from http://homepage.tudelft.nl/1w5b5/IDRS/manual.pdf 94 // A peak in the residual is considered dangerously high if‖ri‖/‖b‖> C(tol/epsilon). 95 // With epsilon the 96 // relative machine precision. The factor tol/epsilon corresponds to the size of a 97 // finite precision number that is so large that the absolute round-off error in 98 // this number, when propagated through the process, makes it impossible to 99 // achieve the required accuracy.The factor C accounts for the accumulation of 100 // round-off errors. This parameter has beenset to 10−3. 101 // mp is epsilon/C 102 // 10^3 * eps is very conservative, so normally no residual replacements will take place. 103 // It only happens if things go very wrong. Too many restarts may ruin the convergence. 104 const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon(); 105 106 107 108 //Compute initial residual 109 const RealScalar tolb = tol * normb; //Relative tolerance 110 VectorType r = b - A * x; 111 112 VectorType x_s, r_s; 113 114 if (smoothing) 115 { 116 x_s = x; 117 r_s = r; 118 } 119 120 RealScalar normr = r.norm(); 121 122 if (normr <= tolb) 123 { 124 //Initial guess is a good enough solution 125 iter = 0; 126 relres = normr / normb; 127 return true; 128 } 129 130 DenseMatrixType G = DenseMatrixType::Zero(N, S); 131 DenseMatrixType U = DenseMatrixType::Zero(N, S); 132 DenseMatrixType M = DenseMatrixType::Identity(S, S); 133 VectorType t(N), v(N); 134 Scalar om = 1.; 135 136 //Main iteration loop, guild G-spaces: 137 iter = 0; 138 139 while (normr > tolb && iter < maxit) 140 { 141 //New right hand size for small system: 142 VectorType f = (r.adjoint() * P).adjoint(); 143 144 for (Index k = 0; k < S; ++k) 145 { 146 //Solve small system and make v orthogonal to P: 147 //c = M(k:s,k:s)\f(k:s); 148 lu_solver.compute(M.block(k , k , S -k, S - k )); 149 VectorType c = lu_solver.solve(f.segment(k , S - k )); 150 //v = r - G(:,k:s)*c; 151 v = r - G.rightCols(S - k ) * c; 152 //Preconditioning 153 v = precond.solve(v); 154 155 //Compute new U(:,k) and G(:,k), G(:,k) is in space G_j 156 U.col(k) = U.rightCols(S - k ) * c + om * v; 157 G.col(k) = A * U.col(k ); 158 159 //Bi-Orthogonalise the new basis vectors: 160 for (Index i = 0; i < k-1 ; ++i) 161 { 162 //alpha = ( P(:,i)'*G(:,k) )/M(i,i); 163 Scalar alpha = P.col(i ).dot(G.col(k )) / M(i, i ); 164 G.col(k ) = G.col(k ) - alpha * G.col(i ); 165 U.col(k ) = U.col(k ) - alpha * U.col(i ); 166 } 167 168 //New column of M = P'*G (first k-1 entries are zero) 169 //M(k:s,k) = (G(:,k)'*P(:,k:s))'; 170 M.block(k , k , S - k , 1) = (G.col(k ).adjoint() * P.rightCols(S - k )).adjoint(); 171 172 if (internal::isApprox(M(k,k), Scalar(0))) 173 { 174 return false; 175 } 176 177 //Make r orthogonal to q_i, i = 0..k-1 178 Scalar beta = f(k ) / M(k , k ); 179 r = r - beta * G.col(k ); 180 x = x + beta * U.col(k ); 181 normr = r.norm(); 182 183 if (replacement && normr > tolb / mp) 184 { 185 trueres = true; 186 } 187 188 //Smoothing: 189 if (smoothing) 190 { 191 t = r_s - r; 192 //gamma is a Scalar, but the conversion is not allowed 193 Scalar gamma = t.dot(r_s) / t.norm(); 194 r_s = r_s - gamma * t; 195 x_s = x_s - gamma * (x_s - x); 196 normr = r_s.norm(); 197 } 198 199 if (normr < tolb || iter == maxit) 200 { 201 break; 202 } 203 204 //New f = P'*r (first k components are zero) 205 if (k < S-1) 206 { 207 f.segment(k + 1, S - (k + 1) ) = f.segment(k + 1 , S - (k + 1)) - beta * M.block(k + 1 , k , S - (k + 1), 1); 208 } 209 }//end for 210 211 if (normr < tolb || iter == maxit) 212 { 213 break; 214 } 215 216 //Now we have sufficient vectors in G_j to compute residual in G_j+1 217 //Note: r is already perpendicular to P so v = r 218 //Preconditioning 219 v = r; 220 v = precond.solve(v); 221 222 //Matrix-vector multiplication: 223 t = A * v; 224 225 //Computation of a new omega 226 om = internal::omega(t, r, angle); 227 228 if (om == RealScalar(0.0)) 229 { 230 return false; 231 } 232 233 r = r - om * t; 234 x = x + om * v; 235 normr = r.norm(); 236 237 if (replacement && normr > tolb / mp) 238 { 239 trueres = true; 240 } 241 242 //Residual replacement? 243 if (trueres && normr < normb) 244 { 245 r = b - A * x; 246 trueres = false; 247 replacements++; 248 } 249 250 //Smoothing: 251 if (smoothing) 252 { 253 t = r_s - r; 254 Scalar gamma = t.dot(r_s) /t.norm(); 255 r_s = r_s - gamma * t; 256 x_s = x_s - gamma * (x_s - x); 257 normr = r_s.norm(); 258 } 259 260 iter++; 261 262 }//end while 263 264 if (smoothing) 265 { 266 x = x_s; 267 } 268 relres=normr/normb; 269 return true; 270 } 271 272 } // namespace internal 273 274 template <typename _MatrixType, typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> > 275 class IDRS; 276 277 namespace internal 278 { 279 280 template <typename _MatrixType, typename _Preconditioner> 281 struct traits<Eigen::IDRS<_MatrixType, _Preconditioner> > 282 { 283 typedef _MatrixType MatrixType; 284 typedef _Preconditioner Preconditioner; 285 }; 286 287 } // namespace internal 288 289 290 /** \ingroup IterativeLinearSolvers_Module 291 * \brief The Induced Dimension Reduction method (IDR(s)) is a short-recurrences Krylov method for sparse square problems. 292 * 293 * This class allows to solve for A.x = b sparse linear problems. The vectors x and b can be either dense or sparse. 294 * he Induced Dimension Reduction method, IDR(), is a robust and efficient short-recurrence Krylov subspace method for 295 * solving large nonsymmetric systems of linear equations. 296 * 297 * For indefinite systems IDR(S) outperforms both BiCGStab and BiCGStab(L). Additionally, IDR(S) can handle matrices 298 * with complex eigenvalues more efficiently than BiCGStab. 299 * 300 * Many problems that do not converge for BiCGSTAB converge for IDR(s) (for larger values of s). And if both methods 301 * converge the convergence for IDR(s) is typically much faster for difficult systems (for example indefinite problems). 302 * 303 * IDR(s) is a limited memory finite termination method. In exact arithmetic it converges in at most N+N/s iterations, 304 * with N the system size. It uses a fixed number of 4+3s vector. In comparison, BiCGSTAB terminates in 2N iterations 305 * and uses 7 vectors. GMRES terminates in at most N iterations, and uses I+3 vectors, with I the number of iterations. 306 * Restarting GMRES limits the memory consumption, but destroys the finite termination property. 307 * 308 * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix. 309 * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner 310 * 311 * \implsparsesolverconcept 312 * 313 * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations() 314 * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations 315 * and NumTraits<Scalar>::epsilon() for the tolerance. 316 * 317 * The tolerance corresponds to the relative residual error: |Ax-b|/|b| 318 * 319 * \b Performance: when using sparse matrices, best performance is achied for a row-major sparse matrix format. 320 * Moreover, in this case multi-threading can be exploited if the user code is compiled with OpenMP enabled. 321 * See \ref TopicMultiThreading for details. 322 * 323 * By default the iterations start with x=0 as an initial guess of the solution. 324 * One can control the start using the solveWithGuess() method. 325 * 326 * IDR(s) can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink. 327 * 328 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner 329 */ 330 template <typename _MatrixType, typename _Preconditioner> 331 class IDRS : public IterativeSolverBase<IDRS<_MatrixType, _Preconditioner> > 332 { 333 334 public: 335 typedef _MatrixType MatrixType; 336 typedef typename MatrixType::Scalar Scalar; 337 typedef typename MatrixType::RealScalar RealScalar; 338 typedef _Preconditioner Preconditioner; 339 340 private: 341 typedef IterativeSolverBase<IDRS> Base; 342 using Base::m_error; 343 using Base::m_info; 344 using Base::m_isInitialized; 345 using Base::m_iterations; 346 using Base::matrix; 347 Index m_S; 348 bool m_smoothing; 349 RealScalar m_angle; 350 bool m_residual; 351 352 public: 353 /** Default constructor. */ 354 IDRS(): m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {} 355 356 /** Initialize the solver with matrix \a A for further \c Ax=b solving. 357 358 This constructor is a shortcut for the default constructor followed 359 by a call to compute(). 360 361 \warning this class stores a reference to the matrix A as well as some 362 precomputed values that depend on it. Therefore, if \a A is changed 363 this class becomes invalid. Call compute() to update it with the new 364 matrix A, or modify a copy of A. 365 */ 366 template <typename MatrixDerived> 367 explicit IDRS(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_S(4), m_smoothing(false), 368 m_angle(RealScalar(0.7)), m_residual(false) {} 369 370 371 /** \internal */ 372 /** Loops over the number of columns of b and does the following: 373 1. sets the tolerence and maxIterations 374 2. Calls the function that has the core solver routine 375 */ 376 template <typename Rhs, typename Dest> 377 void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const 378 { 379 m_iterations = Base::maxIterations(); 380 m_error = Base::m_tolerance; 381 382 bool ret = internal::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S,m_smoothing,m_angle,m_residual); 383 384 m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence; 385 } 386 387 /** Sets the parameter S, indicating the dimension of the shadow space. Default is 4*/ 388 void setS(Index S) 389 { 390 if (S < 1) 391 { 392 S = 4; 393 } 394 395 m_S = S; 396 } 397 398 /** Switches off and on smoothing. 399 Residual smoothing results in monotonically decreasing residual norms at 400 the expense of two extra vectors of storage and a few extra vector 401 operations. Although monotonic decrease of the residual norms is a 402 desirable property, the rate of convergence of the unsmoothed process and 403 the smoothed process is basically the same. Default is off */ 404 void setSmoothing(bool smoothing) 405 { 406 m_smoothing=smoothing; 407 } 408 409 /** The angle must be a real scalar. In IDR(s), a value for the 410 iteration parameter omega must be chosen in every s+1th step. The most 411 natural choice is to select a value to minimize the norm of the next residual. 412 This corresponds to the parameter omega = 0. In practice, this may lead to 413 values of omega that are so small that the other iteration parameters 414 cannot be computed with sufficient accuracy. In such cases it is better to 415 increase the value of omega sufficiently such that a compromise is reached 416 between accurate computations and reduction of the residual norm. The 417 parameter angle =0.7 (”maintaining the convergence strategy”) 418 results in such a compromise. */ 419 void setAngle(RealScalar angle) 420 { 421 m_angle=angle; 422 } 423 424 /** The parameter replace is a logical that determines whether a 425 residual replacement strategy is employed to increase the accuracy of the 426 solution. */ 427 void setResidualUpdate(bool update) 428 { 429 m_residual=update; 430 } 431 432 }; 433 434 } // namespace Eigen 435 436 #endif /* EIGEN_IDRS_H */ 437