lu.c 1.6 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394
  1. #include <math.h>
  2. #include <grass/gis.h>
  3. #define TINY 1.0e-20;
  4. int G_ludcmp(double **a, int n, int *indx, double *d)
  5. {
  6. int i, imax = 0, j, k;
  7. double big, dum, sum, temp;
  8. double *vv;
  9. vv = G_alloc_vector(n);
  10. *d = 1.0;
  11. for (i = 0; i < n; i++) {
  12. big = 0.0;
  13. for (j = 0; j < n; j++)
  14. if ((temp = fabs(a[i][j])) > big)
  15. big = temp;
  16. if (big == 0.0) {
  17. *d = 0.0;
  18. return 0; /* Singular matrix */
  19. }
  20. vv[i] = 1.0 / big;
  21. }
  22. for (j = 0; j < n; j++) {
  23. for (i = 0; i < j; i++) {
  24. sum = a[i][j];
  25. for (k = 0; k < i; k++)
  26. sum -= a[i][k] * a[k][j];
  27. a[i][j] = sum;
  28. }
  29. big = 0.0;
  30. for (i = j; i < n; i++) {
  31. sum = a[i][j];
  32. for (k = 0; k < j; k++)
  33. sum -= a[i][k] * a[k][j];
  34. a[i][j] = sum;
  35. if ((dum = vv[i] * fabs(sum)) >= big) {
  36. big = dum;
  37. imax = i;
  38. }
  39. }
  40. if (j != imax) {
  41. for (k = 0; k < n; k++) {
  42. dum = a[imax][k];
  43. a[imax][k] = a[j][k];
  44. a[j][k] = dum;
  45. }
  46. *d = -(*d);
  47. vv[imax] = vv[j];
  48. }
  49. indx[j] = imax;
  50. if (a[j][j] == 0.0)
  51. a[j][j] = TINY;
  52. if (j != n) {
  53. dum = 1.0 / (a[j][j]);
  54. for (i = j + 1; i < n; i++)
  55. a[i][j] *= dum;
  56. }
  57. }
  58. G_free_vector(vv);
  59. return 1;
  60. }
  61. #undef TINY
  62. void G_lubksb(double **a, int n, int *indx, double b[])
  63. {
  64. int i, ii, ip, j;
  65. double sum;
  66. ii = -1;
  67. for (i = 0; i < n; i++) {
  68. ip = indx[i];
  69. sum = b[ip];
  70. b[ip] = b[i];
  71. if (ii >= 0)
  72. for (j = ii; j < i; j++)
  73. sum -= a[i][j] * b[j];
  74. else if (sum)
  75. ii = i;
  76. b[i] = sum;
  77. }
  78. for (i = n - 1; i >= 0; i--) {
  79. sum = b[i];
  80. for (j = i + 1; j < n; j++)
  81. sum -= a[i][j] * b[j];
  82. b[i] = sum / a[i][i];
  83. }
  84. }