380 lines
8.9 KiB
C
380 lines
8.9 KiB
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);
|
|
|
|
int readMatrixLine(FILE* file, double* matrixLine, int maxCols);
|
|
|
|
int main(int argc, char* argv[]) {
|
|
char filePath[MAX_FILE_PATH_LENGTH];
|
|
|
|
if(argc >= 2) {
|
|
strncpy(filePath, argv[1], MAX_FILE_PATH_LENGTH);
|
|
} else {
|
|
puts("Please enter the path of the file you'd like to open");
|
|
int result = scanf("%s" STR(MAX_FILE_PATH_LENGTH) "[^\n]", filePath);
|
|
if(result == EOF || result == 0 || filePath[0] == 0) {
|
|
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 == '0' || algorithm == '1') {
|
|
algorithm -= '0';
|
|
break;
|
|
}
|
|
if(algorithm == '\n' || algorithm == 0) {
|
|
puts("Defaulting to Jacobi");
|
|
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;
|
|
}
|
|
if(scan == 1) {
|
|
break;
|
|
}
|
|
flushStdin();
|
|
fputs("Invalid input - please enter a valid floating point number!\n", stderr);
|
|
}
|
|
|
|
Vector* result = solve(algorithm, matrix, b, x, e);
|
|
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) {
|
|
Vector* vectors = malloc(sizeof(Vector) * (MAX_ITERATION_STEPS + 1));
|
|
|
|
|
|
|
|
initVector(&vectors[0], b->n);
|
|
memcpy(vectors[0].data, x->data, b->n * sizeof(double));
|
|
|
|
int vectorCount = 1;
|
|
double temp;
|
|
double delta = 0.0;
|
|
|
|
if (method == 0){
|
|
|
|
do {
|
|
delta= 0.0;
|
|
initVector(&vectors[vectorCount],b->n);
|
|
for(int i=0; i<b->n; i++){
|
|
temp=0.0;
|
|
for(int j=0; j<b->n; j++){
|
|
if(j!=i){
|
|
temp= temp+A->data[i][j]*vectors[vectorCount-1].data[j];
|
|
}
|
|
}
|
|
vectors[vectorCount].data[i]=(b->data[i]-temp)/A->data[i][i];
|
|
delta=fmax(fabs(vectors[vectorCount].data[i]-vectors[vectorCount-1].data[i]),delta);
|
|
}
|
|
printVector(&vectors[vectorCount]);
|
|
vectorCount++;
|
|
} while (delta > e && vectorCount < MAX_ITERATION_STEPS);
|
|
|
|
}
|
|
|
|
if (method ==1){
|
|
|
|
|
|
do{
|
|
delta = 0.0;
|
|
if (vectorCount==1){
|
|
initVector(&vectors[vectorCount], b->n);
|
|
memcpy(vectors[vectorCount].data, x->data, b->n * sizeof(double));}
|
|
else {
|
|
initVector(&vectors[vectorCount], b->n);
|
|
memcpy(vectors[vectorCount].data, vectors[vectorCount-1].data, b->n * sizeof(double));}
|
|
|
|
for (int i = 0; i<b->n; i++){
|
|
temp=0.0;
|
|
for (int j=0; j<b->n; j++){
|
|
if (j!=i){
|
|
temp = temp + A->data[i][j]*vectors[vectorCount].data[j];
|
|
|
|
}
|
|
}
|
|
vectors[vectorCount].data[i] = (b->data[i]-temp)/A->data[i][i];
|
|
delta=fmax(fabs(vectors[vectorCount].data[i]-vectors[vectorCount-1].data[i]),delta);
|
|
|
|
}
|
|
printVector(&vectors[vectorCount]);
|
|
vectorCount++;
|
|
}while(delta > e && vectorCount < MAX_ITERATION_STEPS);
|
|
}
|
|
|
|
// MAX_ITERATION_STEPS enthält die maximal zulässige Anzahl an Iterationsschritten (100)
|
|
// Die einzelnen Vektoren müssen noch mit initVector initialisiert werden
|
|
|
|
|
|
// HIER kommt der Code hin ;)
|
|
|
|
// on success
|
|
// Sei x die Anzahl der durchgeführten Iterationschritte. Dann setzt vectors[x+1].n = 0. Damit weiß das folgende Programm wie viele Schritte getätigt wurden.
|
|
vectors[vectorCount].n=0;
|
|
|
|
return vectors;
|
|
}
|
|
|
|
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("%lf, ", 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("%lf, ", vector->data[i]);
|
|
}
|
|
puts("");
|
|
}
|