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

#define N 4
#define MAX_ITER 150
#define EPS 1e-7

void print_system(double A[N][N], double b[N]);
void jacobi_solver(double A[N][N], double b[N]);
void gauss_seidel_solver(double A[N][N], double b[N]);
void verify_solution(double A[N][N], double b[N], double x[N]);

int main() {
// Preprocessed diagonally dominant system matrices
double A[N][N] = {
{ 5.0, -1.0, -2.0, 0.0},
{ 0.0, 8.0, 1.0, 7.0},
{-2.0, -7.0, -4.0, 0.0},
{ 4.0, 0.0, -3.0, -25.0}
};
double b[N] = {-1.32, 14.41, 0.70, -0.03};

print_system(A, b);

printf("\n======================================================\n");
printf(" 1. JACOBI METHOD \n");
printf("======================================================\n");
jacobi_solver(A, b);

printf("\n======================================================\n");
printf(" 2. GAUSS-SEIDEL METHOD \n");
printf("======================================================\n");
gauss_seidel_solver(A, b);

return 0;
}

void print_system(double A[N][N], double b[N]) {
printf("Preprocessed Diagonally Dominant Linear System:\n");
for (int i = 0; i < N; i++) {
printf("[ ");
for (int j = 0; j < N; j++) {
printf("%6.1f ", A[i][j]);
}
printf("] [x%d] = [ %6.2f ]\n", i + 1, b[i]);
}
}

void jacobi_solver(double A[N][N], double b[N]) {
double x[N] = {0.0, 0.0, 0.0, 0.0}; // Initial guess x^(0)
double x_new[N];
int converged = 0;

printf("Iter\tResidual SS\tx1\t\tx2\t\tx3\t\tx4\n");
printf("------------------------------------------------------------------------\n");

for (int k = 1; k <= MAX_ITER; k++) {
for (int i = 0; i < N; i++) {
double sum = 0.0;
for (int j = 0; j < N; j++) {
if (j != i) {
sum += A[i][j] * x[j];
}
}
x_new[i] = (b[i] - sum) / A[i][i];
}

// Calculate Residual Sum of Squares
double rss = 0.0;
for (int i = 0; i < N; i++) {
rss += pow(x_new[i] - x[i], 2);
}

// Display iteration row
printf("%d\t%.2e\t%.5f\t%.5f\t%.5f\t%.5f\n", k, rss, x_new[0], x_new[1], x_new[2], x_new[3]);

// Update values
for (int i = 0; i < N; i++) {
x[i] = x_new[i];
}

if (rss < EPS) {
printf("\n>> Converged successfully in %d iterations.\n", k);
converged = 1;
break;
}
}

if (!converged) {
printf("\nWarning: Did not reach convergence criteria.\n");
}

// --- Explicit Solution Display ---
printf("\n=== FINAL SOLUTION (JACOBI) ===\n");
for (int i = 0; i < N; i++) {
printf("x[%d] = %10.6f\n", i + 1, x[i]);
}

verify_solution(A, b, x);
}

void gauss_seidel_solver(double A[N][N], double b[N]) {
double x[N] = {0.0, 0.0, 0.0, 0.0}; // Initial guess x^(0)
double x_old[N];
int converged = 0;

printf("Iter\tResidual SS\tx1\t\tx2\t\tx3\t\tx4\n");
printf("------------------------------------------------------------------------\n");

for (int k = 1; k <= MAX_ITER; k++) {
for (int i = 0; i < N; i++) {
x_old[i] = x[i];
}

for (int i = 0; i < N; i++) {
double sum = 0.0;
for (int j = 0; j < N; j++) {
if (j != i) {
sum += A[i][j] * x[j];
}
}
x[i] = (b[i] - sum) / A[i][i];
}

// Calculate Residual Sum of Squares
double rss = 0.0;
for (int i = 0; i < N; i++) {
rss += pow(x[i] - x_old[i], 2);
}

// Display iteration row
printf("%d\t%.2e\t%.5f\t%.5f\t%.5f\t%.5f\n", k, rss, x[0], x[1], x[2], x[3]);

if (rss < EPS) {
printf("\n>> Converged successfully in %d iterations.\n", k);
converged = 1;
break;
}
}

if (!converged) {
printf("\nWarning: Did not reach convergence criteria.\n");
}

// --- Explicit Solution Display ---
printf("\n=== FINAL SOLUTION (GAUSS-SEIDEL) ===\n");
for (int i = 0; i < N; i++) {
printf("x[%d] = %10.6f\n", i + 1, x[i]);
}

verify_solution(A, b, x);
}

void verify_solution(double A[N][N], double b[N], double x[N]) {
printf("\n--- In-Code Verification Check (Ax == b) ---\n");
for (int i = 0; i < N; i++) {
double lhs_calc = 0.0;
for (int j = 0; j < N; j++) {
lhs_calc += A[i][j] * x[j];
}
double error = lhs_calc - b[i];
printf("Row %d: Computed Ax = %10.5f | Target b = %10.5f | Delta Error = %e\n",
i + 1, lhs_calc, b[i], error);
}
}