prg-1/main.c

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>

#define _STR(X) #X
#define STR(X) _STR(X)

#define MAX_FILE_PATH_LENGTH FILENAME_MAX

#define MAX_MATRIX_SIZE 500
#define MAX_MATRIX_IN_FILE_SIZE (MAX_MATRIX_SIZE + 2)

#define MAX_ITERATION_STEPS 100

typedef struct {
	int n;
	double** data;
} Matrix;
Matrix* createMatrix(void);
void freeMatrix(Matrix* matrix);
void createMatrixRows(Matrix* matrix, int rows);
void printMatrix(Matrix* matrix);

typedef struct {
	int n;
	double* data;
} Vector;
Vector* createVector(void);
void initVector(Vector* vector, int size);
void freeVector(Vector* vector);
void printVector(Vector* vector);

typedef enum {
	JACOBI = 0, GAUSS_SEIDEL = 1
} Method;

void flushStdin(void);

bool load(const char* filename, Matrix* A, Vector* b, Vector* x);
Vector* solve(Method method, Matrix* A, Vector* b, Vector* x, double e);
Vector* solveJacobi(Matrix* A, Vector* b, Vector* x, double e);
Vector* solveGaussSeidel(Matrix* A, Vector* b, Vector* x, double e);

int readMatrixLine(FILE* file, double* matrixLine, int maxCols);

int main(int argc, char* argv[]) {
	char filePath[MAX_FILE_PATH_LENGTH];

	if(argc >= 2) { // If arguments specified
		strncpy(filePath, argv[1], MAX_FILE_PATH_LENGTH); // use first argument as filePath
	} else { // Otherwise read from stdin
		puts("Please enter the path of the file you'd like to open");
		int result = scanf("%" STR(MAX_FILE_PATH_LENGTH) "[^\n]", filePath);
		if(result == EOF || result == 0 || filePath[0] == 0) { // If an error occured or no text has been entered -> exit
			return 0;
		}
		if(result != 1) {
			fputs("Couldn't read file path - stopping\n", stderr);
			return 1;
		}
		flushStdin();
	}

	Matrix* matrix = createMatrix();
	Vector* b = createVector();
	Vector* x = createVector();

	int returnCode = 0;

	while(true) {
		if(load(filePath, matrix, b, x)) {

			// Debug outputs
			puts("Data successfully loaded\nMatrix A:");
			printMatrix(matrix);
			puts("Vector b:");
			printVector(b);
			puts("Vector x:");
			printVector(x);
			
			puts("Please enter the algorithm to use:\n\t0: Jacobi (default)\n\t1: Gauss-Seidel");
			int algorithm;

			while(true) {
				algorithm = getchar();
				if(algorithm == EOF) {
					goto end;
				}
				if(algorithm == '\n' || algorithm == 0) {
					puts("Defaulting to Jacobi");
					algorithm = 0;
					break;
				}
				flushStdin();
				if(algorithm == '0' || algorithm == '1') {
					algorithm -= '0';
					break;
				}
				fputs("Please enter 0, 1 or leave empty to exit!\n", stderr);
			}

			puts("Please enter the precision to use:");
			double e;
			int scan;
			while(true) {
				scan = scanf("%lf", &e);
				if(scan == EOF) {
					goto end;
				}
				flushStdin();
				if(scan == 1) {
					break;
				}
				fputs("Invalid input - please enter a valid floating point number!\n", stderr);
			}

			Vector* result = solve(algorithm, matrix, b, x, e);
			if(result == NULL) {
				fputs("Failed to run the algorithm!\n", stderr);
				returnCode = 100;
			} else {
				int input;
				while(true) {
					puts("Choose what to output:\n\t0: Only the final result (default)\n\t1: All the iteration steps");
					input = getchar();
					if(input == '\n') {
						puts("Defaulting to final result");
						input = 0;
					} else {
						flushStdin();
					}
					if(input == 0 || input == '\n' || input == '0') {
						puts("Showing final result:");
						int i = 0;
						while(result[++i].n != 0 && i < MAX_ITERATION_STEPS)
							free(result[i-1].data);
						printVector(&result[i-1]);
						free(result[i-1].data);
						break;
					}
					if(input == '1') {
						puts("Showing all iterations:");
						for(int i = 0; i < MAX_ITERATION_STEPS; i++) {
							printVector(&result[i]);
							if(result[i].n == 0) {
								break;
							}
							free(result[i].data);
						}
						break;
					}
					fputs("Please enter 0, 1 or leave empty!", stderr);
				};
				free(result);
			}
			break;
		} else {
			fputs("Failed to load data from file.\nEnter new file path or leave empty to exit.\n", stderr);

			int result = scanf("%" STR(MAX_FILE_PATH_LENGTH) "[^\n]", filePath);
			if(result == EOF || result == 0 || filePath[0] == 0) {
				goto end;
			}
			if(result != 1) {
				fputs("Couldn't read file path - stopping\n", stderr);
				returnCode = 1;
				goto end;
			}
			flushStdin();
		}
	}

end:
	freeMatrix(matrix);
	freeVector(b);
	freeVector(x);

	printf("\nReturning with code %i\n", returnCode);
	return returnCode;
}

bool load(const char* filename, Matrix* matrix, Vector* b, Vector* x) {
	FILE* file = fopen(filename, "r");
	if(file == NULL) {
		fprintf(stderr, "Failed to open file \"%s\"\n", filename);
		return false;
	} else {
		double* firstLineBuffer = malloc(sizeof(double) * (MAX_MATRIX_IN_FILE_SIZE));

		int colsInFile = readMatrixLine(file, firstLineBuffer, MAX_MATRIX_IN_FILE_SIZE);

		if(colsInFile == -1) {
			fputs("Unexpected input on line 0", stderr);

			free(firstLineBuffer);
			goto failure;
		} else if(colsInFile == -2) {
			fprintf(stderr, "Exceeded maximum matrix size of %i\n", MAX_MATRIX_SIZE);

			free(firstLineBuffer);
			goto failure;
		} else {
			// success
			int cols = colsInFile - 1;
			initVector(b, cols);
			initVector(x, cols);

			b->data[0] = firstLineBuffer[cols];
			createMatrixRows(matrix, cols);
			matrix->data[0] = firstLineBuffer;

			int colsInLine;
			for(int i = 1; i < cols; i++) {
				matrix->data[i] = malloc(sizeof(double) * (colsInFile));
				colsInLine = readMatrixLine(file, matrix->data[i], colsInFile);
				if(colsInLine < 0) {
					if(i == cols - 1) {
						puts("Optional parameters are being used");

						// Matrix is one smaller than assumed
						cols--;
						b->n--;
						x->n--;
						matrix->n--;
						free(matrix->data[i]);

						x->data[0] = 0.4;

						// Copy b to x
						memcpy(x->data, b->data, b->n * sizeof(double));

						// extract b
						for(int j = 0; j < cols; j++) {
							b->data[j] = matrix->data[j][cols];
						}
						goto success;
					} else {
						fprintf(stderr, "Line %i contains illegal formatting - please fix!\n", i + 1);
						goto failure;
					}
				} else if(colsInLine != colsInFile) {
					fprintf(stderr, "Illegal line length found in line %i\n", i + 1);
					goto failure;
				} else {
					b->data[i] = matrix->data[i][cols];
				}
			}
			// successful but with no optional parameter

			// initialize vector x with zeros
			memset(x->data, 0, cols * sizeof(double));
			
		}
	}
success:
	fclose(file);
	return true;

failure:
	fclose(file);
	return false;
}


int readMatrixLine(FILE* file, double* matrixLine, int maxCols) {
	char nextChar = 0;
	int col = 0;
	double buffer;

	while(true) {
		if(fscanf(file, "%lf", &buffer) != 1) {
			return -1;
		}
		
		matrixLine[col] = buffer;

		col++;
		if(col > maxCols) {
			return -2;
		}
		nextChar = fgetc(file);
		if(nextChar == EOF || nextChar == '\n') break;
	}
	return col;
}

void flushStdin(void) {
	int c;
	while((c = getchar()) != '\n' && c != EOF && c != 0);
}

Vector* solve(Method method, Matrix* A, Vector* b, Vector* x, double e) {
	if(method == JACOBI) {
		return solveJacobi(A, b, x, e);
	} else {
		return solveGaussSeidel(A, b, x, e);
	}
}

Vector* solveJacobi(Matrix* A, Vector* b, Vector* x, double e) {
	Vector* history = malloc(sizeof(Vector) * (MAX_ITERATION_STEPS + 1));

	initVector(history, x->n);
	memcpy(history[0].data, b->data, b->n * sizeof(double));
	double delta;
	int i, j, steps = 1;
	double value;
	double* data;

	do {
		delta = 0.;

		initVector(&history[steps], x->n);
		data = history[steps].data;

		for(i = 0; i < x->n; i++) {
			value = b->data[i];
			for(j = 0; j < x->n; j++) {
				if(i != j)
					value -= history[steps-1].data[i] * A->data[i][j];
			}
			data[i] = value / A->data[i][i];

			delta = fmax(delta, fabs(data[i] - history[steps - 1].data[i]));
		}

		steps++;
	} while(delta > e && steps < MAX_ITERATION_STEPS);

	if(steps < MAX_ITERATION_STEPS) {
		history[steps].n = 0;
	}

	return history;
}

Vector* solveGaussSeidel(Matrix* A, Vector* b, Vector* x, double e) {
	Vector* history = malloc(sizeof(Vector) * (MAX_ITERATION_STEPS + 1)); 

	initVector(history, x->n);
	memcpy(history[0].data, b->data, x->n * sizeof(double));

	double delta;
	int i, k, steps = 1;
	double* data;
	double value;

	do {
		delta = 0.;

		initVector(&history[steps], x->n);
		data = history[steps].data;

		for(k = 0; k < x->n; k++) {
			value = b->data[k];
			for(i = 0; i < k; i++) {
				value -= A->data[k][i] * data[i];
			}
			for(i = k + 1; i < x->n; i++) {
				value -= A->data[k][i] * history[steps - 1].data[i];
			}
			data[k] = value / A->data[k][k];

			delta = fmax(delta, fabs(data[k] - history[steps - 1].data[k]));
		}

		steps++;
	} while(delta > e && steps < MAX_ITERATION_STEPS);

	if(steps < MAX_ITERATION_STEPS) {
		history[steps].n = 0;
	}

	return history;
}

inline Matrix* createMatrix(void) {
	Matrix* matrix = malloc(sizeof(Matrix));
	matrix->n = 0;
	return matrix;
}

inline void createMatrixRows(Matrix* matrix, int rows) {
	matrix->n = rows;
	matrix->data = malloc(sizeof(double*) * rows);
}

inline void freeMatrix(Matrix* matrix) {
	for(int i = 0; i < matrix->n; i++) {
		free(matrix->data[i]);
	}
	free(matrix->data);
	free(matrix);
}

void printMatrix(Matrix* matrix) {
	for(int i = 0; i < matrix->n; i++) {
		for(int j = 0; j < matrix->n; j++) {
			printf("%le, ", matrix->data[i][j]);
		}
		puts("");
	}
}

inline Vector* createVector(void) {
	Vector* vector = malloc(sizeof(Vector));
	vector->n = 0;
	return vector;
}

inline void initVector(Vector* vector, int size) {
	vector->n = size;
	vector->data = malloc(sizeof(double) * size);
}

inline void freeVector(Vector* vector) {
	free(vector->data);
	free(vector);
}

void printVector(Vector* vector) {
	for(int i = 0; i < vector->n; i++) {
		printf("%le, ",  vector->data[i]);
	}
	puts("");
}