precond_sgs.m

采用matlab编写的数字图像恢复程序

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

%  [u,Qresidnormvec] = ...
%    precond_sgs(b,AH,p0,kstark_hat,dcoef_xh,dcoef_yh,alpha,Qpcg_iter,nu);
%
%  Apply Multilevel Block Symmetric Gauss Seidel preconditioner 
%  to system Ah*u = b. This is based on the decomposition
%    Ah = P*Ah*P + P*Ah*Q + Q*Ah*P + Q*Ah*Q 
%  where P is the projection onto the coarse subspace and Q=I-P is
%  the projection onto its orthogonal complement. In block form,
%    [Pu]   [P*Ah*P   P*Ah*Q] [Pb]
%    [  ] = [               ] [  ]
%    [Qu]   [Q*Ah*P   Q*Ah*Q] [Qb]
%  where u = Pu + Qu and b = Pb + Qb. The SGS preconditioning yields
%  the approximation
%    [Pu Qu]' = inv(D + U) * (I - L*inv(L + D)) * [Pb Qb]'
%
%  Step 1: Compute Pv  = pinv(P*Ah*P) * Pb
%  Step 2: Compute Qw  = Qb - (Q*Ah*P) * Pv
%  Step 3: Compute Qu  = pinv(Q*Ah*Q) * Qw
%  Step 4: Compute Px  = Pb - (P*Ah*Q) * Qu
%  Step 5: Compute Pu  = pinv(P*Ah*P) * Px
%  The following additional approximations are made:
%  In step 1, Pv is approximated by
%        prolong( inv(AH) * restrict(b)),
%  where AH is a coarse grid approximation to Ah, and prolong and 
%  restrict denote intergrid transfer operators.
%  In step 3, Q*Ah*Q  is approximated by Q*(alpha*L)*Q.
%  In step 5, an approximation analogous to that in step 1 is used.
%
%  Steps 2 and 4 require the application of the operator Ah. 
%  If p_opt=1, Ah*v is computed approximately using a single fft.
%  If p_opt=2, Ah*v is computed exactly using a pair of fft's.

  n = max(size(b));
  n0 = 2^p0;
  p = log2(n);
  
  Lh = get_L(dcoef_xh,dcoef_yh);
  bH = restrict(b,p-p0);
  Pb = prolong(bH,p-p0);
  Qb = b - Pb;
  
%  Step 1.
  
  vH = reshape(AH \ bH(:),n0,n0);
  Pv = prolong(vH,p-p0);
  
%  Step 2.

%  AhPv = Ah * Pv;

  LhPv = reshape(Lh*Pv(:),n,n);
  AhPv = integral_op(Pv,kstark_hat) + alpha*LhPv;
  LPv = AhPv - prolong(restrict(AhPv,p-p0),p-p0);
  Qw = Qb - LPv;
 
%  Step 3. Solve Q*(alpha*L)*Qu = Qw using a projected multigrid/PCG scheme.
%  

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

  
%  Step 4.

  LQu = reshape(Lh*Qu(:),n,n);
  AhQu = integral_op(Qu,kstark_hat) + alpha*LQu;
  xH = bH - restrict(AhQu,p-p0);
  
%  Step 5.

  Pu = prolong(reshape(AH\xH(:),n0,n0),p-p0);
  
  u = Pu + Qu;