//Comment ca marche?
//	On fait un getf de ce fichier, et soit on fait tout à la main, soit on tape prepare()

//     Sans test du Epsilon (nombre d'iterations à effectuer)
// Gauss(A,x0,b) ou x0 est le vecteur de debut d'iterations   <- Algo. de Gauss-Seidel
// Jacobi(A,x0,b) ou x0 est le vecteur de debut d'iterations  <- Algo. de Jacobi

//     Avec le test du Epsilon
// Gauss2(A,x0,b) ou x0 est le vecteur de debut d'iterations  <- Algo. de Gauss-Seidel
// Jacobi2(A,x0,b) ou x0 est le vecteur de debut d'iterations <- Algo. de Jacobi

//genTriDiag(taille,a,b) genere une matrice carrée et tridiagonale ou
//	a est la valeur des elements de la diagonale de la matrice, et 
//	b celle des autres (qui touchent la diagonale) 








//************Fonction Gauss sans test du Epsilon***************/
function [L]=Gauss(A,x0,b)
L=list()

[n,m] = size(A);
//Verifie que cette matrice est une matrice carree
if n~=m | det(A)==0 then
	printf(" La matrice n''est pas carree, a un determinant nul");

else 
//verifie que ca va converger
	M=tril(A)
	N=M-A
	[valp,vecp]=bdiag(M^(-1)*N)      //Calcul du rayon spectral
	valeursp=diag(valp) 
	rayonspec=norm(abs(valeursp),'inf')
	if rayonspec > 1 
		printf(" La matrice A ne permet pas une convergence par cet algorithme");
	else
//Verifie que cette matrice est de bonne taille
		[p]=size(b);
		if p~=n then
			printf(" Second membre incompatible (taille incompatible) ");

//Verifie que cette matrice est de bonne taille
		else 
			[q]=size(x0);
			if q~=n then
				printf(" Vecteur de depart incompatible (taille incompatible) ");
			else
				nb = input("Combien d iterations : ");
				x=x0
				for d=1:nb
					x=M^(-1)*N*x+M^(-1)*b
					L(size(L)+1)=x
				end
			end	
		end
	end
end





//************Fonction Jacobi sans test du Epsilon***************/
function [L]=Jacobi(A,x0,b)
L=list()

[n,m] = size(A);
//Verifie que cette matrice est une matrice carree
if n~=m | det(A)==0 then
	printf(" La matrice n''est pas carree ou a un determinant nul");

else 
//verifie que ca va converger
	D=triu(A)-(A-tril(A))
	N=D-A
	[valp,vecp]=bdiag(D^(-1)*N);      //Calcul du rayon spectral
	valeursp=diag(valp) 
	rayonspec=norm(abs(valeursp),'inf')
	if rayonspec > 1 
		printf(" La matrice A ne permet pas une convergence par cet algorithme");
	else
//Verifie que cette matrice est de bonne taille
		[p]=size(b);
		if p~=n then
			printf(" Second membre incompatible (taille incompatible) ");

//Verifie que cette matrice est de bonne taille
		else 
			[q]=size(x0);
			if q~=n then
				printf(" Vecteur de depart incompatible (taille incompatible) ");
			else
				nb = input("Combien d iterations : ");
				x=x0
				for d=1:nb
					x=D^(-1)*N*x+D^(-1)*b
					L(size(L)+1)=x
				end
			end
		end
	end
end


//************Fonction Gauss avec test du Epsilon***************/	
function [L]=Gauss2(A,x0,b)
L=list()

[n,m] = size(A);
//Verifie que cette matrice est une matrice carree
if n~=m | det(A)==0 then
	printf(" La matrice n''est pas carree ou a un determinant nul");

else 
//verifie que ca va converger
	M=tril(A)
	N=M-A
	[valp,vecp]=bdiag(M^(-1)*N)      //Calcul du rayon spectral
	valeursp=diag(valp) 
	rayonspec=norm(abs(valeursp),'inf')
	if rayonspec > 1 
		printf(" La matrice A ne permet pas une convergence par cet algorithme");
	else
//Verifie que cette matrice est de bonne taille
		[p]=size(b);
		if p~=n then
			printf(" Second membre incompatible (taille incompatible) ");

//Verifie que cette matrice est de bonne taille
		else 
			[q]=size(x0);
			if q~=n then
				printf(" Vecteur de depart incompatible (taille incompatible) ");
			else
				epsilon = input("Quelle valeur de Epsilon: ");
				maxim = epsilon + 1;
   				while epsilon < maxim
    		 			x1 = M^(-1)*N*x0+M^(-1)*b
     					L(size(L)+1)=x1
     					for i = 1 : m
       				maxim = 0;
       					if abs(x0(i,1) - x1(i,1)) > maxim then
         						maxim = abs(x0(i,1) - x1(i,1));
       					end
     					end
     					x0 = x1;
				end
   			end             
		end
	end
end



//************Fonction Jacobi avec test du Epsilon***************/
function [L]=Jacobi2(A,x0,b)
L=list()

[n,m] = size(A);
//Verifie que cette matrice est une matrice carree
if n~=m | det(A)==0 then
	printf(" La matrice n''est pas carree ou a un determinant nul");
	
else 
//verifie que ca va converger
	D=triu(A)-(A-tril(A))
	N=D-A
	[valp,vecp]=bdiag(D^(-1)*N);      //Calcul du rayon spectral
	valeursp=diag(valp) 
	rayonspec=norm(abs(valeursp),'inf')
	if rayonspec > 1 
		printf(" La matrice A ne permet pas une convergence par cet algorithme");
	else
//Verifie que cette matrice est de bonne taille
		[p]=size(b);
		if p~=n then
			printf(" Second membre incompatible (taille incompatible) ");

//Verifie que cette matrice est de bonne taille
		else 
			[q]=size(x0);
			if q~=n then
				printf(" Vecteur de depart incompatible (taille incompatible) ");
			else
				epsilon = input("Quelle valeur de Epsilon: ");
				maxim = epsilon + 1;
   				while epsilon < maxim
    		 			x1 =  D^(-1)*N*x0+D^(-1)*b
     					L(size(L)+1)=x1
     					for i = 1 : m
       					maxim = 0;
       					if abs(x0(i,1) - x1(i,1)) > maxim then
         						maxim = abs(x0(i,1) - x1(i,1));
       					end
     					end
     					x0 = x1;
				end
   			end             
		end
	end
end




//************Fonction de creation d'une matrice tridiagonale***************/
function A=genTriDiag(taille,a,b)
	A=zeros(taille,taille)
	for i=1:taille	
		A(i,i)=a
		if i>=2
			A(i-1,i)=b
			A(i,i-1)=b
		end
	end




//************Fonction d'execution rapide***************/
function [L]=prepare()
printf('\nCreation et resolution d un systeme d equations Ax=b avec A matrice tridiagonale '); 

printf('\nCreation de la matrice A '); 
taille=input('Combien de variables (taille de la matrice): ')
diag1=input('Valeur de l unique veritable diagonale: ')
diag2=input('Valeur des diagonales adjacentes: ')
A=genTriDiag(taille,diag1,diag2)    

printf('\nCreation de la matrice b (dans l ordre)'); 
b=zeros(taille)
for i=1:taille
	b(i)=input('Resultat:');
end

printf('\nCreation de la matrice x0 (dans l ordre)'); 
deb=zeros(taille)
for i=1:taille
	deb(i)=input('Valeur initiale:');
end

	printf('\nChoix de l algorithme de resolution\n  1->Gauss sans Epsilon    \n  2->Jacobi sans Epsilon    \n  3->Gauss avec Epsilon    \n  4->Jacobi avec Epsilon')
 
	lequel=input('\nQuel algorithme?')
	printf('\nIterations: ');

		//Sans le test du Epsilon
	if lequel==1 then         
		L=Gauss(A,deb,b)         
	end
	if lequel==2 then         
		L=Jacobi(A,deb,b)         
	end
		//Avec le test du Epsilon
	if lequel==3 then         
		L=Gauss2(A,deb,b)         
	end
	if lequel==4 then         
		L=Jacobi2(A,deb,b)         
	end