/* prototypes of standard functions */
#include <math.h>
#include <stdio.h>
#include <string.h>


/* DPTSV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) */
void dptsv_(int *N, int *NRHS, double *D, double *E, double *B , int *LDB, int *INFO);
void dpttrf_ (int *N, double *D, double *E, int *Info);
void dpttrs_ (int *N, int *NRHS, double *D, double *E, double *B, int *LDB, int *Info); 

void  aufgabe1_neu()
{

	double D[500];
	double Db[500];
	double E[499];
	double Eb[500];
	double x0[500];
	double f1;
	double f2;
	double f3;

	int i;
	const int s=51000;
	double h=1.0/s;
#ifdef FULLOUTPUT
	int j;
#endif
	
	for (i=0;i<499;i++){
		D[i]=1+502002*h;
		E[i]=-251001*h;
		x0[i]=1;
	}
	D[499]=1+h*502002;
	x0[499]=1;
  	int ldA=500, n=500;

  int one=1;
  int info;
  //Matrix LU zerlegen
  dpttrf_(&n, D, E, &info);


  /* der Funktionsaufruf zum loesen des LGS */
  for (i = 0; i < s; i++){
//	memcpy(Db, D, 500*sizeof(double));
//	memcpy(Eb, E, 499*sizeof(double));
//	dptsv_(&n,  &one, Db,Eb,x0,&ldA,&info);
	dpttrs_(&n, &one, D, E, x0, &ldA, &info);
	if ( fabs((i+1)*h-0.1)<0.0000000001 ){
		printf("t = 0.1\t n = %d\n", i+1);
		f1=fabs(0.002976799850281-x0[0])/0.002976799850281;
		f2=fabs(0.410969184652936-x0[166])/0.410969184652936;
		f3=fabs(0.002976799850260-x0[499])/0.002976799850260;
		printf("x(1)=\t%11.4e\trel.Fehler: %g\n",x0[0], f1);
		printf("x(167)=\t%11.4e\trel.Fehler: %g\n",x0[166], f2);
		printf("x(500)=\t%11.4e\trel.Fehler: %g\n",x0[499], f3);
	}
#ifdef FULLOUTPUT
		for ( j = 0; j < 500; j++) printf("%11.4e\t",x0[j]);
		printf("\n");
#endif
  }
	f1=fabs(0.000000412968866-x0[0])/0.000000412968866;
	f2=fabs(0.000057034631830-x0[166])/0.000057034631830;
	f3=fabs(0.000000412968866-x0[499])/0.000000412968866;
	printf("\nt = 1.0\n");
	printf("x(1)=\t%11.4e\trel.Fehler: %g\n",x0[0], f1);
	printf("x(167)=\t%11.4e\trel.Fehler: %g\n",x0[166], f2);
	printf("x(500)=\t%11.4e\trel.Fehler: %g\n",x0[499], f3);

}