% graconju04.m
% fd, 10, 13 novembre 2006.
% fd, 28, 29 juin 2026
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   clear , format long
   
printf("graconju04 ; gradient conjugue preconditionne \n"); 
printf("             matrice penta-diagonale du laplacien a 2 dimensions\n"); 
printf("             preconditionnement par la matrice tridiagonale associee\n"); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   
n = 40 ; nbiter = n ; dx = 1 / n ; 
%                                            grille 
xx = zeros(1,n+1) ;  yy = zeros(1,n+1) ; 
for i = 1 : n+1 ;  xx(1,i:i) = (i-1) * dx ; end ;  
for j = 1 : n+1 ;  yy(1,j:j) = (j-1) * dx ; end ; 
%
uu=zeros(n+1,n+1) ;  bb=uu ; dx = 1 / n ; 
sgn = -1 ; %  sgn = 0 ; 
for i = 1 : n+1 ;  for j = 1 : n+1 ; 
bb(i:i,j:j) = 1 ; 
uu(i:i,j:j) = sgn ; sgn = -sgn ;  end ;  end ;  
for j = 1 : n+1 ; 
uu(1:1,j:j) = 0 ; uu(n+1:n+1,j:j) = 0 ;  
uu(j:j,1:1) = 0 ; uu(j:j,n+1:n+1) = 0 ;  end ; 
%
mesh ( xx , yy , uu ) ; 

gr = uu ;  for i = 2 : n ; for j = 2 : n ;
gr(i:i,j:j) = ( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2  - bb(i:i,j:j) ; 
end ; end ; 

ww = gr ;                                  %      preconditionnement 
aaaa = zeros(1,n-1) ; bbbb = aaaa ; cccc = aaaa ; 
for j=1:n-1 ; 
aaaa(1,j:j) = 4 ;  bbbb(1,j:j) = -1 ;   cccc(1,j:j) = -1 ; end ; 
sm = zeros(1,n-1) ; %  le second membre 
for j=2:n ; 
for i=1:n-1 ; sm(1,i:i) = gr(i+1:i+1,j:j)  ; end ;  
sol1d = zeros(1,n-1) ; sol1d=tridia(n-1,aaaa,bbbb,cccc,sm) ; 
for i=2:n ;  ww(i:i,j:j)  = sol1d (1,i-1:i-1) ; end ; 
end ;                    % mesh ( xx , yy , ww ) ; pause 
gg = 0 ; for i = 2 : n ; for j = 2 : n ;
gg = gg + ww(i:i,j:j) * gr(i:i,j:j) ; end ; end ; 

zco = zeros(1,nbiter) ;  
for iter = 1:nbiter  ; 
waw = 0 ; gw = 0 ;   for i = 2 : n ;  for j = 2 : n ; 
waw = waw +  ww(i:i,j:j) *  \
( 4*ww(i:i,j:j) - ww(i:i,j-1:j-1) - ww(i:i,j+1:j+1) \
      - ww(i-1:i-1,j:j) - ww(i+1:i+1,j:j) ) / dx**2 ; 
gw = gw + gr(i:i,j:j) * ww(i:i,j:j) ;   end ; end ;  
zco(1,iter:iter) = log(gw) / (2*log(10)) ;  
ro = - gw / waw ; 
for i = 2 : n ;   for j = 2 : n ; 
uu(i:i,j:j) = uu(i:i,j:j) + ro*ww(i:i,j:j) ; 
end  ; end ;

hold off 
mesh ( xx , yy , uu ) ; ; zlim([0 0.08])
pause (.1) 


%                               calcul de la fonction cout ; 
cout = 0 ;  for i = 2 : n ;  for j = 2 : n ;  
cout = cout + .5 * uu(i:i,j:j) * \
( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2 ; 
cout = cout -  bb(i:i,j:j) * uu(i:i,j:j) ; end ; end ; cout ; 
% 
grn = gr ;  for i = 2 : n ; for j = 2 : n ;
grn(i:i,j:j) = ( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2  - bb(i:i,j:j) ; 
end ; end ;  
     
smgr = grn ;   for j=2:n ;  %             preconditionnement 
for i=1:n-1 ; sm(1,i:i) = grn(i+1:i+1,j:j)  ; end ;  
sol1d = zeros(1,n-1) ; sol1d=tridia(n-1,aaaa,bbbb,cccc,sm) ; 
for i=2:n ;  smgr(i:i,j:j)  = sol1d (1,i-1:i-1) ; end ; 
end ;        %     mesh ( xx , yy , smgr ) ; pause 
ggn = 0 ; for i = 2 : n ;  for j = 2 : n ;  
ggn = ggn + smgr(i:i,j:j) * grn(i:i,j:j) ; end ; end ; 
alpha = ggn / gg ; 
gr = grn ; gg = ggn ; 
for i = 2 : n ;  for j = 2 : n ; 
ww(i:i,j:j) = smgr(i:i,j:j)  + alpha *  ww(i:i,j:j) ; end ;  end ;  

% mesh ( xx , yy , ww  ) ; pause
end ; 


hold off 
title ('solution finale de l algorithme de gradient conjugue preconditionne') ;
print('04-a-gradient-conjugue-preconditionne-matrice-penta-diagonale', '-djpg','-S891,630'); 
mesh ( xx , yy , uu ) ;
pause (1.)

plot ( zco  , '*b;logarithme de la norme du gradient;' ) ; 


title ('convergence de l algorithme') ;
print('04-b-gradient-conjugue-preconditionne-matrice-penta-diagonale', '-djpg','-S891,630'); 


