#include <stdio.h>

void tridiaglu(int n,double a[n],
	double b[n-1],double c[n-1]){
	int k;
	for(k=1;k<n;k++){
		b[k-1]=b[k-1]/a[k-1];
		a[k]=a[k]-b[k-1]*c[k-1];
	}
}
void lusolve(int n,double a[n],
	double b[n-1],double c[n-1],
	double f[n],double y[n]){
	double z[n];
	z[0]=f[0];
	int k;
	for(k=1;k<n;k++) z[k]=f[k]-b[k-1]*z[k-1];
	y[n-1]=z[n-1]/a[n-1];
	for(k=n-2;k>=0;k--) y[k]=(z[k]-c[k]*y[k+1])/a[k];
}
void printvec(int n,double x[n]){
	int k;
	for(k=0;k<n;k++) printf("%g ",x[k]);
	printf("\n");
}
int main(){
	double b[4]={ 1, 1, 1, 1 };
	double a[5]={ 4, 4, 4, 4, 4 };
	double c[4]={ 2, 2, 2, 2 };
	double f[5]={-1, 2, 2, 2, -1};
	double y[5];
	printf("input b="); printvec(4,b);
	printf("input a="); printvec(5,a);
	printf("input c="); printvec(4,c);
	printf("input f="); printvec(5,f);
	tridiaglu(5,a,b,c);
	printf("b="); printvec(4,b);
	printf("a="); printvec(5,a);
	printf("c="); printvec(4,c);
	lusolve(5,a,b,c,f,y);
	printf("solution to Ay=f where A is tridiagonal\n"
		"with input a on the diagonal, input b lower\n"
		"diagonal and input c upper diagonal\n");
	printf("y="); printvec(5,y);
	return 0;
}