#include <iostream.h>
#include <stdlib.h>
#include <math.h>

/*樓競(jìng)網(wǎng)站www.LouJing.com
擁有該程序的版權(quán),轉(zhuǎn)載請(qǐng)保留該版權(quán).
謝謝合作!*/
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<<"<<雅可比迭代開始，請(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ù)
  www.LouJing.com
--------------------*/
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的初始向量值，采用高斯列主元素消去法
  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++)
    {
        //首先每列選主元，即最大的一個(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);
    }

    /* 從最后一行開始，回代求解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，準(zhǔn)備下一次迭代
        for(i=0;i<num;i++)
            *(x+i)=*(xk+i);
        t++;
    }
}