/* 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&);
}


// my code appears to generate the correct orders of convergence.

int explicitCenteredDiffs(int size)
{


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



        // 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;
		}
	}
	d[n-1] = 1;

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

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

	double pi2 = PI*PI*.25*h*h;
	for(int j = 0; j < n; j++){
	  forcing[j] = pi2*sin(PI*.5*h*(j+1));
	  //forcing[j] = 0;
	}
	forcing[n-1] += h;

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

	double solution[n];
	for(i = 0; i < n; i++){
		solution[i] = sin(PI*.5*h*(i+1)) + h*(i+1);
	}


        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]);
        }
        fclose(finalState);


	filename = "./meanSquaredError.dat";
	FILE * finalState2 = fopen(filename.c_str(),"a+");

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

	fclose(finalState2);

	filename = "./soln1.dat";
	FILE * finalState3 = fopen(filename.c_str(),"a+");

	// compute error at x = 1
	
	error = fabs(forcing[n-1] - solution[n-1]);

	fprintf(finalState3,"%f %f\n",-log(h),log(error));

	fclose(finalState3);

	filename = "./deriv0.dat";
	FILE * finalState4 = fopen(filename.c_str(),"a+");

	// compute deriv error at x = 0
	
	// exact deriv(0) = PI/2
	// forcing[0] is our first data point, so (forcing[0]-0)/h is ~ the deriv at x=0
	error = fabs(forcing[0]/h - (PI*.5+1));

	fprintf(finalState4,"%f %f\n",-log(h),log(error));

	fclose(finalState4);

          return 0;
} 


int main (int argc,char* argv[])
{
	
	string filename = "./meanSquaredError.dat";
	FILE * finalState2 = 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);
	
	}


	
	return 0;
}