 close 
 addpath ../general;
 A = fd3d(30, 40, 10, 0, 0, 0, 0) ;
 n = size(A,1);
%%-------------------- create artificial rhs
%% rhs = A*kron(ones(n/2,1),[-1;1]);
 rhs = A * (ones(n,1)) ; 
%% random initial guess
 x0 = zeros(n,1); 
 im=100; maxits=200; tolIts=1.e-08;
%%-------------------------------------------------
%% MR - ssor preconditioning with omega = 1.

 PRE = struct('L',sparse(0,0),'U',sparse(0,0)) ; 
 [sol,res0] = pcg(A, PRE,'noPRE', rhs, x0,...
                       maxits,tolIts) ; 
 disp(' ** done noPre                  ** ') 
%%-------------------- ssor preconditioning 
 d = spdiags(A,0);
 E  = tril(A,-1);
 F  = triu(A,1);
%% 
 om = 1.0;
 it = 1;
 PRE.L  = spdiags(d,0,om*E);
 PRE.U  = spdiags(d,0,om*F);
 PRE.d2 = (2-om)*d;
 PRE.om = om;
 PRE.iter = it; 

 [sol,res1] = pcg(A,PRE,'precSSOR1', rhs, x0,...
                       maxits,tolIts) ; 
 %%
  disp(' ** done ssor 1                  ** ' )
 it = 1;
 om = 1.7;
 PRE.L  = spdiags(d,0,om*E);
 PRE.U  = spdiags(d,0,om*F);
 PRE.d2 = (2-om)*d;
 PRE.om = om;
 PRE.iter = it; 
 [sol,res2] = pcg(A,PRE,'precSSOR1', rhs, x0,...
                       maxits,tolIts) ; 
  disp(' ** done ssor 2                  ** ') 
 it = 3;
 om = 1.7;
 PRE.L  = spdiags(d,0,om*E);
 PRE.U  = spdiags(d,0,om*F);
 PRE.d2 = (2-om)*d;
 PRE.om = om;
 PRE.iter = it; 

 opts.type = 'crout';
 opts.droptol = 0.001;
 [L, U ] = ilu(A,opts);
 PRE.L = L;
 PRE.U = U;

 [sol,res3] = pcg(A,PRE,'precLU', rhs, x0,...
                  maxits,tolIts) ;
   disp(' ** done ssor 3                  ** ') 
 
%%-------------------- plot the 4 curves 
 its0 = [1:length(res0)];
 semilogy(its0,res0,'r--','LineWidth',3) 
 hold on 
 its1 = [1:length(res1)];
 semilogy(its1,res1,'b-','LineWidth',3) 
 its2 = [1:length(res2)];
 semilogy(its2,res2,'k-','LineWidth',3) 
 its3 = [1:length(res3)];
 semilogy(its3,res3,'b--','LineWidth',3) 
 h=legend('noPre','om=1,it=1','om=1.7,it=1 ',...
          'om=1.7,it=3','location','NorthEast') ;
 set(h,'fontsize',24)
