cg.m

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

  function [uh,residnormvec] = pcg_nash(khat,bh,Active,alpha,L,cg_params)

%  [uh,residnormvec] = pcg_nash(khat,bh,Active,alpha,L,cg_params)
%
%  Solve linear system Hbar*u = b using conjugate gradient iteration.

global cgplotfigsize;
%  CG initialization.

  [nx,ny] = size(bh);
  uh = zeros(nx,ny);
  resid = bh;
  residnormvec = norm(resid(:));
  stepnormvec = [];
  cgiter = 0;
  stop_flag = 0;
  
%  Get CG iteration parameters.
  
  maxiter = cg_params(1);
  steptol = cg_params(2);
  residtol= cg_params(3);
  cg_tab_flag = cg_params(4);
  cg_fig = cg_params(5);

%  CG iterations.
  
  while  stop_flag == 0
    cgiter = cgiter + 1;
    rd = residnormvec(cgiter)^2;

    %  Compute conjugate direction ph.

    if cgiter == 1,
       ph = resid; 
     else
       betak = rd / rdlast;
       ph = resid + betak * ph;
    end

    %  Update uh and residual.
    
    Hbar_ph = mult_hbar(ph,khat,alpha,L,Active);
    alphak = rd / (ph(:)'*Hbar_ph(:));
    uh = uh + alphak*ph;
    resid = resid - alphak*Hbar_ph;
    rdlast = rd;

    residnorm = norm(resid(:));
    stepnorm = abs(alphak)*norm(ph(:))/norm(uh(:));
    residnormvec = [residnormvec; residnorm];
    stepnormvec = [stepnormvec; stepnorm];
 
    %  Check stopping criteria.
    
    if cgiter >= maxiter
      stop_flag = 1;
    elseif stepnorm < steptol
      stop_flag = 2;
    elseif residnorm / residnormvec(1) < residtol
      stop_flag = 3;
    end
    
    %  Output convergence information.
   
    if cg_tab_flag == 1
      fprintf('  CG iter%3.0f, ||resid||=%6.4e, ||step||=%6.4e \n', ... 
         cgiter, residnormvec(cgiter), stepnormvec(cgiter));
    end
    if cg_fig > 0
      figure(cg_fig)
      set(cg_fig,...
        'units','pixels',...
        'pos',cgplotfigsize ,...
        'numbertitle','off',...
        'name','"CG Performance"');

        subplot(221)
          semilogy(residnormvec/residnormvec(1),'o')
          xlabel('CG iteration')
          title('CG residual norm')
        subplot(222)
          semilogy(stepnormvec,'o')
          xlabel('CG iteration')
          title('CG relative step norm')
      drawnow
    end

  end

global figure_list;
uicontrol(cg_fig, ...
        'style','push', ...
        'callback','close(gcf),figure_list(7) = -1;', ...
        'string','Close');