/*
 *  Main
 *  for LU Decomposition
 */

#include <Dialogs.h>
#include <Fonts.h>
#include <MacWindows.h>
#include <Menus.h>
#include <QuickDraw.h>
#include <TextEdit.h>

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

int LUdecomp(float **,int,int *,int *);
void LUsubstitute(float **,int,int *,float *);

#define DEBUG
/* #define DEBUG_LU */

#define SWAP
#define TEST_INVERSE
/* #define TEST_SOLVING */
#define DEBUG_SWAPS


enum
{
	kAlertID = 128
};

static void Initialize(void)
{
	InitGraf(&qd.thePort);
	InitFonts();
	InitWindows();
	InitMenus();
	TEInit();
	InitDialogs(nil);
	InitCursor();
}

void main(void)
{
	float	**A,**L,**U,**Ainv,**Aorig;
	float	*x;
	float	*b;
	int		dim;
	int		i,j,k;
	int		*rowswaps,val;
	float	temp;
	int		seed;
	
	/* Initialize(); */
	
	printf("this program generates a random matrix of nxn");
	printf(" input a random number seed: ");
	scanf("%d",&seed);
	
	srand(seed);

	printf("input dimension: ");
	scanf("%d",&dim);
	rowswaps = (int *)malloc(sizeof(int)*dim);
	A = (float **)malloc(sizeof(float *)*dim);
	for (i=0; i<dim; i++) {
		A[i] = (float *)malloc(sizeof(float)*dim);
	}
	Aorig = (float **)malloc(sizeof(float *)*dim);
	for (i=0; i<dim; i++) {
		Aorig[i] = (float *)malloc(sizeof(float)*dim);
	}
	
	#ifdef DEBUG_LU
		L = (float **)malloc(sizeof(float *)*dim);
		for (i=0; i<dim; i++) {
			L[i] = (float *)malloc(sizeof(float)*dim);
		}
		U = (float **)malloc(sizeof(float *)*dim);
		for (i=0; i<dim; i++) {
			U[i] = (float *)malloc(sizeof(float)*dim);
		}
	#endif
	
	b = (float *)malloc(sizeof(float)*dim);

	#ifdef TEST_SOLVING
		x = (float *)malloc(sizeof(float)*dim);
	#endif

	#ifdef TEST_INVERSE
		Ainv = (float **)malloc(sizeof(float *)*dim);
		for (i=0; i<dim; i++) {
			Ainv[i] = (float *)malloc(sizeof(float)*dim);
		}
	#endif

	
	/* fill A */
	for (i=0; i<dim; i++) {
		for (j=0; j<dim; j++) {
			Aorig[i][j] = A[i][j] = rand()/(float)RAND_MAX;
		}
	}
	
	#ifdef DEBUG
	/* print A */
	printf("given A\n");
	for (i=0; i<dim; i++) {
		for (j=0; j<dim; j++) {
			printf("%f ",A[i][j]);
		}
		printf("\n");
	}
	#endif
	
	#ifdef TEST_SOLVING
		/* fill x */
		for (i=0; i<dim; i++) {
			x[i] = rand()/(float)RAND_MAX;
		}

		#ifdef DEBUG
			/* print x */
			printf("given x\n");
			for (i=0; i<dim; i++) {
				printf("%f ",x[i]);
			}
			printf("\n");
		#endif
		
		/* compute b */
		for (i=0; i<dim; i++) {
			b[i]=0.0;
			for (j=0; j<dim; j++) {
				b[i] += A[i][j]*x[j];
			}
		}
	#endif
	
	LUdecomp(A,dim,rowswaps,&val);

	#ifdef DEBUG_LU
		/* print decomposition of A */
		printf("\nLU decomposition of A\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				printf("%f ",A[i][j]);
			}
			printf("\n");
		}
	#endif
		
	#ifdef DEBUG_SWAPS
		printf("rowswaps\n");
		for (i=0; i<dim; i++) printf("%d ",rowswaps[i]);
		printf("\n");
	#endif
	
	#ifdef DEBUG_LU
		/* extract L */
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++){
				if (i==j) L[i][j] = 1.0;
				if (i<j) L[i][j] = 0.0;
				if (i>j) L[i][j] = A[i][j];
			}
		}
		
		printf("\nL\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				printf("%f ",L[i][j]);
			}
			printf("\n");
		}
		
		/* extract U */
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++){
				if (i==j) U[i][j] = A[i][j];
				if (i<j) U[i][j] = A[i][j];
				if (i>j) U[i][j] =0.0;
			}
		}
		
		printf("\nU\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				printf("%f ",U[i][j]);
			}
			printf("\n");
		}
		
		/* reconstruct A */
		printf("\nreconstructed A\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				temp = 0.0;
				for (k=0; k<dim; k++) {
					temp += L[i][k]*U[k][j];
				}
				printf("%f ",temp);
			}
			printf("\n");
		}
	#endif

	#ifdef TEST_SOLVING
		/* print computed b */
		printf("computed b = Ax\n");
		for (i=0; i<dim; i++) {
			printf("%f ",b[i]);
		}
		printf("\n");
		
		
		LUsubstitute(A,dim,rowswaps,b);
		
		/* print computed x */
		printf("computed x\n");
		for (i=0; i<dim; i++) {
			printf("%f ",b[i]);
		}
		printf("\n");
	
	#endif
	
	#ifdef TEST_INVERSE
		printf("\n Testing Inverse \n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) b[j] = (i==j) ? 1:0;
			#ifdef DEBUG
				printf("column %d of identity = ");
				for (k=0; k<dim; k++) printf("%f ",b[k]);
				printf("\n");
			#endif
			LUsubstitute(A,dim,rowswaps,b);
			#ifdef DEBUG
				printf("computed column = ");
				for (k=0; k<dim; k++) printf(" %f",b[k]);
				printf("\n");
				printf("computed b[%d] = ",i);
				for (j=0; j<dim; j++) {
					temp = 0.0;
					for (k=0; k<dim; k++) temp += Aorig[j][k]*b[k];
					printf(" %f",temp);
				}
				printf("\n");
			#endif
			for (j=0; j<dim; j++) Ainv[j][i] = b[j];
		}
		
		/* print A-inverse */
		printf("A inverse\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				printf("%f ",Ainv[i][j]);
			}
			printf("\n");
		}
		
		/* print A times A-inverse */
		printf("A times A-inverse\n");
		for (i=0; i<dim; i++) {
			for (j=0; j<dim; j++) {
				temp = 0.0;
				for (k=0; k<dim; k++)  
					temp += Aorig[i][k]*Ainv[k][j];
				printf("%f  ",temp);
			}
			printf("\n");
		}
		
	#endif
	
	printf("done");

	/* NoteAlert(kAlertID, NULL); */
}