/* Daniel Summerhays, 9-14-06 */

#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <math.h>

using namespace std;

double PI = 3.1415926535897932;

extern "C" {
  void dgttrf_(const int&,double*,double*,double*,double*,int*,int&);
  void dgttrs_(const char&,const int&,const int&,double*,double*,double*,double*,int*,double*,int&,int&);
}


int explicitCenteredDiffs(int size)
{


        int n = size-2;
	double h = 1./(size-1.);



        // data vector

	int order = n;
	double dl[n-1];
	double d[n];
	double du[n-1];
	double du2[n-2];
	int ipiv[n];
	int info;

	int i;
	for(i = 0; i < n; i++){
		d[i] = 2;
		if(i < n-1){
			dl[i] = -1;
			du[i] = -1;
		}
		if(i < n-2){
			du2[i] = 0;
		}
	}

	dgttrf_(order,dl,d,du,du2,ipiv,info);

	char trans = 'N';
	int nrhs = 1;
	double forcing[n]; // B in dgttrs

	if(n%2 == 1){ // odd number of elements
		for(int j = 0; j < n; j++){
			forcing[j] = 0;
		}
		forcing[n/2] = h;
	}
	if(n%2 == 0){ // even number of elements
		for(int j = 0; j < n; j++){
			forcing[j] = 0;
		}
		forcing[n/2-1] = h/2;
		forcing[n/2] = h/2;
	}

// 	DGTTRS_(&trans,&order,&nrhs,dl,d,du,du2,ipiv,forcing,&n,&info);
	dgttrs_(trans,order,nrhs,dl,d,du,du2,ipiv,forcing,n,info);


        cout << "PRINTING DATA:::::" << endl << endl;
        //char* filename = "./finalState.dat";
	ostringstream oss;

	string filename;
	oss << "./finalStateSize" << size << ".dat";

	filename = oss.str();


        FILE * finalState = fopen(filename.c_str(),"w");

	fprintf(finalState,"0.0 0.0\n");
        for(int k = 0; k < n; k++)
        {
		fprintf(finalState,"%f %f\n",h*(k+1),forcing[k]);
        }
	fprintf(finalState,"1.0 0.0\n");
        fclose(finalState);

	if(n%2 == 0)
		filename = "./meanErrorEven.dat";
	else
		filename = "./meanErrorOdd.dat";
	FILE * finalState2 = fopen(filename.c_str(),"a+");

	// compute mean squared error
	double error = 0;
	for(i = 0; i < n; i++){
		error += forcing[i]*h;
	}
	
	error -= .125;
	error = fabs(error);
	
	fprintf(finalState2,"%f %f\n",log(size-1),log(error));

	fclose(finalState2);

          return 0;
} 


int main (int argc,char* argv[])
{
	// wiping out the data files.
	string filename = "./meanErrorEven.dat";
	FILE * finalState2 = fopen(filename.c_str(),"w");
	filename = "./meanErrorOdd.dat";
	FILE * finalState5 = fopen(filename.c_str(),"w");
	filename = "./soln1.dat";
	FILE * finalState3 = fopen(filename.c_str(),"w");
	filename = "./deriv0.dat";
	FILE * finalState4 = fopen(filename.c_str(),"w");
	
	fclose(finalState2);
	fclose(finalState3);
	fclose(finalState4);
	
        
    int size;
    

	for(size = 10; size < 320; size*=2){
        
		explicitCenteredDiffs(size);
		explicitCenteredDiffs(size+1);
	
	}


	
	return 0;
}