//類的實現
#include"linequ.h"

//the function member of Class Matrix                   //基類成員函數的實現
void Matrix::setMatrix(double *rmatr)                   //設置矩陣
{
	for(int i=0;i<index*index;i++){
		*(MatrixA+i)=rmatr[i];                          //矩陣成員賦初值
	}
}

Matrix::Matrix(int dims)                               //矩陣Matrix類的構造函數
{
	index=dims;                                         //保護數據賦值
	MatrixA=new double[index*index];                    //動態分配內存
}

Matrix::~Matrix()                                       //矩陣Matrix類的析構函數
{
	delete[] MatrixA;                                   //內存釋放
}

void Matrix::printM()                                        //顯示矩陣元素
{
	cout<<"The Matrix is:"<<endl;
	for(int i=0;i<index;i++){
		for(int j=0;j<index;j++)
			cout<<*(MatrixA+i*index+j)<<"  ";
		cout<<endl;
	}
}

//the function member of Class Linequ
Linequ::Linequ(int dims):Matrix(dims)                   //派生類Linequ的構造函數
{                                                       //使用參數調用基類構造函數
	sums=new double[dims];                              //動態內存分配
	solu=new double[dims];
}

Linequ::~Linequ()                                       //派生類Linequ的析構函數
{                                                       //調用基類析構函數
	delete[] sums;                                      //釋放內存
	delete[] solu;
}

void Linequ::setLinequ(double *a,double *b)             //設置線性方程組
{
	setMatrix(a);                                       //調用基類函數
	for(int i=0;i<index;i++)
		sums[i]=b[i];
}

void Linequ::printL()                                   //顯示線性方程組
{
	cout<<"The Line eqution is:"<<endl;
	for(int i=0;i<index;i++){
		for(int j=0;j<index;j++)
			cout<<*(MatrixA+i*index+j)<<"  ";
		cout<<"  "<<sums[i]<<endl;
	}
}

void Linequ::showX()                                     //輸出方程的解
{
	cout<<"The Result is:"<<endl;
	for(int i=0;i<index;i++){
		cout<<"X["<<i+1<<"]="<<solu[i]<<endl;
	}
}

int Linequ::Solve()                                     //全選主元高斯消去法求解方程
{
	int *js,l,k,i,j,is,p,q;
	double d,t;
	js=new int[index];
	l=1;
	for(k=0;k<=index-2;k++){                             //消去過程
		d=0.0;
		for(i=k;i<=index-1;i++)
			for(j=k;j<=index-1;j++){
				t=fabs(MatrixA[i*index+j]);
				if(t>d)
				{d=t;js[k]=j;is=i;}
			}
		if(d+1.0==1.0) l=0;
		else
		{   
			if(js[k]!=k)
				for(i=0;i<=index-1;i++){
					p=i*index+k;q=i*index+js[k];
					t=MatrixA[p];MatrixA[p]=MatrixA[q];MatrixA[q]=t;
				}
			if(is!=k)
			{
				for(j=k;j<=index-1;j++){
					p=k*index+j;q=is*index+j;
					t=MatrixA[p];MatrixA[p]=MatrixA[q];MatrixA[q]=t;
				}
			t=sums[k];sums[k]=sums[is];sums[is]=t;
			}
		}
		if(l==0)
		{
			delete[] js;cout<<"fail"<<endl;
			return(0);
		}
		d=MatrixA[k*index+k];
		for(j=k+1;j<=index-1;j++){
			p=k*index+j;MatrixA[p]=MatrixA[p]/d;}
		sums[k]=sums[k]/d;
		for(i=k+1;i<=index-1;i++){
			for(j=k+1;j<=index-1;j++){
				p=i*index+j;
				MatrixA[p]=MatrixA[p]-MatrixA[i*index+k]*MatrixA[k*index+j];
			}
			sums[i]=sums[i]-MatrixA[i*index+k]*sums[k];
		}
	}
	d=MatrixA[(index-1)*index+index-1];
	if(fabs(d)+1.0==1.0)
	{
		delete[] js;cout<<"fail"<<endl;
		return(0);
	}
	solu[index-1]=sums[index-1]/d;                               //回代過程
	for(i=index-2;i>=0;i--){
		t=0.0;
		for(j=i+1;j<=index-1;j++)
			t=t+MatrixA[i*index+j]*solu[j];
		solu[i]=sums[i]-t;
	}
	js[index-1]=index-1;
	for(k=index-1;k>=0;k--)
		if(js[k]!=k)
		{t=solu[k];solu[k]=solu[js[k]];solu[js[k]]=t;}
	delete[] js;
	return(1);
}