用矩阵直接三角分解法求解方程组

 1 //用矩阵直接三角分解法求解方程组
 2 #include<iostream>
 3 #include<stdio.h>
 4 using namespace std;
 5 void main()
 6 {
 7     int n,i;
 8     float *a,*x;
 9     cout<<"请输入方程组元的个数:";
10     cin>>n;
11 
12     a=new float[n*(n+1)];       //开辟增广矩阵的存储空间
13     x=new float[n];             //开辟方程组解的存储空间
14 
15     cout<<"请输入方程组的系数矩阵："<<endl;
16     for(i=0;i<n*(n+1);i++)
17     {
18          if ((i+1)%(n+1)==0) continue;
19          cin>>*(a+i);
20     }
21 
22     cout<<"请输入方程组的常数项："<<endl;
23     for(i=n;i<n*(n+1);i+=n+1)
24       cin>>*(a+i);
25 
26     void DirectLU(float*,int,float*);
27     DirectLU(a,n,x);
28     cout<<"方程组的解是："<<endl;
29     for(i=0;i<n;i++)
30     cout<<"x["<<i<<"]="<<float(*(x+i))<<endl;
31 }
32 void DirectLU(float*u,int n,float*x)
33 {
34     int i,r,k;
35     for(r=0;r<=n-1;r++)
36     {
37         for(i=r;i<=n;i++)
38         for(k=0;k<=r-1;k++)
39             *(u+r*(n+1)+i)-=*(u+r*(n+1)+k)*(*(u+k*(n+1)+i));
40         for(i=r+1;i<=n-1;i++)
41         {
42             for(k=0;k<=r-1;k++)
43                 *(u+i*(n+1)+r)-=*(u+i*(n+1)+k)*(*(u+k*(n+1)+r));
44             *(u+i*(n+1)+r)/=*(u+r*(n+1)+r);
45         }
46     }
47 
48     for(i=n-1;i>=0;i--)
49     {
50         for(r=n-1;r>=i+1;r--)
51             *(u+i*(n+1)+n)-=*(u+i*(n+1)+r)*x[r];
52         x[i]=*(u+i*(n+1)+n)/(*(u+i*(n+1)+i));
53     }
54 
55 }