GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 #include <grass/gis.h> 00003 #include <grass/gmath.h> 00004 00005 00006 #define TINY 1.0e-20; 00007 00008 00009 int G_ludcmp(double **a, int n, int *indx, double *d) 00010 { 00011 int i, imax = 0, j, k; 00012 double big, dum, sum, temp; 00013 double *vv; 00014 00015 vv = G_alloc_vector(n); 00016 *d = 1.0; 00017 for (i = 0; i < n; i++) { 00018 big = 0.0; 00019 for (j = 0; j < n; j++) 00020 if ((temp = fabs(a[i][j])) > big) 00021 big = temp; 00022 if (big == 0.0) { 00023 *d = 0.0; 00024 return 0; /* Singular matrix */ 00025 } 00026 vv[i] = 1.0 / big; 00027 } 00028 for (j = 0; j < n; j++) { 00029 for (i = 0; i < j; i++) { 00030 sum = a[i][j]; 00031 for (k = 0; k < i; k++) 00032 sum -= a[i][k] * a[k][j]; 00033 a[i][j] = sum; 00034 } 00035 big = 0.0; 00036 for (i = j; i < n; i++) { 00037 sum = a[i][j]; 00038 for (k = 0; k < j; k++) 00039 sum -= a[i][k] * a[k][j]; 00040 a[i][j] = sum; 00041 if ((dum = vv[i] * fabs(sum)) >= big) { 00042 big = dum; 00043 imax = i; 00044 } 00045 } 00046 if (j != imax) { 00047 for (k = 0; k < n; k++) { 00048 dum = a[imax][k]; 00049 a[imax][k] = a[j][k]; 00050 a[j][k] = dum; 00051 } 00052 *d = -(*d); 00053 vv[imax] = vv[j]; 00054 } 00055 indx[j] = imax; 00056 if (a[j][j] == 0.0) 00057 a[j][j] = TINY; 00058 if (j != n) { 00059 dum = 1.0 / (a[j][j]); 00060 for (i = j + 1; i < n; i++) 00061 a[i][j] *= dum; 00062 } 00063 } 00064 G_free_vector(vv); 00065 00066 return 1; 00067 } 00068 00069 #undef TINY 00070 00071 void G_lubksb(double **a, int n, int *indx, double b[]) 00072 { 00073 int i, ii, ip, j; 00074 double sum; 00075 00076 ii = -1; 00077 for (i = 0; i < n; i++) { 00078 ip = indx[i]; 00079 sum = b[ip]; 00080 b[ip] = b[i]; 00081 if (ii >= 0) 00082 for (j = ii; j < i; j++) 00083 sum -= a[i][j] * b[j]; 00084 else if (sum) 00085 ii = i; 00086 b[i] = sum; 00087 } 00088 for (i = n - 1; i >= 0; i--) { 00089 sum = b[i]; 00090 for (j = i + 1; j < n; j++) 00091 sum -= a[i][j] * b[j]; 00092 b[i] = sum / a[i][i]; 00093 } 00094 }