GRASS Programmer's Manual  6.4.2(2012)
lu.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines