gs.c
上传用户:bjtelijie
上传日期:2010-01-01
资源大小:87k
文件大小:2k
源码类别:

数学计算

开发平台:

Visual C++

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <malloc.h>
  4. #include <math.h>
  5. int GS(int,double**,double *,double);
  6. double **TwoArrayAlloc(int,int);
  7. void TwoArrayFree(double **);
  8. void main()
  9. {
  10. int i,n;
  11. double ep,**a,*b;
  12. n = 3;
  13. ep = 1e-4;
  14. a = TwoArrayAlloc(n,n);
  15. b = (double *)calloc(n,sizeof(double));
  16. if(b == NULL)
  17. {
  18. printf("内存分配失败n");
  19. exit(1);
  20. }
  21. a[0][0]= 2; a[0][1]= 6; a[0][2]=-1;
  22. a[1][0]= 5; a[1][1]=-1; a[1][2]= 2;
  23. a[2][0]=-3; a[2][1]=-4; a[2][2]= 1;
  24. b[0] = -12; b[1] = 29;  b[2] = 5;
  25. if(!GS(n,a,b,ep))
  26. {
  27. printf("不可以用高斯消去法求解n");
  28. exit(0);
  29. }
  30. printf("该方程组的解为:n");
  31. for(i=0;i<3;i++)
  32. printf("x%d = %.2fn",i,b[i]);
  33. TwoArrayFree(a);
  34. free(b);
  35. }
  36. int GS(n,a,b,ep)
  37. int n;
  38. double **a;
  39. double *b;
  40. double ep;
  41. {
  42. int i,j,k,l;
  43. double t;
  44. for(k=1;k<=n;k++)
  45. {
  46. for(l=k;l<=n;l++)
  47. if(fabs(a[l-1][k-1])>ep)
  48. break;
  49. else if(l==n)
  50. return(0);
  51. if(l!=k)
  52. {
  53. for(j=k;j<=n;j++)
  54. {
  55. t = a[k-1][j-1];
  56. a[k-1][j-1]=a[l-1][j-1];
  57. a[l-1][j-1]=t;
  58. }
  59. t=b[k-1];
  60. b[k-1]=b[l-1];
  61. b[l-1]=t;
  62. }
  63. t=1/a[k-1][k-1];
  64. for(j=k+1;j<=n;j++)
  65. a[k-1][j-1]=t*a[k-1][j-1];
  66. b[k-1]*=t;
  67. for(i=k+1;i<=n;i++)
  68. {
  69. for(j=k+1;j<=n;j++)
  70. a[i-1][j-1]-=a[i-1][k-1]*a[k-1][j-1];
  71. b[i-1]-=a[i-1][k-1]*b[k-1];
  72. }
  73. }
  74. for(i=n-1;i>=1;i--)
  75. for(j=i+1;j<=n;j++)
  76. b[i-1]-=a[i-1][j-1]*b[j-1];
  77. return(1);
  78. }
  79. double **TwoArrayAlloc(int r,int c)
  80. {
  81. double *x,**y;
  82. int n;
  83. x=(double *)calloc(r*c,sizeof(double));
  84. y=(double **)calloc(r,sizeof(double*));
  85. for(n=0;n<=r-1;++n)
  86. y[n]=&x[c*n];
  87. return (y);
  88. }
  89. void TwoArrayFree(double **x)
  90. {
  91. free(x[0]);
  92. free(x);
  93. }