quicwind/multigrid.m

   1 function [x,myrelres]=multigrid(A,b,x0,M,P,R,p)
   2 % x=multigrid(A,b,P,R)
   3 % solve A{1}*x = b by multigrid method
   4 % In:
   5 %   A   cell vector length N of matrix-vector multiplication functions
   6 %   b   right hand side
   7 %   x0  initial approximation
   8 %   M   cell vector length N of preconditioner functions
   9 %   P   cell vector length N-1 of prolongation functions (interpolation)
  10 %   R   cell vector length N-1 of restriction functions (averaging)
  11 %   p   structure with parameters
  12
  13 N=length(A);
  14 r=cell(1,N+1);
  15 r{1}=b(:);
  16 if isempty(x0),
  17     u{1}=zeros(size(r{1}));
  18 else
  19     u{1}=x0(:);
  20 end
  21 for it=1:p.maxit
  22     myrelres=norm(r{1}-A{1}(u{1}))/max(norm(r{1}),realmin);
  23     fprintf('multigrid iteration %i start relative residual %g\n',it,myrelres)
  24     mglevel(1)
  25     myrelres=norm(r{1}-A{1}(u{1}))/max(norm(r{1}),realmin);
  26     fprintf('multigrid iteration %i end relative residual %g\n',it,myrelres)
  27     if myrelres(end)<p.tol
  28         fprintf('Multigrid iterations converged to given tolerance %g\n',p.tol)
  29         break
  30     elseif it>=p.maxit
  31         fprintf('Multigrid iterations did not converge to given tolerance %g\n',p.tol)
  32     end
  33 end
  34 x = reshape(u{1},size(b));
  35
  36 function mglevel(l)
  37     % multigrid on level l, all matrices and vectors are global
  38     [u{l},flag,relres,iter,resvec] = pcg(A{l},r{l},p.tolsm1,p.maxsm1,M{l},[],u{l});% pre-smoothing
  39     fprintf('level %i relative residuals ',l);  print_in_line(resvec); fprintf('\n')
  40     r{l+1} = R{l}(r{l} - A{l}(u{l}));  % set up coarse problem
  41     u{l+1} = zeros(size(r{l+1}));
  42     if l < N-1
  43         for i=1:p.ncoarse
  44             mglevel(l+1);  % iterate recursively on coarse level
  45         end
  46     else    % coarsest level
  47         [u{l+1},flag,relres,iter,resvec] = pcg(A{l+1},r{l+1},p.tolcr,p.maxcr,M{l+1},[],u{l+1});
  48         fprintf('level %i relative residuals ',l+1); print_in_line(resvec); fprintf('\n')
  49     end
  50     u{l} = u{l} + P{l}(u{l+1}); % apply coarse correction
  51     [u{l},flag,relres,iter,resvec] = pcg(A{l},r{l},p.tolsm2,p.maxsm2,M{1},[],u{l});% post-smoothing
  52     fprintf('level %i relative residuals ',l); print_in_line(resvec); fprintf('\n')
  53 end
  54 end
  55
  56 function print_in_line(a)
  57 for i=1:length(a)
  58     fprintf(' %6.3g',a(i))
  59 end
  60 end
  61