用LDLT分解求解方程组c++

#include "iostream.h" void main() { int i,j,k,m,n ; cout<<"请输入未知数的个数："<<endl; cin>>n; cout<<"请输入方程的个数："<<endl; cin>>m; double** a=new double*[m+1]; //a存放系数 double *b =new double[n+1];//方程的值 double** l=new double*[m+1];//L double *d =new double[n+1];//D double *x =new double[n+1];//结果 double *y =new double[n+1]; for(i = 1; i<m+1; i++) { l[i]=new double[n+1]; a[i]=new double[n+1]; } //初始化 cout<<"请依次输入"<<endl<<endl; for(i=1;i<m+1;i++) { cout<<"第"<<i<<"个方程的系数"<<endl; for(j=1;j<n+1;j++) { l[i][j]=0; cin>>a[i][j]; } cin>>b[i]; } double temp; for(i=1;i<n+1;i++) { temp=0; for(k=1;k<i;k++) temp=temp+l[i][k]*l[i][k]*d[k]; d[i]=a[i][i]-temp; for(j=i+1;j<n+1;j++) //分解A=LDLT { temp=0; for(k=1;k<i;k++) temp=temp+l[j][k]*d[k]*l[i][k]; l[j][i]=(a[j][i]-temp)/d[i]; } } //前代解方程组 y[1]=b[1]; for(i=2;i<n+1;i++) { temp=0; for(k=1;k<i;k++) temp=temp+l[i][k]*y[k]; y[i]=b[i]-temp; } //回代 x[n]=y[n]/d[n]; for(i=n-1;i>0;i--) { temp=0; for(k=i+1;k<n+1;k++) temp=temp+l[k][i]*x[k]; x[i]=y[i]/d[i]-temp; } cout<<"方程组的解为："; for (i = 1; i < n+1; i++) { cout<<x[i]<<" "; } cout<<endl; }