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

void tridag(double a[101],double b[101],double c[101],double r[101],double u[101],int n)
{
	int nmax,j;
	double bet,gam[101];
    nmax = 100; 
    if (b[1] == 0.0) _c_exit();
    bet = b[1];
    u[1] = r[1] / bet;
    for (j = 2; j<= n; j++)
	{
        gam[j] = c[j - 1] / bet;
        bet = b[j] - a[j] * gam[j];
        if (bet ==0.0) _c_exit();
        u[j] = (r[j] - a[j] * u[j - 1]) / bet;
    }
    for (j = n - 1 ; j>=1; j--)
        u[j] = u[j] - gam[j + 1] * u[j + 1];
}

void adi(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 g[12][12],double u[12][12],int jmax,int k,
		 double alpha,double beta,double eps)
{
	int next1,jj,kk,k1,j,n,l,kits,maxits;
	double nrr,zero,two,half,nr,ab,disc,anormg,resid,nits,anrom,anorm;  
    double aa[101],bb[101],cc[101],rr[101],uu[101];
    double psi[101][101],alph[7],bet[7],r[33],s[33][7],rfact;
	jj = 100;
    kk = 6;
    nrr = pow(2 , (kk - 1));
    maxits = 100;
    zero = 0.0;
    two = 2.0;
    half = 0.5;
    if (jmax > jj ) cout<< " 'increase jj'"<<endl;
    if( k > kk - 1) cout<< " 'increase kk'"<<endl;
    k1 = k + 1;
    nr = pow(2 , k);
    alph[1] = alpha;
    bet[1] = beta;
    for (j = 1; j<=k; j++)
	{
        alph[j + 1] = sqrt(alph[j] * bet[j]);
        bet[j + 1] = half * (alph[j] + bet[j]);
    }
    s[1][1] = sqrt(alph[k1] * bet[k1]);
    for (j = 1; j<= k; j++)
	{
        ab = alph[k1 - j] * bet[k1 - j];
        for (n = 1; n<=pow( 2 , (j - 1)); n++)
		{
            disc = sqrt(pow(s[n][j] , 2) - ab);
            s[2 * n][j + 1] = s[n][j] + disc;
            s[2 * n - 1][j + 1] = ab / s[2 * n][j + 1];
        }
    }
    for (n = 1; n<=nr; n++)
        r[n] = s[n][k1];
    
    anormg = zero;
    for (j = 2; j<=jmax - 1; j++)
	{
      for (l = 2; l<=jmax - 1; l++)
	  {
        anormg = anormg + fabs(g[j][l]);
        psi[j][l] = -d[j][l] * u[j][l - 1] + (r[1] - e[j][l]) * u[j][l];
        psi[j][l] = psi[j][l] - f[j][l] * u[j][l + 1];
      }
    }
    nits = maxits / nr;
    for (kits = 1; kits<= nits; kits++)
	{
        for (n = 1; n<=nr; n++)
		{
            if (n == nr)
                next1 = 1;
            else
                next1 = n + 1;
            rfact = r[n] + r[next1];
            for (l = 2; l<= jmax - 1; l++)
			{
                for (j = 2; j<= jmax - 1; j++)
				{
                    aa[j - 1] = a[j][l];
                    bb[j - 1] = b[j][l] + r[n];
                    cc[j - 1] = c[j][l];
                    rr[j - 1] = psi[j][l] - g[j][l];
                }
                tridag(aa,bb,cc,rr,uu,jmax - 2);
                for (j = 2; j<=jmax - 1; j++)
                    psi[j][l] = -psi[j][l] + two * r[n] * uu[j - 1];
               
            }
              for (j = 2; j<= jmax - 1; j++)
			  {
					for (l = 2 ;l<= jmax - 1;l++)
					{
                    aa[l - 1] = d[j][l];
                    bb[l - 1] = e[j][l] + r[n];
                    cc[l - 1] = f[j][l];
                    rr[l - 1] = psi[j][l];
					}
                tridag(aa,bb,cc,rr,uu,jmax - 2);
                for (l = 2; l<= jmax - 1; l++)
				{
                    u[j][l] = uu[l - 1];
                    psi[j][l] = -psi[j][l] + rfact * uu[l - 1];
                }
			  }
        }
        anorm = zero;
        for (j = 2; j<= jmax - 1; j++)
		{
          for (l = 2; l<= jmax - 1; l++)
		  {
            resid = a[j][l] * u[j - 1][l] + (b[j][l] + e[j][l]) * u[j][l];
            resid = resid + c[j][l] * u[j + 1][l];
            resid = resid + d[j][l] * u[j][l - 1] * u[j][l - 1];
            resid = resid + f[j][l] * u[j][l + 1] + g[j][l];
            anrom = anorm + fabs(resid);
          }
        }
        if (anorm < eps * anormg) return;
    }
    cout<< " maxits exceeded"<<endl;
}

void main()
{
    //program d16r2
    //driver for routine ADI
	int jmax,i,j,mid1,k;
	double pi,aaa,alpha,beta,alim,eps;
    jmax = 11;
    pi = 3.1415926;
    double a[12][12],b[12][12],c[12][12],d[12][12];
    double e[12][12],f[12][12],u[12][12],g[12][12];
    for (i=1; i<=jmax; i++)
	{
        for (j = 1; j<=jmax; j++)
		{
            a[i][j] = -1.0;
            b[i][j] = 2.0;
            c[i][j] = -1.0;
            d[i][j] = -1.0;
            e[i][j] = 2.0;
            f[i][j] = -1.0;
            g[i][j] = 0.0;
            u[i][j] = 0.0;
        }
    }
    mid1 = jmax / 2 + 1;
    g[mid1][mid1] = 2.0;
    alpha = 2.0 * (1.0 - cos(pi / jmax));
    beta = 2.0 * (1.0 - cos((jmax - 1) * pi / jmax));
    alim = log(4.0 * jmax / pi);
    k = 0;
    do
        k = k + 1;
    while (pow(2 , k) < alim);
    eps = 0.0001;
    adi(a,b,c,d,e,f,g,u,jmax,k,alpha,beta,eps);
    cout<<endl;
    cout<<"ADI 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-1 ; j++)
		{
			aaa =  -4.0 * u[i][j] + u[i + 1][j];
			g[i][j] =aaa + u[i - 1][j] + u[i][j - 1] + u[i][j + 1];
		}
		for (j = 2; j<=jmax-1 ; j++)
			cout<<setw(8)<<g[i][j];
		cout<<endl;
    }
}