clear


//////////////////////////
//
//
//  fonction preconditionnement SSOR: 
//  en entree: n\in N*, r\in R^n , A = D-E-F \in M_n, w \in ]0,2[ 
//
//  On suppose Aii non nuls pour tous i=1,...,n
//
//  en sortie p tel que (D/w - E)D-1(D/w - F)p = r 
//
function p=SSOR(n,A,w,r)
    
    // descente (D/w-E)p1 = r
    
    p1 = zeros(n,1)
    p= zeros(n,1)
    
    
    p1(1,1) = w*r(1,1)/A(1,1);
    for i=2:n
        p1(i,1) = w*(r(i,1) - A(i,1:i-1)*p1(1:i-1,1))/A(i,i);
    end
    
    // p1 <- D*p1 
    
    for i=1:n
        p1(i,1) = A(i,i)*p1(i,1)
    end
    
    // remontee (D/w-F)p = p1
    
    p(n,1) = w*p1(n,1)/A(n,n); 
    for i=n-1:-1:1
       p(i,1) = w*(p1(i,1)- A(i,i+1:n)*p(i+1:n,1))/A(i,i); 
    end
    
endfunction



//////////////////////////

// DONNEES 

// EXEMPLE DU LAPLACIEN 2D 

nx = 20;
ny = nx; 

n = nx*ny;

dx = 1/(nx+1);
dy = 1/(ny+1);

//A=sparse([],[],[n,n]); // matrice creuse

u = (2/dx**2 + 2/dy**2)*ones(n,1)
ux = -1/dx**2*ones(n-1,1)
uy = -1/dy**2*ones(n-nx,1) 

ux1 = -1./dx**2*ones(n-1,1)
uy1 = -1/dy**2*ones(n-nx,1) 

A=diag(u) + diag(ux1,1) + diag(ux,-1) + diag(uy1,nx) + diag(uy,-nx)



for k=1:n
    xbar(k,1) = rand(k)
end

b = A*xbar;

///////////////////////////




// Richardson a pas variable

eps = 1.0E-10

w = 1.5


x = zeros(n,1);
r = b-A*x;
nr0 = norm(r);
nr = nr0;
k = 0;


while (nr/nr0>eps)&(k<1000) then 
        
    k = k + 1;    
    res(k) = log(nr/nr0)/log(10);
    
    p = SSOR(n,A,w,r); 
//    p = r;    
    q = A*p;
    alpha = (p'*r)/(p'*q);
    x = x + alpha*p;
    r = r - alpha*q;
    nr = norm(r);   
        
end


plot(res,'-xr')









