/*
-----------------------------------------------
假设有如下方程组：
Ax=b
用Jacobi迭代法求解方程组的解
方法：将A分裂为A=D-L-U，等价的迭代方程组为x=Bx+f。
有关算法的详细说明，参看http://www.loujing.com/mywork/c++/project/Jacobi.pdf
-----------------------------------------------
*/
#include <iostream.h>
#include <stdlib.h>
#include <math.h>
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;
}
/*
----------------------
分配内存空间函数
----------------------
*/
double* allocMem(int num)
{
double *head;
if((head=new double[num])==NULL)
{
cout<<"内存空间分配失败，程序终止运行！"<<endl;
exit(0);
}
return head;
}
/*
-----------------------------------------------
计算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++)
{
//首先每列选主元，即最大的一个
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的精确值
迭代公式为：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++;
}
}
//程序编写结束
