#include <math.h>
#include <iomanip.h>
#include <iostream.h>
#include <process.h>

void sor(double a[12][12], double b[12][12], double c[12][12],
		 double d[12][12], double e[12][12], double f[12][12],
		 double u[12][12], int jmax, double rjac)
{
	int n,j,l;
	double maxits,eps,zero,half,qtr,one,anormf,omega,aaa,bbb;
	double resid,anorm;
    maxits = 1000;
    eps = 0.00001;
    zero = 0.0;
    half = 0.5;
    qtr = 0.25;
    one = 1.0;
    anormf = zero;
    for (j = 2; j<=jmax - 1; j++)
        for (l = 2; l<=jmax - 1; l++)
            anormf = anormf + fabs(f[j][l]);
    omega = one;
    for (n = 1; n<=maxits; n++)
	{
		anorm = zero;
		for (j = 2; j<=jmax - 1; j++)
		{
			for (l = 2; l<=jmax - 1; l++)
			{
				if (((j + l) % 2) ==( n % 2))
				{
					aaa = a[j][l] * u[j + 1][l] + b[j][l] * u[j - 1][l];
					bbb = c[j][l] * u[j][l + 1] + d[j][l] * u[j][l - 1];
					resid = aaa + bbb + e[j][l] * u[j][l] - f[j][l];
					anorm = anorm + fabs(resid);
					u[j][l] = u[j][l] - omega * resid / e[j][l];
				}
			}
		}
        if (n == 1 )
            omega = one / (one - half * rjac * rjac);
        else
            omega = one / (one - qtr * rjac * rjac* omega);
        if ((n > 1) && (anorm < eps * anormf) ) return;
    }
    cout <<" maxits exceeded"<<endl;
}

void main()
{
    //program d16r1
    //driver for routine sor
	int jmax=11,i,j,midl;
	double pi= 3.1415926,rjac,aaa;
    double a[12][12],b[12][12],c[12][12],d[12][12];
    double e[12][12],f[12][12],u[12][12];
   
    for (i = 1; i<=jmax; i++)
	{
        for (j = 1; j<=jmax; j++)
		{
            a[i][j] = 1.0;
            b[i][j] = 1.0;
            c[i][j] = 1.0;
            d[i][j] = 1.0;
            e[i][j] = -4.0;
            f[i][j] = 0.0;
            u[i][j] = 0.0;
        }
    }
    midl = jmax / 2 + 1;
    f[midl][midl] = 2.0;
    rjac = cos(pi / jmax);
    sor(a,b,c,d,e,f,u,jmax,rjac);
    cout<<endl;
    cout <<"sor Solution:"<<endl;
    cout<<endl;
	cout<<setprecision(2)<<setiosflags(ios::fixed);
    for (i = 1; i<=jmax; i++)
	{
        for (j = 1; j<=jmax; j++)
            cout<<setw(7)<<u[i][j];
        cout<<endl;
	}
    cout<<endl;
    cout <<"Test that sulotion satisfies Difference Eqns:"<<endl;
    cout<<endl;
    for (i = 2; i<=jmax - 1; i++)
	{
		for (j = 2; j<=jmax; j++)
		{
			aaa = u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1];
			f[i][j] = aaa - 4.0 * u[i][j];
		}
		for (j = 2; j<=jmax-1 ; j++)
			cout <<setw(8)<<f[i][j];
		cout<<endl;
    }
}