/* To Solve a differential equation 
 *
 *    2
 *   d y
 *   --- = y + x
 *     2
 *   dx 
 *
 */
// Includes
#include <stdio.h>
#include <math.h>

// Defines
#define		MAX		100

// External Routines From Lapack
void sgesv_(long* n, long* nrhs, float* a, long* lda,
			long* ipiv, float* b, long* ldb, long* info);

float exactSolution(float x)
{
	float a = 0.8509181282;
	return -x + a * ( exp(x) - exp(-x) );
}

void printMat(int n, int m, float* a, int lda)
{
	int i, j;
	for  ( i = 0; i < n; i++ )
	{
		for ( j = 0; j < m; j++ )
			printf("%7.4f\t",a[i*lda+j]);
		printf("\n");
	}
}

processArguments(int argc, char* argv[], long* n, long* debuglevel)
{
	switch(argc) {
		case 3:
			*debuglevel = atoi(argv[2]);
			*n = atoi(argv[1]);
			return;
		case 2:
			*debuglevel = 0;
			*n = atoi(argv[1]);
			return;
		default:
			*n			= 10;
			*debuglevel	= 0;
	}
	return;
}

int main(int argc, char* argv[])
{
	long	n;
	long	nrhs;
	float	a[MAX][MAX];
	long	lda = MAX;
	long	ipiv[MAX];
	float	b[MAX][MAX];
	long	ldb;
	long	info;
	long	debuglevel;
	float	h, h2;
	float	x;
	int		i, j;

	processArguments(argc, argv, &n, &debuglevel);

	if ( n > MAX ) return -1;

	// Setup A matrix
	h = 1.0 / (n+1);
	h2 = h * h;
	for ( i = 0; i < n; i++ ) {
		for ( j = i; j < n; j++ ) {
			if ( i == j )
				a[i][i] = - 2.0 * ( 1.0 + h2 );
			else if ( j == i + 1 ) {
				a[i][j] = 1.0;
				a[j][i] = 1.0;
			}
			else {
				a[i][j] = 0.0;
				a[j][i] = 0.0;
			}
		}
	}
	if ( debuglevel ) printMat(n,n,(float*)a,lda/*,"A matrix"*/);

	// Setup B matrix;
	for ( i = 0; i < n ; i++ )
		b[0][i] = 2 * (i + 1) * h2 * h;
	b[0][n-1] -= 1.0;		// Boundary condition at x = 1;
	if ( debuglevel ) printMat(n, 1, (float*)b, 1/*, "B matrix"*/);

	// Setup other data
	nrhs = 1;
	lda = MAX;
	ldb = MAX;

	// Call the lapack subroutine
	sgesv_(&n, &nrhs, (float*)a, &lda, ipiv, (float*)b, &ldb, &info);

	// If everything is all-right, print out the result
	if ( info == 0 ) {
		printf("Solution\n");
		printf("x\texact\t\tnumerical\n");
		printf("0.000\t%10.6f\t%10.6f\n",0,0);
		for ( i = 0; i < n; i++ ) {
			x = (i + 1) * h;
			printf("%5.3f\t%10.6f\t%10.6f\n",x,exactSolution(x),b[0][i]);
		}
		printf("1.000\t%10.6f\t%10.6f\n",1.0,1.0);
	}
}