#include "iostream.h"
#include "math.h"



void asub(double xin[], double xout[])
{
    int i,n = 20;
    xout[1] = xin[1] + 2.0 * xin[2];
    xout[n] = -2.0 * xin[n - 1] + xin[n];
    for (i = 2; i<=n-1; i++)
	{
        xout[i] = -2.0 * xin[i - 1] + xin[i] + 2.0 * xin[i + 1];
    }
}

void atsub(double xin[], double xout[])
{
    int i,n = 20;
    xout[1] = xin[1] - 2.0 * xin[2];
    xout[n] = 2.0 * xin[n - 1] + xin[n];
    for (i = 2; i<=n-1; i++)
	{
        xout[i] = 2.0 * xin[i - 1] + xin[i] - 2.0 * xin[i + 1];
	}
}
void sparse(double b[], int n, double x[], double rsq)
{
    const int nmax = 500;
    double eps2,eps = 0.000001;
    double g[500], h[500], xi[500], xj[500];
	int j;
    eps2 = n * pow(eps, 2);
    int irst = 0;
	r1:   irst = irst + 1;
    asub(x, xi);
    double rp = 0.0;
    double bsq = 0.0;
    for (j = 1; j<=n; j++)
	{
        bsq = bsq + pow(b[j],2);
        xi[j] = xi[j] - b[j];
        rp = rp + pow(xi[j],2);
    }
    atsub(xi, g);
    for (j = 1; j<=n; j++)
	{
        g[j] = -g[j];
        h[j] = g[j];
    }
    for (int iter = 1; iter<=10 * n; iter++)
	{
        asub(h, xi);
        double anum = 0.0;
        double aden = 0.0;
        for (j = 1; j<=n; j++)
		{
            anum = anum + g[j] * h[j];
            aden = aden + pow(xi[j] ,2);
        }
        if (aden == 0.0)
		{
			cout<<"very singular matrix"<<endl;
		}
        anum = anum / aden;
        for (j = 1; j<=n; j++)
		{
            xi[j] = x[j];
            x[j] = x[j] + anum * h[j];
        }
        asub(x, xj);
        double rsq = 0.0;
        for (j = 1; j<=n; j++)
		{
            xj[j] = xj[j] - b[j];
            rsq = rsq + pow(xj[j], 2);
        }
        if (rsq == rp || rsq <= bsq * eps2)
		{
			return;
		}
        if (rsq > rp)
		{
            for (j = 1; j<=n; j++)
			{
                x[j] = xi[j];
            }
            if (irst >= 3)
			{
                return;
			}
            goto r1;
        }
        rp = rsq;
        atsub(xj, xi);
        double gg = 0.0;
        double dgg = 0.0;
        for (j = 1; j<=n; j++)
		{
            gg = gg + pow(g[j],2);
            dgg = dgg + (xi[j] + g[j]) * xi[j];
        }
        if (gg == 0.0)
		{
			return;
		}
        double gam = dgg / gg;
        for (j = 1; j<=n; j++)
		{
            g[j] = -xi[j];
            h[j] = g[j] + gam * h[j];
        }
    }
    cout<<"too many iterations"<<endl;
}

void main()
{
    //program d1r9
	//driver program for routine sparse
    int i,j,n = 20;
    double b[21], x[21], bcmp[21];
    //輸入已知的方程組的右端向量b
    for (i = 1; i<=n; i++)
	{
        x[i] = 0.0;
        b[i] = 1.0;
    }
    b[1] = 3.0;
    b[n] = -1.0;
    //利用下列子程序asub輸入已知的方程組的系數矩陣
	double rsq=0;
    sparse(b, n, x, rsq);
    cout<<endl;
    cout<<"sum - squared residual:"<<rsq<<endl;
    cout<<endl;
	cout.setf(ios::fixed|ios::right);
	cout.precision(5);
    cout<<"計算出的方程組的解"<<endl;
    for (i = 1; i<=4; i++)
	{
        for (j = i * 5 - 4; j<=i * 5; j++)
		{
			cout.width(12);
            cout<<x[j];
        }
        cout<<endl;
    }
    //將計算出的方程組的解乘以系數矩陣，以驗證計算結果正確
    asub(x, bcmp);
    cout<<endl;
    cout<<"計算出的方程組的解乘以系數矩陣的結果"<<endl;
    cout<<endl;
    cout<<"解乘以系數矩陣"<<"        "<<"方程組的右端向量"<<endl;
    for (i = 1; i<=n; i++)
	{
		cout.width(10);
		cout<<bcmp[i]<<"           ";
		cout.width(12);
		cout<<b[i]<<endl;
    }
}