/*
  Numerical Anaysis (E) 2016
  Tokyo Institute of Technology
 */

// To compile:
// gcc -o inverseiteration inverseiteration.c -lm 

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define	N		4 // matrix size
#define	SIGMA		1.75
//#define	SIGMA		3.21
//#define	SIGMA		5.67
#define	EPSILON		1.0E-4
#define	LOOP_COUNT	10

void setting(double (*)[N], double (*)[N], double (*));
void Gaussian_elimination(double (*)[N], double (*), double (*));
double inner_product(double (*), double (*));

int
main()
{
	int	loop, i;
	double	A[N][N], B[N][N], x[N], v[N], 
		mu, lambda, temp;

	setting(A, B, x);
	for(loop=0; loop<LOOP_COUNT; loop++) {
		Gaussian_elimination(B, x, v);
		mu = inner_product(x, v);
		lambda = SIGMA + 1.0/mu;
		printf(" This program is incomplete. You have to add code. See the source code for the details.\n");
		exit(1);
		// Add code here
		// to make the program work
		printf("%3d\t%10.8f\n", loop, lambda);
	}
	return 0;
}


void setting(double (*A)[N], double (*B)[N], double (*x))
{
	int i, j;

	printf("A=\n");
	printf("[\n");
	for(i=0; i<N; i++) {
		for(j=0; j<N; j++) {
			if (j==i) {
				A[i][j] = (double) (i+1);
			} else {
				A[i][j] = 0.0;
			}
			printf(" %.2f", A[i][j]);
		}
		printf("\n");
	}
	printf("]\n");

	for(i=0; i<N; i++) {
		for(j=0; j<N; j++) {
			if (j==i) {
				B[i][j] = A[i][j] - SIGMA;
			} else {
				B[i][j] = A[i][j];
			}
		}
	}

	for(i=0; i<N; i++) x[i] = 1.0/sqrt((double) N);
}

// solve Bx = b using the Gaussian elimination method
void Gaussian_elimination(double (*B)[N], double (*b), double (*x))
{
	int i, j, m, p;
	double r, temp, C[N][N];


	for(i=0; i<N; i++) {
		for(j=0; j<N; j++) {
			C[i][j] = B[i][j];
		}
	}

	for(m=0; m<N-1; m++) {
		p = m;
		for(i=m+1; i<N; i++) {
			if (fabs(C[p][m]) < fabs(C[i][m])) p = i;
		}

		if (fabs(C[p][m]) < EPSILON) {
			printf(" warning: C is a singular matrix \n");
		}

		if (p != m) {
			for(j=m; j<N; j++) {
				temp = C[m][j];
				C[m][j] = C[p][j];
				C[p][j] = temp;
			}
			
			temp = b[m];
			b[m] = b[p];
			b[p] = temp;
		}

		for(i=m+1; i<N; i++) {
			r = C[i][m]/C[m][m];
			for(j=m+1; j<N; j++) C[i][j] -= r*C[m][j];
			b[i] -= r*b[m];
		}
	}

	for(i=N-1; i>=0; i--) {
		if (i==N-1) {
			if (fabs(C[i][i]) < EPSILON) {
				printf(" warning: C is a singular matrix \n");
			}
			x[i] = b[i]/C[i][i];
		} else {
			temp = 0.0;
			for(j=i+1; j<N; j++) temp += C[i][j]*x[j];
			x[i] = (b[i] - temp)/C[i][i];
		}
	}
}


double inner_product(double (*a), double (*b)) 
{
	int		i;
	double	sum;

	sum = 0.0;
	for(i=0; i<N; i++) sum += a[i] * b[i];
	return(sum);
}