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

// Defines
#define		MAX		100
#define		PI2		9.86960440

// External Routines From Lapack
void ssyev_(char* jobz, char* uplo, long* n, float* a, long* lda,
			float* w, float* work, long* lwork, long* info);

float exactSolution(int i)
{
	return PI2 * i * i;
}

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[])
{
	char	jobz = 'N';
	char	uplo = 'U';
	long	n;
	float	a[MAX][MAX];
	long	lda = MAX;
	float	w[MAX];
	long	lwork = 3*MAX;
	float	work[3*MAX];
	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;
			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 other data
	lda = MAX;

	// Call the lapack subroutine
	ssyev_(&jobz, &uplo, &n, (float*)a, &lda, w, work, &lwork, &info);

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