1######################################################################## 2## 3## Copyright (C) 2016-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{x} =} tfqmr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0}, @dots{}) 28## @deftypefnx {} {@var{x} =} tfqmr (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M}, [], @var{x0}, @dots{}) 29## @deftypefnx {} {[@var{x}, @var{flag}, @var{relres}, @var{iter}, @var{resvec}] =} tfqmr (@var{A}, @var{b}, @dots{}) 30## Solve @code{A x = b} using the Transpose-Tree qmr method, based on the cgs. 31## 32## The input parameters are: 33## 34## @itemize @minus 35## 36## @item @var{A} is the matrix of the linear system and it must be square. 37## @var{A} can be passed as a matrix, function handle, or inline 38## function @code{Afun} such that @code{Afun(x) = A * x}. Additional 39## parameters to @code{Afun} are passed after @var{x0}. 40## 41## @item @var{b} is the right hand side vector. It must be a column vector 42## with the same number of rows as @var{A}. 43## 44## @item @var{tol} is the relative tolerance, if not given or set to [] the 45## default value 1e-6 is used. 46## 47## @item @var{maxit} the maximum number of outer iterations, if not given or 48## set to [] the default value @code{min (20, numel (b))} is used. To be 49## compatible, since the method as different behaviors in the iteration 50## number is odd or even, is considered as iteration in @code{tfqmr} the 51## entire odd-even cycle. That is, to make an entire iteration, the algorithm 52## performs two sub-iterations: the odd one and the even one. 53## 54## @item @var{M1}, @var{M2} are the preconditioners. The preconditioner 55## @var{M} is given as @code{M = M1 * M2}. 56## Both @var{M1} and @var{M2} can be passed as a matrix or as a function 57## handle or inline function @code{g} such that @code{g(x) = M1 \ x} or 58## @code{g(x) = M2 \ x}. 59## The technique used is the right-preconditioning, i.e., it is solved 60## @code{A*inv(M)*y = b} and then @code{x = inv(M)*y}, instead of 61## @code{A x = b}. 62## 63## @item @var{x0} the initial guess, if not given or set to [] the default 64## value @code{zeros (size (b))} is used. 65## 66## @end itemize 67## 68## The arguments which follow @var{x0} are treated as parameters, and passed in 69## a proper way to any of the functions (@var{A} or @var{M}) which are passed 70## to @code{tfqmr}. 71## 72## The output parameters are: 73## 74## @itemize @minus 75## 76## @item @var{x} is the approximation computed. If the method doesn't 77## converge then it is the iterated with the minimum residual. 78## 79## @item @var{flag} indicates the exit status: 80## 81## @itemize @minus 82## @item 0: iteration converged to the within the chosen tolerance 83## 84## @item 1: the maximum number of iterations was reached before convergence 85## 86## @item 2: the preconditioner matrix is singular 87## 88## @item 3: the algorithm reached stagnation 89## 90## @item 4: the algorithm can't continue due to a division by zero 91## @end itemize 92## 93## @item @var{relres} is the relative residual obtained as 94## @code{(@var{A}*@var{x}-@var{b}) / @code{norm (@var{b})}}. 95## 96## @item @var{iter} is the iteration which @var{x} is 97## computed. 98## 99## @item @var{resvec} is a vector containing the residual at each iteration 100## (including @code{norm (b - A x0)}). 101## Doing @code{length (@var{resvec}) - 1} is possible to see the 102## total number of iterations performed. 103## 104## @end itemize 105## 106## Let us consider a trivial problem with a tridiagonal matrix 107## 108## @example 109## @group 110## n = 20; 111## A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... 112## toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... 113## sparse (1, 2, 1, 1, n) * n / 2); 114## b = A * ones (n, 1); 115## restart = 5; 116## [M1, M2] = ilu (A); # in this tridiag case it corresponds to chol (A)' 117## M = M1 * M2; 118## Afun = @@(x) A * x; 119## Mfun = @@(x) M \ x; 120## M1fun = @@(x) M1 \ x; 121## M2fun = @@(x) M2 \ x; 122## @end group 123## @end example 124## 125## @sc{Example 1:} simplest usage of @code{tfqmr} 126## 127## @example 128## x = tfqmr (A, b, [], n) 129## @end example 130## 131## @sc{Example 2:} @code{tfqmr} with a function which computes 132## @code{@var{A} * @var{x}} 133## 134## @example 135## x = tfqmr (Afun, b, [], n) 136## @end example 137## 138## @sc{Example 3:} @code{tfqmr} with a preconditioner matrix @var{M} 139## 140## @example 141## x = tfqmr (A, b, [], 1e-06, n, M) 142## @end example 143## 144## @sc{Example 4:} @code{tfqmr} with a function as preconditioner 145## 146## @example 147## x = tfqmr (Afun, b, 1e-6, n, Mfun) 148## @end example 149## 150## @sc{Example 5:} @code{tfqmr} with preconditioner matrices @var{M1} 151## and @var{M2} 152## 153## @example 154## x = tfqmr (A, b, [], 1e-6, n, M1, M2) 155## @end example 156## 157## @sc{Example 6:} @code{tfmqr} with functions as preconditioners 158## 159## @example 160## x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun) 161## @end example 162## 163## @sc{Example 7:} @code{tfqmr} with as input a function requiring an argument 164## 165## @example 166## @group 167## function y = Ap (A, x, z) # compute A^z * x 168## y = x; 169## for i = 1:z 170## y = A * y; 171## endfor 172## endfunction 173## Apfun = @@(x, string, p) Ap (A, x, string, p); 174## x = tfqmr (Apfun, b, [], [], [], [], [], 2); 175## @end group 176## @end example 177## 178## @sc{Example 8:} explicit example to show that @code{tfqmr} uses a 179## right preconditioner 180## 181## @example 182## @group 183## [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed 184## M = M1 * M2; 185## 186## ## reference solution computed by tfqmr after one iteration 187## [x_ref, fl] = tfqmr (A, b, [], 1, M) 188## 189## ## right preconditioning 190## [y, fl] = tfqmr (A / M, b, [], 1) 191## x = M \ y # compare x and x_ref 192## 193## @end group 194## @end example 195## 196## Reference: 197## 198## @nospell{Y. Saad}, @cite{Iterative Methods for Sparse Linear Systems}, 199## Second edition, 2003, SIAM 200## 201## @seealso{bicg, bicgstab, cgs, gmres, pcg, qmr, pcr} 202## 203## @end deftypefn 204 205function [x_min, flag, relres, iter_min, resvec] = ... 206 tfqmr (A, b, tol = [], maxit = [], M1 = [], M2 = [], ... 207 x0 = [], varargin) 208 209 [Afun, M1fun, M2fun] = __alltohandles__ (A, b, M1, M2, "tfqmr"); 210 211 [tol, maxit, x0] = __default__input__ ({1e-06, 2 * min(20, rows (b)), ... 212 zeros(rows (b), 1)}, tol, ... 213 maxit, x0); 214 215 maxit = 2 * maxit; # To be compatible, since iteration = odd+even ones 216 217 norm_b = norm (b, 2); 218 if (norm_b == 0) 219 if (nargout < 2) 220 printf("The right hand side vector is all zero so tfqmr \n") 221 printf ("returned an all zero solution without iterating.\n") 222 endif 223 x_min = zeros (numel (b), 1); 224 iter_min = 0; 225 flag = 0; 226 resvec = 0; 227 relres = 0; 228 return 229 endif 230 231 x = x_pr = x_min = x0; 232 iter = iter_min = m = 0; 233 resvec = zeros (maxit, 1); 234 flag = 1; 235 236 w = u = r = r_star = b - feval (Afun, x0, varargin{:}); 237 rho_1 = (r_star' * r); 238 d = 0; 239 tau = norm (r, 2); 240 theta = eta = 0; 241 resvec (1, 1) = norm (r, 2); 242 it = 1; 243 244 try 245 warning("error", "Octave:singular-matrix", "local"); 246 u_hat = feval (M1fun, u, varargin{:}); 247 u_hat = feval (M2fun, u_hat, varargin{:}); 248 v = feval (Afun, u_hat, varargin{:}); 249 catch 250 flag = 2; 251 end_try_catch 252 while ((flag != 2) && (iter < maxit) && ... 253 (resvec (iter + 1, 1) >= norm_b * tol)) 254 if (it > 0) # iter is even 255 v_r = r_star' * v; # inner prod between r_star and v 256 if (v_r == 0) 257 ## Essentially the next iteration doesn't change x, 258 ## and the iter after this will have a division by zero 259 flag = 4; 260 break 261 endif 262 alpha = rho_1 / v_r; 263 u_1 = u - alpha * v; # u at the after iteration 264 endif 265 u_hat = feval (M1fun, u, varargin{:}); 266 u_hat = feval (M2fun, u_hat, varargin{:}); 267 w -= alpha * feval (Afun, u_hat, varargin{:}); 268 d = u_hat + ((theta * theta) / alpha) * eta * d; 269 theta = norm (w, 2) / tau; 270 c = 1 / sqrt (1 + theta * theta); 271 tau *= theta * c; 272 eta = (c * c) * alpha; 273 x += eta * d; 274 r -= eta * feval (Afun, d, varargin{:}); 275 if (it < 0) # iter is odd 276 rho_2 = rho_1; 277 rho_1 = (r_star' * w); 278 if (rho_1 == 0) 279 ## Essentially the next iteration doesn't change x, 280 ## and the iter after this will have a division by zero 281 flag = 4; 282 break 283 endif 284 beta = rho_1 / rho_2; 285 u_1 = w + beta * u; # u at the after iteration 286 u1_hat = feval (M1fun, u_1, varargin{:}); 287 u1_hat = feval (M2fun, u1_hat, varargin{:}); 288 v = feval (Afun, u1_hat, varargin{:}) + ... 289 beta * (feval (Afun, u_hat, varargin{:}) + beta * v); 290 endif 291 u = u_1; 292 iter += 1; 293 resvec (iter + 1, 1) = norm (r, 2); 294 if (resvec (iter + 1, 1) <= resvec (iter_min + 1, 1)) 295 ## iter with min residual 296 x_min = x; 297 iter_min = iter; 298 endif 299 if (norm (x_pr - x) <= norm (x) * eps) 300 flag = 3; # Stagnation 301 break 302 endif 303 x_pr = x; 304 it = -it; 305 endwhile 306 resvec = resvec (1: (iter + 1)); 307 308 relres = resvec (iter_min + 1) / norm (b); 309 iter_min = floor(iter_min / 2); # compatibility, since it 310 # makes two times the effective iterations 311 312 if (relres <= tol) 313 flag = 0; 314 endif 315 316 if (nargout < 2) # Output strings 317 switch (flag) 318 case {0} 319 printf ("tfqmr converged at iteration %i ", iter_min); 320 printf ("to a solution with relative residual %e\n", relres); 321 case {1} 322 printf ("tfqmr stopped at iteration %i ", iter); 323 printf ("without converging to the desired tolerance %e\n", tol); 324 printf ("because the maximum number of iterations was reached.\n"); 325 printf ("The iterate returned (number %i) ", iter_min); 326 printf ("has relative residual %e\n", relres); 327 case {2} 328 printf ("tfqmr stopped at iteration %i ", iter); 329 printf ("without converging to the desired tolerance %e\n", tol); 330 printf ("because the preconditioner matrix is singular.\n"); 331 printf ("The iterate returned (number %i) ", iter_min); 332 printf ("has relative residual %e\n", relres); 333 case {3} 334 printf ("tfqmr stopped at iteration %i ", iter); 335 printf ("without converging to the desired tolerance %e\n", tol); 336 printf ("because the method stagnated.\n"); 337 printf ("The iterate returned (number %i) ", iter_min); 338 printf ("has relative residual %e\n", relres); 339 case {4} 340 printf ("tfqmr stopped at iteration %i ", iter); 341 printf ("without converging to the desired tolerance %e\n", tol); 342 printf ("because the method can't continue.\n"); 343 printf ("The iterate returned (number %i) ", iter_min); 344 printf ("has relative residual %e\n", relres); 345 endswitch 346 endif 347endfunction 348%!test 349%! ## Check that all type of inputs work 350%! A = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); 351%! b = sum (A, 2); 352%! M1 = diag (sqrt (diag (A))); 353%! M2 = M1; 354%! maxit = 10; 355%! Afun = @(z) A * z; 356%! M1_fun = @(z) M1 \ z; 357%! M2_fun = @(z) M2 \ z; 358%! [x, flag] = tfqmr (A,b); 359%! assert (flag, 0); 360%! [x, flag] = tfqmr (A, b, [], maxit, M1, M2); 361%! assert (flag, 0); 362%! [x, flag] = tfqmr (A, b, [], maxit, M1_fun, M2_fun); 363%! assert (flag, 0); 364%! [x, flag] = tfqmr (A, b, [], maxit, M1_fun, M2); 365%! assert (flag, 0); 366%! [x, flag] = tfqmr (A, b, [], maxit, M1, M2_fun); 367%! assert (flag, 0); 368%! [x, flag] = tfqmr (Afun, b); 369%! assert (flag, 0); 370%! [x, flag] = tfqmr (Afun, b, [], maxit, M1, M2); 371%! assert (flag, 0); 372%! [x, flag] = tfqmr (Afun, b, [], maxit, M1_fun, M2); 373%! assert (flag, 0); 374%! [x, flag] = tfqmr (Afun, b, [], maxit, M1, M2_fun); 375%! assert (flag, 0); 376%! [x, flag] = tfqmr (Afun, b, [], maxit, M1_fun, M2_fun); 377%! assert (flag, 0); 378 379%!shared A, b, n, M1, M2 380%! 381%!test 382%! n = 100; 383%! A = spdiags ([-2*ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); 384%! b = sum (A, 2); 385%! tol = 1e-8; 386%! maxit = 15; 387%! M1 = spdiags ([ones(n,1)/(-2) ones(n,1)],-1:0, n, n); 388%! M2 = spdiags ([4*ones(n,1) -ones(n,1)], 0:1, n, n); 389%! [x, flag, relres, iter, resvec] = tfqmr (A, b, tol, maxit, M1, M2); 390%! assert (x, ones (size (b)), 1e-7); 391%! 392%!test 393%!function y = afun (x, a) 394%! y = a * x; 395%!endfunction 396%! 397%! tol = 1e-8; 398%! maxit = 15; 399%! 400%! [x, flag, relres, iter, resvec] = tfqmr (@(x) afun (x, A), b, 401%! tol, maxit, M1, M2); 402%! assert (x, ones (size (b)), 1e-7); 403 404%!test 405%! ## Jacobi preconditioner works 406%! n = 10; 407%! tol = 1e-8; 408%! A = hilb (n) + 1i * hilb (n); 409%! A(1,1) = 100; 410%! A(n, n) = 100; 411%! b = sum (A, 2); 412%! [x, flag, relres, iter, resvec] = tfqmr (A, b, tol); 413%! assert (x, ones (size (b)), 0.005); 414%! assert (iter, 8); 415%! [x, flag, relres, iter, resvec] = tfqmr (A, b, tol, [], diag (diag (A))); 416%! assert (x, ones (size (b)), 0.002); 417%! assert (iter, 6); 418 419%!test 420%! ## Solve complex linear system 421%! A = [1 + 1i, 1 + 1i; 2 - 1i, 2 + 1i]; 422%! b = A * [1; 1]; 423%! [x, flag, relres, iter, resvec] = tfqmr (A, b, [], 3); 424%! assert (x, [1; 1], 1e-6); 425 426%!test 427%! A = diag (1:50); 428%! A (1,50) = 10000; 429%! b = ones (50,1); 430%! [x, flag, relres, iter, resvec] = tfqmr (A, b, [], 100); 431%! assert (flag, 0) 432%! assert (x, A \ b, 1e-05) 433%! ## Detects a singular preconditioner 434%! M = ones (50); 435%! M(1, 1) = 0; 436%! [x, flag] = tfqmr (A, b, [], 100, M); 437%! assert (flag, 2) 438 439%!test 440%! A = single (1); 441%! b = 1; 442%! [x, flag] = tfqmr (A, b); 443%! assert (class (x), "single") 444 445%!test 446%! A = 1; 447%! b = single (1); 448%! [x, flag] = tfqmr (A, b); 449%! assert (class (x), "single") 450 451%!test 452%! A = single (1); 453%! b = single (1); 454%! [x, flag] = tfqmr (A, b); 455%! assert (class (x), "single") 456 457%!test 458%!function y = Afun (x) 459%! A = toeplitz ([2, 1, 0, 0], [2, -1, 0, 0]); 460%! y = A * x; 461%!endfunction 462%! [x, flag] = tfqmr ("Afun", [1; 2; 2; 3]); 463%! assert (x, ones(4, 1), 1e-6) 464 465%!test # unpreconditioned residual 466%! A = toeplitz (sparse ([2, 1, 0, 0, 0]), sparse ([2, -1, 0, 0, 0])); 467%! b = sum (A, 2); 468%! M = magic (5); 469%! [x, flag, relres] = tfqmr (A, b, [], 3, M); 470%! assert (relres, norm (b - A * x) / norm (b), 8 * eps) 471 472%!demo # simplest use 473%! n = 20; 474%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... 475%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... 476%! sparse (1, 2, 1, 1, n) * n / 2); 477%! b = A * ones (n, 1); 478%! [M1, M2] = ilu (A + 0.1 * eye (n)); 479%! M = M1 * M2; 480%! x = tfqmr (A, b, [], n); 481%! Afun = @(x) A * x; 482%! x = tfqmr (Afun, b, [], n); 483%! x = tfqmr (A, b, 1e-6, n, M); 484%! x = tfqmr (A, b, 1e-6, n, M1, M2); 485%! Mfun = @(z) M \ z; 486%! x = tfqmr (Afun, b, 1e-6, n, Mfun); 487%! M1fun = @(z) M1 \ z; 488%! M2fun = @(z) M2 \ z; 489%! x = tfqmr (Afun, b, 1e-6, n, M1fun, M2fun); 490%! function y = Ap (A, x, z) # compute A^z * x or (A^z)' * x 491%! y = x; 492%! for i = 1:z 493%! y = A * y; 494%! endfor 495%! endfunction 496%! Afun = @(x, p) Ap (A, x, p); 497%! x = tfqmr (Afun, b, [], 2*n, [], [], [], 2); # solution of A^2 * x = b 498 499%!demo 500%! n = 10; 501%! A = toeplitz (sparse ([1, 1], [1, 2], [2, 1] * n ^ 2, 1, n)) + ... 502%! toeplitz (sparse (1, 2, -1, 1, n) * n / 2, ... 503%! sparse (1, 2, 1, 1, n) * n / 2); 504%! b = A * ones (n, 1); 505%! [M1, M2] = ilu (A + 0.3 * eye (n)); # factorization of A perturbed 506%! M = M1 * M2; 507%! 508%! ## reference solution computed by tfqmr after one iteration 509%! [x_ref, fl] = tfqmr (A, b, [], 1, M); 510%! x_ref 511%! 512%! ## right preconditioning 513%! [y, fl] = tfqmr (A / M, b, [], 1); 514%! x = M \ y # compare x and x_ref 515