%|^^^^^^^^^^^/^^^^^^|_|_|_|^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^|
%|richa_nr = nr_1.m |_|_|_| Preconditioned linear solver 24_6_2000   |
%|___________/______|_|_|_|__________________________________________|
%|Features: 1/ Repeated 1-level DWT for k times (S=DWT)              |
%|          2/ Schur Complement Method  (rich for Richardson)        |
%|          3/ Multilevel Setting Direct Solver on the Coarsest Level|
%|___________/______|_|_|_|__________________________________________|
% USAGE: [xo,r1,mit] = nr_richard(R1,x1,b1,max_it,tol)
function [xo,res,mit] = nr_richard(R1,x1,b1,max_it,tol)
  if nargin<1,help richa_nr; return; end
  warning off MATLAB:flops:UnavailableFunction
% Input (see nr_gmres.m)
%
%      A=R1:      Linear system to be solved
%      x=x1:      Initial guess
%      b=b1:      Right-hand side
%
% subroutine :    richa_up.m

        global mu n_s levs order prob
iout=0; res=[];
N =length(R1); Lev=floor(log(N)/log(2)); space='               '; 
maxi = floor(log(N)/log(10));               k = levs; 
n_tot=max_it;

%%% ________ %%% _______________________________________________________%
%%% Step (1) %%% Input ++++++++++++++++++++++++++++++++++++++++++++++++++

disp(['---Richardson richa_nr.m:  N = ' int2str(N) ' = 2^Lev = 2^' int2str(Lev) ... \
      '------------'])

%%% ________ %%% _______________________________________________________%
%%% Step (2) %%% Make variables +++++++++++++++++++++++++++++++++++++++++

    R=[]; A=[];L=[];U=[];B=[];C=[]; b=[];x=[];z=[];  %M-level matrices/vectors
for i=1:k+2
    id = floor(log(i)/log(10));
  if i<k+2
    s=sprintf('A%d%s',i,space(1:maxi-id)); A=[A;s];
    s=sprintf('B%d%s',i,space(1:maxi-id)); B=[B;s];
    s=sprintf('C%d%s',i,space(1:maxi-id)); C=[C;s];
  end
    s=sprintf('R%d%s',i,space(1:maxi-id)); R=[R;s];
    s=sprintf('L%d%s',i,space(1:maxi-id)); L=[L;s];
    s=sprintf('U%d%s',i,space(1:maxi-id)); U=[U;s];
    s=sprintf('b%d%s',i,space(1:maxi-id)); b=[b;s];
    s=sprintf('x%d%s',i,space(1:maxi-id)); x=[x;s];
    s=sprintf('z%d%s',i,space(1:maxi-id)); z=[z;s];
end % No more making arrays

%%% ________ %%% _______________________________________________________%
%%% Step (3) %%% Matrix one-off transformations +++++++++++++++++++++++++


  thr=1.0/N^2; n_e=[];n_o=[]; n_e(1)=N;  n_o(1)=N;  % thr=thresh
for i=2:k+2
  n_e(i) = floor( n_o(i-1)/2 );  %% n x n is the full matrix at LEV = i
  n_o(i) = n_o(i-1) - n_e(i);
end

for i=1:k+1
    ni = n_e(i+1); nn = n_o(i);   % nxn the full matrix & ni=nn/2
 eval([R(i,:) ' = fwts(' R(i,:) ',order,1,1,1,1);']) 
    if i>1
      eval([x(i,:) ' = zeros(nn,1);']) % Initialize x_j's
    end
    if i==1
      M0 = R1;  %% extra --- saved to compute residual  \tilde{A} 
      clear ba
      ba = banded(R1(1:ni,1:ni),thr); ba=min(ba,16); ba=max(ba,1);
      fprintf('    Lev=%d mu=%d Thresh=%e  give  bandw=ba=%d\n',k,mu,thr,ba);
    end
 eval([A(i,:) ' = ba_cut(' R(i,:) '(1:ni,1:ni),    ba);'])  % A_bar
 eval([B(i,:) ' = ba_cut(' R(i,:) '(1:ni,ni+1:nn), ba);'])  % B_bar
 eval([C(i,:) ' = ba_cut(' R(i,:) '(1+ni:nn,1:ni), ba);'])  % C_bar
 eval([R(i+1,:) '=' R(i,:) '(1+ni:nn,1+ni:nn);'])  % Mr T_i+1 (full)
 eval([R(i,:) '= [' A(i,:) B(i,:) ';' C(i,:) R(i+1,:) ']-' R(i,:) ';'])
 eval(['[' L(i,:) U(i,:) ']=lu(' A(i,:) ');'])  % A's factorised
 eval(['clear ' A(i,:) ])  % A_i no longer needed cos of L U's
end % done i
 i = levs+2;
 eval(['[' L(i,:) U(i,:) ']=lu(' R(i,:) ');'])  % A's factorised
 eval(['clear ' R(i,:) ])    % A_i no longer needed cos of L U's

%%% ________ %%% ____________
%%% Step (4) %%% Main Loops +

      iter = 0; check=zeros(k+2,1);    

while iter<n_tot % Mr Global ********************************************

      i = 1; 
%% richa_up;         %%%  up to i=k+1 (Coarsest Lev) one-off
%++++++++++++++++++++++++++++++
%%%%  richa_up;     % up to i=k+1 (Coarsest Lev) recursive now (see test.m)
%%% Step (1) %%% Upward to i=k+1 the coarsest level

while i <= k+1 % *****************************************************
   if iout==1
   fprintf('Upward i k+1 = %d %d [%d / %d @ %d]\n', i,k+1, check(i),mu,iter)
   end

   if check(i)==0
      eval([b(i,:) '=fwts(' b(i,:) ',order,1,0,1,1);']);
      if i==1 & iter==0, 
       r0=eval(b(i,:));  %% extra --- saved \tilde{b} with \tilde{A}
       res=norm(r0);
      end %% 
   end
    ni = n_e(i+1); nn = n_o(i);

   eval([z(i,:) '=' R(i,:) '*' x(i,:) '+' b(i,:) ';']);
   eval([z(i,:) '(1:ni)=' U(i,:) '\(' L(i,:) '\' ... \
                             z(i,:) '(1:ni) );']);
   eval([z(i,:) '(1+ni:nn)=' z(i,:) '(1+ni:nn) -' ... \
                  C(i,:) '*' z(i,:) '(1:ni);']);
   eval([b(i+1,:) '=' C(i,:) '*(' U(i,:) '\ (' ... \
                      L(i,:) '\(' ... \
                      B(i,:) '*' x(i,:) '(1+ni:nn) ))) +' ...\
                      z(i,:) '(1+ni:nn) ;'])
   i = i + 1;
end % while i ********************************************************
%%% ________ %%% __________________________________
%%% Step (2) %%% Direct Solver on Level i=k+2 %%%
   eval([x(i,:) '=' U(i,:) '\(' L(i,:) '\' b(i,:) ');']);
   if iout==1
   fprintf('DIRECT i k+1 = %d %d [%d / %d @ %d]\n', i,k+1, check(i),mu,iter)
   end
%++++++++++++++++++++++++++++++ end of richa_up

  while i>1  %%% Level Schur Method starts  * * * * * * * * * * * * * * *

   i = i - 1;
    ni = n_e(i+1); nn = n_o(i);   

      if i>=k+1                    % i=k+1 from no DWT
        eval([x(i,:) '(1+ni:nn)=' x(i+1,:) ';'])
      else
        eval([x(i,:) '(1+ni:nn)=iwts(' x(i+1,:) ',order,1,0,1,1);']);
      end
      eval([x(i,:) '(1:ni)=' z(i,:) '(1:ni)-' U(i,:) '\(' ... \
            L(i,:) '\(' B(i,:) '*' x(i,:) '(1+ni:nn)) );'])

   if iout==1
   fprintf(' MAIN  i k+1 = %d %d [%d / %d @ %d]\n', i,k+1, check(i),mu,iter)
   end

   if i==1
      check(i) = mu;  % treat i=1 specially (not necessory)
      iter = iter + 1;
   else
      check(i) = check(i) + 1;
   end

   if check(i)<mu %  *   *   *   *   *   *   *   *   *   *   *   *   *   *
%++++++++++++++++++++++++++++++
%%%%  richa_up;     % up to i=k+1 (Coarsest Lev) recursive now (see test.m)
%%% Step (1) %%% Upward to i=k+1 the coarsest level

while i <= k+1 % *****************************************************
   if iout==1
   fprintf('Upward i k+1 = %d %d [%d / %d @ %d]\n', i,k+1, check(i),mu,iter)
   end

   if check(i)==0
      eval([b(i,:) '=fwts(' b(i,:) ',order,1,0,1,1);']);
      if i==1 & iter==0, 
       r0=eval(b(i,:));  %% extra --- saved \tilde{b} with \tilde{A}
      end %% 
   end
    ni = n_e(i+1); nn = n_o(i);

   eval([z(i,:) '=' R(i,:) '*' x(i,:) '+' b(i,:) ';']);
   eval([z(i,:) '(1:ni)=' U(i,:) '\(' L(i,:) '\' ... \
                             z(i,:) '(1:ni) );']);
   eval([z(i,:) '(1+ni:nn)=' z(i,:) '(1+ni:nn) -' ... \
                  C(i,:) '*' z(i,:) '(1:ni);']);
   eval([b(i+1,:) '=' C(i,:) '*(' U(i,:) '\ (' ... \
                      L(i,:) '\(' ... \
                      B(i,:) '*' x(i,:) '(1+ni:nn) ))) +' ...\
                      z(i,:) '(1+ni:nn) ;'])
   i = i + 1;
end % while i ********************************************************
%%% ________ %%% __________________________________
%%% Step (2) %%% Direct Solver on Level i=k+2 %%%
   eval([x(i,:) '=' U(i,:) '\(' L(i,:) '\' b(i,:) ');']);
   if iout==1
   fprintf('DIRECT i k+1 = %d %d [%d / %d @ %d]\n', i,k+1, check(i),mu,iter)
   end
   else
      if i>1, check(i) = 0; end
   end % check  *   *   *   *   *   *   *   *   *   *   *   *   *   *   *
%++++++++++++++++++++++++++++++ end of richa_up

  end % while i * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

%%%%%%% Residual checking

    r1 = r0 - M0*x1;  now_Residual = norm(r1)/norm(r0);
    res=[res now_Residual];
    fprintf('richa_nr.m Steps = %2d  Residual = %10.2e\n', iter,now_Residual)
    iter_stop=iter; 
    if res(iter+1) >= 0.999*res(iter) & iter>5; break;end % Probable Failure
    if now_Residual<=tol | isnan(now_Residual), iter=n_tot; end

end % Mr Global *********************************************************

    xo = iwts(x1,order,1,0,1,1);  % res is not DWT back
    res=res(2:end);mit = length(res);

Flop = round(flops/N^2); TOL=round( log(tol)/log(10) );
s1=sprintf('|P%2d Richard_nr.m  N=2^%2d =%4d | %2d steps to TOL=%2d', ... \
            prob,Lev,N, iter_stop,TOL);
s2=sprintf(' LEV=%2d | mu/ba=%1d/%2d TOL %3d  flop %d__|', ... \
            k,mu, ba, TOL,Flop) ;
fid = fopen('wk_richa.out','a+');
  fprintf(fid,'%s %s\n',s1,s2);
fclose(fid);
%-----------------------------------------------------------------------%
   if iout==1
   Level_dimensions=[1:k+2; n_o;  n_e]
   var=who('R*')', var=who('A*')', var=who('B*')', var=who('C*')'
   end
%<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
% banded.m --- find out the total semi-bandwidth of a banded matrix
function b=banded(A,tol);
   if nargin < 2
      tol=1.0e-6;
   end
    n=max(size(A)); b=0; save=zeros(1,floor(n/2));
for row=2:n/2    % Only search lower half
    flag = 0;
for col=1:row
    if (abs(A(row,col)) >= tol & flag==0),
       b_row = col; flag=1; save(row)=abs(col-row);
    end
end
end 
b=max(save);

% ba-cut.m --- Cut out a semi-band ba out of A plus border ba if id=1
% -------- ||| If ba is -ve, try threshold with 0.01*max(ele)
function B=ba_cut(A,ba,id);
[n m] = size(A); B=speye(n,m);

if ba<0 %% Thresh
B=A;
return
end %%%%%% Thresh

mb=m-ba; nb = n-ba; if mb<0; mb=0; end,  if nb<0; nb=0; end
if n>m
 nb=nb-n+m;
elseif m>n
 mb=mb-m+n;
end
   if nargin < 2, ba=1; end
   if nargin < 3, id=1; end  % border off = 0
for row=1:n
 %_________________________diag bands
 ok1 = row-ba; ok2=row+ba; %% left right - usual  col = row-ba:row+ba;
    if ok1<1, ok1=1; end
    if ok2>m, ok2=m; end
    col = ok1:ok2;
 %_________________________bott rows
 if row>=nb+1
    col = row-ba:m;
    if row-ba<1, col=1:m; end
    if id==1
    colr = 1:m;  B(row,colr)=A(row,colr);
    end
 %_________________________righ bands
 else
    if id==1
    colr = mb+1:m;  B(row,colr)=A(row,colr);
    end
 end
    if ~isempty(row) & ~isempty(col)
    B(row,col)=A(row,col);
    end
end
