qlqsolve.m

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

  function [uh,residnormvec] = QLQsolve(b,p0,dcoef_xh,dcoef_yh,max_cg,nu)

%  [uh,residnormvec] = QLQsolve(b,p0,dcoef_xh,dcoef_yh,max_cg,nu)
%
%  Compute uh = pinv(Q*Lh*Q) * Qb  using multigrid preconditioned conjugate
%  gradient iteration.
%  Q is the projection onto the orthogonal complement of the coarse
%  grid subspace.
%  Lh = Lh(dceof) is a discretization of the diffusion operator
%      L*v = -div( dcoef grad v).

  n = max(size(b ));

  if nu > 0          %  Set up parameters for multigrid preconditioner.
    p = log2(n);
    level = p;
    pcgflag = 1;
  else               %  No MG preconditioning.
    pcgflag = 0;
  end
  
  Lh = get_L(dcoef_xh,dcoef_yh);
  uh = zeros(n,n);
  Qb = Qproject(b,p0);
  resid = Qb;
  residnormvec = norm(resid(:));
  
  for cg_iter = 1:max_cg
    
    if pcgflag == 1       %  Apply multigrid preconditioner
      d = zeros(n,n);
      [d] = vcycle(d,resid,dcoef_xh,dcoef_yh,level,p0,nu);
    else
      d = resid;
    end

    rd = resid(:)'*d(:);
    if cg_iter == 1,
       pvec = d; 
    else
       betak = rd / rdlast;
       pvec = d + betak * pvec;
    end

    Qp = Qproject(pvec,p0);
    Ap = Qproject( reshape(Lh*Qp(:),n,n), p0);
    alphak = rd / (pvec(:)'*Ap(:));
    uh = uh + alphak*pvec;
    resid = resid - alphak*Ap;
    rdlast = rd;
    residnormvec = [residnormvec; rd];
    
  end