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;
	double 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;
}