function w = adi_pre(v,transflag)
%
% Implements ADI preconditioning. 
%
% NOTE: The matrices necessary for ADI preconditioning are 
% set up and stored when either adi_init_flag_GL is nonzero or 
% the problem dimension or rho changes from the previous run. 
% The main program should contain statements 
%     global H_GL V_GL adi_init_flag_GL; 
%     adi_init_flag_GL = <nonzero value>; 
% at the beginning and subsequently initialize the "horizontal" 
% and "vertical" operators H_GL and V_GL before any call to 
% adi_pre. 
%

global H_GL V_GL adi_init_flag_GL; 

persistent oldsize rho oldrho I rhoI Hpls Hmns Vpls Vmns ...
	   HplsT HmnsT VplsT VmnsT; 

nsweeps = 1; 

% NOTE: The choices of rho below are pinhead choices. 
  rho = .05;
%  rho = 1/sqrt(size(H_GL,1));

if adi_init_flag_GL~=0|size(v,1)~=oldsize|rho~=oldrho,  
  adi_init_flag_GL = 0; 
  I = speye(size(H_GL)); 
  rhoI = rho*I; 
  Hpls = rhoI + H_GL; Hmns = rhoI - H_GL;
  Vpls = rhoI + V_GL; Vmns = rhoI - V_GL; 
  HplsT = Hpls'; HmnsT = Hmns';
  VplsT = Vpls';VmnsT = Vmns';
  oldsize = size(v,1); 
  oldrho = rho; 
end

if nargin == 1 | (nargin==2 & strcmp(transflag,'notransp'))
  w = Vpls\(v + Hmns*(Hpls\v)); 
  for k = 2:nsweeps
    w = w + Vpls\(Hmns*(Hpls\(Vmns*w)));
  end
elseif strcmp(transflag,'transp')
  w = v;
  for k = 2:nsweeps
    w = w + VmnsT*(HplsT\(HmnsT*(VplsT\w))); 
  end
  w = VplsT\w;
  w = w + HplsT\(HmnsT*w); 
end
