Jacobi Relaxation.cs

    internal class Program
    {
        private const int size = 500;
        private const int iterations = 5000;

        private static void Main( string[] args )
        {
            TimeSpan elapsed = TimeOriginal();
            Console.WriteLine( elapsed );

            elapsed = TimeModified();
            Console.WriteLine( elapsed );

            elapsed = TimeOptimized();
            Console.WriteLine(elapsed);

            Console.ReadLine();
        }

        private static TimeSpan TimeOptimized()
        {
            var oldGrid = new double[size][];
            var grid = new double[size][];
            for (int i = 0; i < size; i++)
            {
                oldGrid[i] = new double[size];
                grid[i] = new double[size];
            }
            Stopwatch stopwatch = Stopwatch.StartNew();
            Jacobi3(oldGrid, grid);
            stopwatch.Stop();
            return stopwatch.Elapsed;
        }

        private static TimeSpan TimeModified()
        {
            var oldGrid = new double[size][];
            var grid = new double[size][];
            for ( int i = 0; i < size; i++ )
            {
                oldGrid[i] = new double[size];
                grid[i] = new double[size];
            }
            Stopwatch stopwatch = Stopwatch.StartNew();
            Jacobi2( oldGrid, grid );
            stopwatch.Stop();
            return stopwatch.Elapsed;
        }

        private static TimeSpan TimeOriginal()
        {
            var oldGridD = new double[size,size];
            var gridD = new double[size,size];
            Stopwatch stopwatch = Stopwatch.StartNew();
            Jacobi( oldGridD, gridD );
            stopwatch.Stop();
            TimeSpan elaped = stopwatch.Elapsed;
            return elaped;
        }

        private static double[,] Jacobi( double[,] oldGrid, double[,] grid )
        {
            for ( int i = 0; i < iterations; i++ )
            {
                for ( int x = 0; x < size; x++ )
                {
                    for ( int y = 0; y < size; y++ )
                    {
                        int denom = 1;
                        double sum = oldGrid[x, y];
                        //left neighbor 
                        if ( x > 0 )
                        {
                            sum += oldGrid[x - 1, y];
                            denom++;
                        }
                        //right neighbor 
                        if ( x < size - 1 )
                        {
                            sum += oldGrid[x + 1, y];
                            denom++;
                        }
                        //upper neighbor 
                        if ( y > 0 )
                        {
                            sum += oldGrid[x, y - 1];
                            denom++;
                        }
                        //lower neighbor 
                        if ( y < size - 1 )
                        {
                            sum += oldGrid[x, y + 1];
                            denom++;
                        }
                        //left-upper neighbor 
                        if ( ( x > 0 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x - 1, y - 1];
                            denom++;
                        }
                        //right-upper neighbor 
                        if ( ( x < size - 1 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x + 1, y - 1];
                            denom++;
                        }
                        //right-lower neighbor 
                        if ( ( x < size - 1 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x + 1, y + 1];
                            denom++;
                        }
                        //left-lower neighbor 
                        if ( ( x > 0 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x - 1, y + 1];
                            denom++;
                        }

                        grid[x, y] = sum / denom;
                    }
                }

                grid[size / 2, size / 2] = 1; //maintain the center point's value

                //swap the pointers of the first and second array; this will ensure that the next 
                //iteration will rely on the values calculated in this iteration 
                double[,] temp = grid;
                grid = oldGrid;
                oldGrid = temp;
            }
            return grid;
        }

        private static double[][] Jacobi2( double[][] oldGrid, double[][] grid )
        {
            for ( int i = 0; i < iterations; i++ )
            {
                for ( int x = 0; x < size; x++ )
                {
                    for ( int y = 0; y < size; y++ )
                    {
                        int denom = 1;
                        double sum = oldGrid[x][y];
                        //left neighbor 
                        if ( x > 0 )
                        {
                            sum += oldGrid[x - 1][y];
                            denom++;
                        }
                        //right neighbor 
                        if ( x < size - 1 )
                        {
                            sum += oldGrid[x + 1][y];
                            denom++;
                        }
                        //upper neighbor 
                        if ( y > 0 )
                        {
                            sum += oldGrid[x][y - 1];
                            denom++;
                        }
                        //lower neighbor 
                        if ( y < size - 1 )
                        {
                            sum += oldGrid[x][y + 1];
                            denom++;
                        }
                        //left-upper neighbor 
                        if ( ( x > 0 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x - 1][y - 1];
                            denom++;
                        }
                        //right-upper neighbor 
                        if ( ( x < size - 1 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x + 1][y - 1];
                            denom++;
                        }
                        //right-lower neighbor 
                        if ( ( x < size - 1 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x + 1][y + 1];
                            denom++;
                        }
                        //left-lower neighbor 
                        if ( ( x > 0 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x - 1][y + 1];
                            denom++;
                        }

                        grid[x][y] = sum / denom;
                    }
                }

                grid[size / 2][size / 2] = 1; //maintain the center point's value

                //swap the pointers of the first and second array; this will ensure that the next 
                //iteration will rely on the values calculated in this iteration 
                double[][] temp = grid;
                grid = oldGrid;
                oldGrid = temp;
            }
            return grid;
        }

        private static double[][] Jacobi3( double[][] oldGrid, double[][] grid )
        {
            int i = 0;
            do
            {
                int x = 0;
                do
                {
                    int y = 0;
                    do
                    {
                        int denom = 1;
                        double sum = oldGrid[x][y];
                        //left neighbor 
                        if ( x > 0 )
                        {
                            sum += oldGrid[x - 1][y];
                            denom++;
                        }
                        //right neighbor 
                        if ( x < size - 1 )
                        {
                            sum += oldGrid[x + 1][y];
                            denom++;
                        }
                        //upper neighbor 
                        if ( y > 0 )
                        {
                            sum += oldGrid[x][y - 1];
                            denom++;
                        }
                        //lower neighbor 
                        if ( y < size - 1 )
                        {
                            sum += oldGrid[x][y + 1];
                            denom++;
                        }
                        //left-upper neighbor 
                        if ( ( x > 0 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x - 1][y - 1];
                            denom++;
                        }
                        //right-upper neighbor 
                        if ( ( x < size - 1 ) && ( y > 0 ) )
                        {
                            sum += oldGrid[x + 1][y - 1];
                            denom++;
                        }
                        //right-lower neighbor 
                        if ( ( x < size - 1 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x + 1][y + 1];
                            denom++;
                        }
                        //left-lower neighbor 
                        if ( ( x > 0 ) && ( y < size - 1 ) )
                        {
                            sum += oldGrid[x - 1][y + 1];
                            denom++;
                        }

                        grid[x][y] = sum / denom;
                        ++y;
                    } while ( y < size );
                    ++x;
                } while ( x < size );

                grid[size >> 2][size >> 2] = 1; //maintain the center point's value

                //swap the pointers of the first and second array; this will ensure that the next 
                //iteration will rely on the values calculated in this iteration 
                double[][] temp = grid;
                grid = oldGrid;
                oldGrid = temp;

                ++i;
            } while ( i < iterations );
            return grid;
        }
    }