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

void main()
{
    //program d1r8
    //driver program for routine svbksb,svdcmp
    int i,j,k,l,n = 5;
    double a[6][6], b[6], u[26], v[26], w[6], x[6];
    //輸入已知的方程組的系數矩陣
    a[1][1]=1.4; a[1][2]=2.1; a[1][3]=2.1; a[1][4]=7.4; a[1][5]=9.6;
    a[2][1]=1.6; a[2][2]=1.5; a[2][3]=1.1; a[2][4]=0.7; a[2][5]=5.0;
    a[3][1]=3.8; a[3][2]=8.0; a[3][3]=9.6; a[3][4]=5.4; a[3][5]=8.8;
    a[4][1]=4.6; a[4][2]=8.2; a[4][3]=8.4; a[4][4]=0.4; a[4][5]=8.0;
    a[5][1]=2.6; a[5][2]=2.9; a[5][3]=0.1; a[5][4]=9.9; a[5][5]=7.7;
    //輸入已知的方程組的右端向量b
    b[1] = 1.1;
    b[2] = 1.6;
    b[3] = 4.7;
    b[4] = 9.1;
    b[5] = 0.1;
	cout.setf(ios::fixed|ios::right);
	cout.precision(5);
    cout<<endl;
    cout<<"已知的方程組的右端向量"<<endl;
	cout.width(12);    cout<<b[1]<<endl;
	cout.width(12);    cout<<b[2]<<endl;
	cout.width(12);    cout<<b[3]<<endl;
	cout.width(12);    cout<<b[4]<<endl;
	cout.width(12);    cout<<b[5]<<endl;
    for (i = 1; i<=n; i++)
	{
        for (j = 1; j<=n; j++)
		{
           u[(i-1)*n+j] = a[i][j];
		   //a1[(i-1)*n+j]=a[i][j];
        }
    }
    svdcmp(u, n, n, w, v);
    //find maximum singular vaulue
	double wmax, wmin;
    wmax = 0.0;
    for (k = 1; k<=n; k++)
	{
        if (w[k] > wmax)
		{
			wmax = w[k];
		}
    }
    //define"small"
    wmin = wmax * 0.0000001;
    //zero the "small" singular values
    for (k = 1; k<=n; k++)
	{
        if (w[k] < wmin)
		{
			w[k] = 0.0;
		}
    }
    svbksb(u, w, v, n, n, b, x);
    cout<<endl;
    cout<<"計算出的方程組的解"<<endl;
	cout.width(12);    cout<<x[1]<<endl;
	cout.width(12);    cout<<x[2]<<endl;
	cout.width(12);    cout<<x[3]<<endl;
	cout.width(12);    cout<<x[4]<<endl;
	cout.width(12);    cout<<x[5]<<endl;
    //將計算出的解乘以系數矩陣，以驗證計算結果正確
    for (l = 1; l<=n;l++)
	{
        b[l] = 0.0;
        for (j = 1; j<=n; j++)
		{
            b[l] = b[l] + a[l][j] * x[j];
        }
    }
    cout<<endl;
    cout<<"計算出的解乘以系數矩陣的結果"<<endl;
	cout.width(12);    cout<<b[1]<<endl;
	cout.width(12);    cout<<b[2]<<endl;
	cout.width(12);    cout<<b[3]<<endl;
	cout.width(12);    cout<<b[4]<<endl;
	cout.width(12);    cout<<b[5]<<endl;
}