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

void gaussj(double a[], int n, double b[])
{
	int i,j,k,l,ll,irow,icol;
	double big,pivinv,dum;
    int ipiv[50], indxr[50], indxc[50];
    for (j=0;j<=n-1;j++)
	{
        ipiv[j]=0;
    }
    for (i=0;i<=n-1;i++)
	{
        big=0.0;
        for (j=0;j<=n-1;j++)
		{
            if(ipiv[j]!=1)
			{
                for(k=0;k<=n-1;k++)
				{
					if(ipiv[k]==0)
					{
						if(fabs(a[j*n+k])>=big)
						{
							big=fabs(a[j*n+k]);
							irow=j;
							icol=k;
						}
						else if(ipiv[k]>1)
						{
							cout<<"singular matrix";
						}
					}
                }
            }
        }
        ipiv[icol]=ipiv[icol]+1;
        if(irow!=icol)
		{
            for(l=0;l<=n-1;l++)
			{
                dum=(a[irow*n+l]);
                a[irow*n+l]=a[icol*n+l];
                a[icol*n+l]=dum;
            }
            dum=b[irow];
            b[irow]=b[icol];
            b[icol]=dum;
		}
        indxr[i]=irow;
        indxc[i]=icol;
        if(a[icol*n+icol]==0.0)
		{
			cout<< "singular matrix.";
		}
        pivinv=1.0/(a[icol*n+icol]);
        a[icol*n+icol]=1.0;
        for(l=0;l<=n-1;l++)
		{
            a[icol*n+l]=a[icol*n+l]*pivinv;
        }
        b[icol]=b[icol]*pivinv;
        for(ll=0;ll<=n-1;ll++)
        {
            if(ll!=icol)
			{
                dum=a[ll*n+icol];
                a[ll*n+icol]=0.0;
                for(l=0;l<=n-1;l++)
				{
                    a[ll*n+l]=a[ll*n+l]-a[icol*n+l]*dum;
                }
                b[ll]=b[ll]-b[icol]*dum;
            }
        }
	}	
    for(l=n-1;l<=0;l--)
	{
        if(indxr[l]!=indxc[l])
		{
            for(k=0;k<=n-1;k++)
			{
                dum=a[k*n+indxr[l]];
                a[k*n+indxr[l]]=a[k*n+indxc[l]];
                a[k*n+indxr[l]]=dum;
            }
        }
    }
}



void main()
{
    //program d1r1
    //driver program for routine gaussj
	int i,j,l,n;
    n=3;
    double a[3][3], b[3], a1[9], b1[3];
    //輸入已知的方程組的系數矩陣
    a[0][0]=2; a[0][1]=1;  a[0][2]=2;
    a[1][0]=5; a[1][1]=-1; a[1][2]=1;
    a[2][0]=1; a[2][1]=-3; a[2][2]=-4;
    //輸入已知的方程組的右端向量b
    b[0]=5;
    b[1]=8;
    b[2]=-4;
    cout<<endl;
    cout<<"已知的方程組的右端向量"<<endl;
	cout.width(10);
    cout<<b[0]<<endl;
	cout.width(10);
	cout<<b[1]<<endl;
	cout.width(10);
	cout<<b[2]<<endl;
    for(i=0; i<=n-1; i++)
	{
        for(j=0; j<=2; j++)
		{
			a1[i*n+j]=a[i][j];
        }
    }
    gaussj(a1, n, b);
    cout<<"計算出的方程組的解"<<endl;
	cout.width(10);
    cout<<b[0]<<endl;
	cout.width(10);
	cout<<b[1]<<endl;
	cout.width(10);
	cout<<b[2]<<endl;
    //將計算出的解b乘以系數矩陣，以驗證計算結果正確
    for(l=0; l<=n-1; l++)
	{
        b1[l]=0.0;
        for(j=0; j<=n-1; j++)
		{
            b1[l]=b1[l]+a[l][j]*b[j];
        }
    }
    cout<<endl;
    cout<<"計算出的解乘以系數矩陣的結果"<<endl;
	cout.width(10);
    cout<<b1[0]<<endl;
	cout.width(10);
	cout<<b1[1]<<endl;
	cout.width(10);
	cout<<b1[2]<<endl;
}