求助！关于cholesky矩阵分解的代码，但是无法输出结果，望帮忙解决，十分感谢！！！

我爱GM

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. #include <malloc.h>
   
4. #include <math.h>
   
5. double**MonAlloc(int n,int m)
   
6. {
   
7.     int i,j;
   
8.     double *vect,**mat;
   
9.     vect=(double*)malloc((n*m)*sizeof(double));
   
10.     mat=(double**)malloc(n*sizeof(double));
    
11.     for(i=0;i<n;i++)
    
12.     {
    
13.         mat[i]=vect;
    
14.         vect=vect+m;
    
15.     }
    
16.     return(mat);
    
17. }
    
18. 
    
19. 
    
20. double* sol(int N,double** A,double* b)
    
21. {
    
22.     int i,j;
    
23.     double sum,*x;
    
24.     x=(double*)malloc(N*sizeof(double));
    
25.     for(i=0;i<N;i++)
    
26.     {
    
27.         sum=0.0;
    
28.         for(j=0;j<i;j++)
    
29.         {
    
30.             sum=sum + A[i][j]*x[j];
    
31.         }
    
32.         x[i]=(b[i]-sum)/A[i][i];
    
33.     }
    
34.     return(x);
    
35. }
    
36. 
    
37. 
    
38. double* Cholesky(int N,double** A,double* b)
    
39. {
    
40.     int i,j,k;
    
41.     double **L,**temp,*x,*y,sum;
    
42.     L=MonAlloc(N,N);
    
43.     temp=MonAlloc(N,N);
    
44.     L[0][0]=sqrt(A[0][0]);
    
45.     for(k=0;k<N;k++)
    
46.     {
    
47.         sum=0;
    
48.         for(i=0;i<k;i++)
    
49.             sum+=sum+L[k][j]*L[k][i];
    
50.         sum=A[k][k] - sum;
    
51.         L[k][k] = sqrt(sum>0? sum : 0);
    
52.         for(i=k+1;i<N;i++)
    
53.         {
    
54.             sum=0;
    
55.             for(j=0;j<k;j++)
    
56.                 sum += L[i][j]*L[k][j];
    
57.             L[i][k]=(A[i][k]-sum)/L[k][k];
    
58.         }
    
59.         for(int j=0; j<k; j++)
    
60.          L[j][k]=0;
    
61.     }
    
62.     y=sol(N,L,b);
    
63.     for(i = 0;i < N;i++)
    
64.         for(j = 0;j < i;j++)
    
65.         {
    
66.             temp[i][j] = L[i][j];
    
67.             L[i][j] = L[j][i];
    
68.             L[j][i] = temp[i][j];
    
69.         }
    
70.     x=sol(N,L,y);
    
71.     return(x);
    
72. }
    
73. void main()
    
74. {
    
75.     int i,N;
    
76.     double **A,*b,*x;
    
77.     N=4;
    
78.     A=MonAlloc(N,N);
    
79.     x=(double*)malloc(N*sizeof(double));
    
80.     b=(double*)malloc(N*sizeof(double));
    
81.     A[0][0]=10; A[0][1]=7;   A[0][2]=8;   A[0][3]=7;
    
82.     A[1][0]=7;  A[1][1]=5;   A[1][2]=6;   A[1][3]=5;
    
83.     A[2][0]=8;  A[2][1]=6;   A[2][2]=10;  A[2][3]=9;
    
84.     A[3][0]=7;  A[3][1]=5;   A[3][2]=9;   A[3][3]=10;
    
85.     b[0]=32;   b[0]=23;   b[0]=33;   b[0]=31;
    
86.     Cholesky(N,A,b);
    
87. }

复制代码