#include <stdlib.h>
#include <math.h>
/*樓競網站www.LouJing.com
擁有該程序的版權,轉載請保留該版權.
謝謝合作!*/
double* allocMem(int ); //分配內存空間函數
void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值
void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比叠代公式求解x的值
void main()
{
short matrixNum; //矩陣的行數(列數)
double *matrixA; //矩陣A，初始系數矩陣
double *matrixD; //矩陣D為A中的主對角陣
double *matrixL; //矩陣L為A中的下三角陣
double *matrixU; //矩陣U為A中的上三角陣
double *B; //矩陣B為雅可比方法叠代矩陣
double *f; //矩陣f為中間的過渡的矩陣
double *x; //x為壹維數組，存放結果
double *xk; //xk為壹維數組，用來在叠代中使用
double *b; //b為壹維數組，存放方程組右邊系數
int i,j,k;
cout<<"<<請輸入矩陣的行數（列數與行數壹致）>>：";
cin>>matrixNum;
//分別為A、D、L、U、B、f、x、b分配內存空間
matrixA=allocMem(matrixNum*matrixNum);
matrixD=allocMem(matrixNum*matrixNum);
matrixL=allocMem(matrixNum*matrixNum);
matrixU=allocMem(matrixNum*matrixNum);
B=allocMem(matrixNum*matrixNum);
f=allocMem(matrixNum);
x=allocMem(matrixNum);
xk=allocMem(matrixNum);
b=allocMem(matrixNum);
//輸入系數矩陣各元素值
cout<<endl<<endl<<endl<<"<<請輸入矩陣中各元素值（為 "<<matrixNum<<"*"<<matrixNum<<"，***計 "<<matrixNum*matrixNum<<" 個元素）"<<">>："<<endl<<endl;
for(i=0;i<matrixNum;i++)
{
cout<<"請輸入矩陣中第 "<<i+1<<" 行的 "<<matrixNum<<" 個元素：";
for(j=0;j<matrixNum;j++)
cin>>*(matrixA+i*matrixNum+j);
}
//輸入方程組右邊系數b的各元素值
cout<<endl<<endl<<endl<<"<<請輸入方程組右邊系數各元素值，***計 "<<matrixNum<<" 個"<<">>："<<endl<<endl;
for(i=0;i<matrixNum;i++)
cin>>*(b+i);
/* 下面將A分裂為A=D-L-U */
//首先將D、L、U做初始化工作
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
*(matrixD+i*matrixNum+j)=*(matrixL+i*matrixNum+j)=*(matrixU+i*matrixNum+j)=0;
//D、L、U分別得到A的主對角線、下三角和上三角；其中D取逆矩陣、L和U各元素取相反數
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
if(i==j&&*(matrixA+i*matrixNum+j)) *(matrixD+i*matrixNum+j)=1/(*(matrixA+i*matrixNum+j));
else if(i>j) *(matrixL+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
else *(matrixU+i*matrixNum+j)=-*(matrixA+i*matrixNum+j);
//求B矩陣中的元素
for(i=0;i<matrixNum;i++)
for(j=0;j<matrixNum;j++)
{
double temp=0;
for(k=0;k<matrixNum;k++)
temp+=*(matrixD+i*matrixNum+k)*(*(matrixL+k*matrixNum+j)+*(matrixU+k*matrixNum+j));
*(B+i*matrixNum+j)=temp;
}
//求f中的元素
for(i=0;i<matrixNum;i++)
{
double temp=0;
for(j=0;j<matrixNum;j++)
temp+=*(matrixD+i*matrixNum+j)*(*(b+j));
*(f+i)=temp;
}
/* 計算x的初始向量值 */
GaussLineMain(matrixA,x,b,matrixNum);
/* 利用雅可比叠代公式求解xk的值 */
int JacobiTime;
cout<<endl<<endl<<endl<<"<<雅可比叠代開始，請輸入希望叠代的次數>>：";
cin>>JacobiTime;
while(JacobiTime<=0)
{
cout<<"叠代次數必須大於0，請重新輸入：";
cin>>JacobiTime;
}
Jacobi(x,xk,B,f,matrixNum,JacobiTime);
//輸出線性方程組的解 */
cout<<endl<<endl<<endl<<"<<方程組運算結果如下>>"<<endl;
cout.precision(20); //設置輸出精度，以此比較不同叠代次數的結果
for(i=0;i<matrixNum;i++)
cout<<"x"<<i+1<<" = "<<*(xk+i)<<endl;
cout<<endl<<endl<<"求解過程結束..."<<endl<<endl;
//釋放掉所有動態分配的內存
delete [] matrixA;
delete [] matrixD;
delete [] matrixL;
delete [] matrixU;
delete [] B;
delete [] f;
delete [] x;
delete [] xk;
delete [] b;
}
/*--------------------
分配內存空間函數
www.LouJing.com
--------------------*/
double* allocMem(int num)
{
double *head;
if((head=new double[num])==NULL)
{
cout<<"內存空間分配失敗，程序終止運行！"<<endl;
exit(0);
}
return head;
}
/*---------------------------------------------
計算Ax=b中x的初始向量值，采用高斯列主元素消去法
www.LouJing.com
---------------------------------------------*/
void GaussLineMain(double* A,double* x,double* b,int num)
{
int i,j,k;
//***處理num-1行
for(i=0;i<num-1;i++)
{
//首先每列選主元，即最大的壹個
double lineMax=fabs(*(A+i*num+i));
int lineI=i;
for(j=i;j<num;j++)
if(fabs(*(A+j*num+i))>fabs(lineMax)) lineI=j;
//主元所在行和當前處理行做行交換，右系數b也隨之交換
for(j=i;j<num;j++)
{
//A做交換
lineMax=*(A+i*num+j);
*(A+i*num+j)=*(A+lineI*num+j);
*(A+lineI*num+j)=lineMax;
//b中對應元素做交換
lineMax=*(b+i);
*(b+i)=*(b+lineI);
*(b+lineI)=lineMax;
}
if(*(A+i*num+i)==0) continue; //如果當前主元為0，本次循環結束
//將A變為上三角矩陣，同樣b也隨之變換
for(j=i+1;j<num;j++)
{
double temp=-*(A+j*num+i)/(*(A+i*num+i));
for(k=i;k<num;k++)
{
*(A+j*num+k)+=temp*(*(A+i*num+k));
}
*(b+j)+=temp*(*(b+i));
}
}
/* 驗證Ax=b是否有唯壹解，就是驗證A的行列式是否為0；
如果|A|!=0，說明有唯壹解*/
double determinantA=1;
for(i=0;i<num;i++)
determinantA*=*(A+i*num+i);
if(determinantA==0)
{
cout<<endl<<endl<<"通過計算,矩陣A的行列式為|A|=0，即A沒有唯壹解。\n程序退出..."<<endl<<endl;
exit(0);
}
/* 從最後壹行開始，回代求解x的初始向量值 */
for(i=num-1;i>=0;i--)
{
for(j=num-1;j>i;j--)
*(b+i)-=*(A+i*num+j)*(*(x+j));
*(x+i)=*(b+i)/(*(A+i*num+i));
}
}
/*------------------------------------
利用雅可比叠代公式求解x的精確值
www.LouJing.com
叠代公式為：xk=Bx+f
------------------------------------*/
void Jacobi(double* x,double* xk,double* B,double* f,int num,int time)
{
int t=1,i,j;
while(t<=time)
{
for(i=0;i<num;i++)
{
double temp=0;
for(j=0;j<num;j++)
temp+=*(B+i*num+j)*(*(x+j));
*(xk+i)=temp+*(f+i);
}
//將xk賦值給x，準備下壹次叠代
for(i=0;i<num;i++)
*(x+i)=*(xk+i);
t++;
}
}