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