source: CIVL/examples/arithmetic/laplace.cvl@ 7de7353

int NX_BOUND = 2;
int NY_BOUND = 2;
int TIME_BOUND = 2;
int nx = 2;  // number of x coordinates (including boundary)
int ny = 2;  // number of rows including boundary
double epsilon = 0.01; // total error tolerance
$input double initialValues[ny][nx]; // initial values
$output int t;
$output double out[ny][nx];

double grid[ny][nx];    // holds values of current iteration

double square(double x) { return x * x; }

void init() {
        int row;
        int col;
        
        for (row=0; row<ny; row++)
                for (col=0; col<nx; col++)
                        grid[row][col] = initialValues[row][col];
}

void write_result(int time, double** data) {
        int row;
        int col;
        
        t = time;
        for (row=ny-1; row>=0; row--)
                for (col=0; col<nx; col++)
                        out[row][col] = data[row][col];
}

void main() {
        double error = epsilon;
        int row;
        int col;
        int time = 0;
        double result;
        double tmp[ny][nx];
        
        init();
        while (error >= epsilon && time < TIME_BOUND) {
                error = 0.0;
                for (row=1; row<ny-1; row++) {
                        for (col=1; col<nx-1; col++) {
                                tmp[row][col] = (grid[row-1][col]+grid[row+1][col]+
                                                                 grid[row][col-1]+grid[row][col+1])/4.0;
                                result = square(grid[row][col] - tmp[row][col]);
                                error += result;
                        }
                }
                for (row=1; row<ny-1; row++)
                        for (col=1; col<nx-1; col++)
                                grid[row][col] = tmp[row][col];
                time++;
        }
        write_result(time, grid);
}