/*
　　-----------------------------------------------
　　假設(shè)有如下方程組：
　　Ax=b
　　用Jacobi迭代法求解方程組的解
　　方法：將A分裂為A=D-L-U，等價(jià)的迭代方程組為x=Bx+f。
　　有關(guān)算法的詳細(xì)說(shuō)明，參看http://www.loujing.com/mywork/c++/project/Jacobi.pdf
　　-----------------------------------------------
　　*/
　　#include <iostream.h>
　　#include <stdlib.h>
　　#include <math.h>
　　double* allocMem(int ); //分配內(nèi)存空間函數(shù)
　　void GaussLineMain(double*,double*,double*,int );//采用高斯列主元素消去法求解x的初始向量值
　　void Jacobi(double*,double*,double*,double*,int,int);//利用雅可比迭代公式求解x的值
　　void main()
　　{
　　short matrixNum; //矩陣的行數(shù)(列數(shù))
　　double *matrixA; //矩陣A，初始系數(shù)矩陣
　　double *matrixD; //矩陣D為A中的主對(duì)角陣
　　double *matrixL; //矩陣L為A中的下三角陣
　　double *matrixU; //矩陣U為A中的上三角陣
　　double *B; //矩陣B為雅可比方法迭代矩陣
　　double *f; //矩陣f為中間的過(guò)渡的矩陣
　　double *x; //x為一維數(shù)組，存放結(jié)果
　　double *xk; //xk為一維數(shù)組，用來(lái)在迭代中使用
　　double *b; //b為一維數(shù)組，存放方程組右邊系數(shù)
　　int i,j,k;
　　cout<<"<<請(qǐng)輸入矩陣的行數(shù)（列數(shù)與行數(shù)一致）>>：";
　　cin>>matrixNum;
　　//分別為A、D、L、U、B、f、x、b分配內(nèi)存空間
　　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);
　　//輸入系數(shù)矩陣各元素值
　　cout<<endl<<endl<<endl<<"<<請(qǐng)輸入矩陣中各元素值（為 "<<matrixNum<<"*"
　　 <<matrixNum<<"，共計(jì) "<<matrixNum*matrixNum<<" 個(gè)元素）"<<">>："<<endl<<endl;
　　for(i=0;i<matrixNum;i++)
　　{
　　 cout<<"請(qǐng)輸入矩陣中第 "<<i+1<<" 行的 "<<matrixNum<<" 個(gè)元素：";
　　 for(j=0;j<matrixNum;j++)
　　 cin>>*(matrixA+i*matrixNum+j);
　　}
　　//輸入方程組右邊系數(shù)b的各元素值
　　cout<<endl<<endl<<endl<<"<<請(qǐng)輸入方程組右邊系數(shù)各元素值，共計(jì) "<<matrixNum<<
　　 " 個(gè)"<<">>："<<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的主對(duì)角線、下三角和上三角；其中D取逆矩陣、L和U各元素取相反數(shù)
　　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;
　　}
　　/* 計(jì)算x的初始向量值 */
　　GaussLineMain(matrixA,x,b,matrixNum);
　　/* 利用雅可比迭代公式求解xk的值 */
　　int JacobiTime;
　　cout<<endl<<endl<<endl<<"<<雅可比迭代開(kāi)始，請(qǐng)輸入希望迭代的次數(shù)>>：";
　　cin>>JacobiTime;
　　while(JacobiTime<=0)
　　{
　　 cout<<"迭代次數(shù)必須大于0，請(qǐng)重新輸入：";
　　 cin>>JacobiTime;
　　}
　　Jacobi(x,xk,B,f,matrixNum,JacobiTime);
　　//輸出線性方程組的解 */
　　cout<<endl<<endl<<endl<<"<<方程組運(yùn)算結(jié)果如下>>"<<endl;
　　cout.precision(20); //設(shè)置輸出精度，以此比較不同迭代次數(shù)的結(jié)果
　　for(i=0;i<matrixNum;i++)
　　 cout<<"x"<<i+1<<" = "<<*(xk+i)<<endl;
　　cout<<endl<<endl<<"求解過(guò)程結(jié)束..."<<endl<<endl;
　　//釋放掉所有動(dòng)態(tài)分配的內(nèi)存
　　delete [] matrixA;
　　delete [] matrixD;
　　delete [] matrixL;
　　delete [] matrixU;
　　delete [] B;
　　delete [] f;
　　delete [] x;
　　delete [] xk;
　　delete [] b;
　　}
　　/*
　　----------------------
　　 分配內(nèi)存空間函數(shù)
　　----------------------
　　*/
　　double* allocMem(int num)
　　{
　　double *head;
　　if((head=new double[num])==NULL)
　　{
　　 cout<<"內(nèi)存空間分配失敗，程序終止運(yùn)行！"<<endl;
　　 exit(0);
　　}
　　return head;
　　}
　　/*
　　-----------------------------------------------
　　 計(jì)算Ax=b中x的初始向量值，采用高斯列主元素消去法
　　-----------------------------------------------
　　*/
　　void GaussLineMain(double* A,double* x,double* b,int num)
　　{
　　int i,j,k;
　　//共處理num-1行
　　for(i=0;i<num-1;i++)
　　{
　　 //首先每列選主元，即最大的一個(gè)
　　 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;
　　 //主元所在行和當(dāng)前處理行做行交換，右系數(shù)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中對(duì)應(yīng)元素做交換
　　 lineMax=*(b+i);
　　 *(b+i)=*(b+lineI);
　　 *(b+lineI)=lineMax;
　　 }
　　 if(*(A+i*num+i)==0) continue; //如果當(dāng)前主元為0，本次循環(huán)結(jié)束
　　 //將A變?yōu)樯先蔷仃嚕瑯觔也隨之變換
　　 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));
　　 }
　　}
　　/* 驗(yàn)證Ax=b是否有唯一解，就是驗(yàn)證A的行列式是否為0；
　　 如果|A|!=0，說(shuō)明有唯一解
　　*/
　　double determinantA=1;
　　for(i=0;i<num;i++)
　　 determinantA*=*(A+i*num+i);
　　if(determinantA==0)
　　{
　　 cout<<endl<<endl<<"通過(guò)計(jì)算,矩陣A的行列式為|A|=0，即A沒(méi)有唯一解。\n程序退出..."<<endl<<endl;
　　 exit(0);
　　}
　　/* 從最后一行開(kāi)始，回代求解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的精確值
　　迭代公式為：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，準(zhǔn)備下一次迭代
　　 for(i=0;i<num;i++)
　　 *(x+i)=*(xk+i);
　　 t++;
　　}
　　}
　　//程序編寫(xiě)結(jié)束