#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAXN 500

/* GLOBALNI PROMENNE */
double U1[MAXN][MAXN];
double U[MAXN][MAXN];
double b[MAXN][MAXN];

double h;
      

int    gs(int n, double eps);
void   uloz(const char *fname, int n);

main()
{
    int n, n2, num;
    int i, j;
    n = 10;
    h = 1. / n;
    n2 = (n - 1) * (n - 1);
    
    for(i = 0; i <= n; i++)
    {
        U[0][i] = 0; 
        U[n][i] = 0; 
        U[i][n] = 0; 
        U[i][0] = 0; 
    }

    for(i = 0; i <= n; i++)
        for(j = 0; j <= n; j++)
        {
            U[i][j] = 0;
            b[i][j] = h * h * sin(i * h * M_PI) * sin(j * h * M_PI);
        }
    
    num = gs(n, 1e-12);
    printf("GS %d x %d, eps 1e-12, iter %d\n", n2, n2, num);
    uloz("sol.dat", n);
}

int gs(int n, double eps)
{
    int i, j, iter;
    double rez, rezi; 
    iter = 0;
    rez = 1 + eps;
    while (eps < rez)
    {
        rez = 0;
        for(i = 1; i < n; i++)
            for(j = 1; j < n; j++)
            {
                rezi = (b[i][j] + U[i - 1][j] + U[i + 1][j]  + U[i][j -1] + U[i][j + 1] - 4 * U[i][j]) / 4.;
                U[i][j] += rezi;
                rezi = fabs(rezi);
                rez  = (rezi > rez) ? rezi : rez;
            }
        if (iter % n == 0)
            printf("iter %d, reziduum %g \n", iter, rez);
        iter++;
    }
    return iter;
}

void uloz(const char *fname, int n)
{
    FILE *fid;
    int i, j;

    printf("Reseni ukladam do souboru: %s \n", fname);
    fid = fopen(fname, "w");
    if (fid == NULL)
        exit(1);

    for(i = 0; i <= n; i++)
    {
        for(j = 0; j <= n; j++)
            fprintf(fid, "%20.10g\t", U[i][j]);
        fprintf(fid, "\n");
    }
    fclose(fid);
}