precond_jac.m

来自「采用matlab编写的数字图像恢复程序」· M 代码 · 共 38 行

  function [u,Qresidnormvec] = ...
    precond_jac(b,AH,p0,dcoef_xh,dcoef_yh,alpha,Qpcg_iter,nu)

%  [u,Qresidnormvec]=precond_jac(b,AH,p0,dcoef_xh,dcoef_yh,alpha,Qpcg_iter,nu) 
%
%  Apply Multilevel Block Jacobi preconditioner to system Ah*u = b,
%  Conceptually,
%     u =    Pu      +    Qu
%       = inv(P*Ah*P + Q*Ah*Q) * b
%       = pinv(P*Ah*P)*(P*b) + pinv(Q*Ah*Q)*(Q*b),
%  where P is the projection onto the coarse subspace and Q=I-P is
%  the projection onto its orthogonal complement.
%  The first term Pu is computed by
%        prolong( inv(AH) * restrict(b)),
%  where AH is a coarse grid approximation to Ah, and prolong and 
%  restrict denote intergrid transfer operators.
%  The second term Qu is approximated by
%        pinv(Q*(alpha*L)*Q) * (Q*b)

  n = max(size(b));
  h = 1/n;
  n0 = 2^p0;
  p = log2(n);
  bH = restrict(b,p-p0);
  Qb = b - prolong(bH,p-p0);

  uH = reshape(AH \ bH(:),n0,n0);
  Pu = prolong(uH,p-p0);

%  Solve Q*(alpha*L)*Qu = Q*b using a projected multigrid/PCG scheme.

  [Qu,Qresidnormvec] = QLQsolve(Qb,p0,dcoef_xh,dcoef_yh,Qpcg_iter,nu); 
  Qu = (1/alpha) * Qu;

%

  u = Pu + Qu;