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