/*
-----------------------------------------------
假设有如下方程组：
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;
