~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/numeric/solve.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)
  2 **/
  3 #include <math.h>
  4 #include <tina/sys.h>
  5 #include <tina/sysfuncs.h>
  6 #include <tina/math.h>
  7 #include <tina/mathfuncs.h>
  8 
  9 #define TINY 1.0e-20;
 10 
 11 void    ludcmp(double **a, int n, int *indx, double *d)
 12 {
 13     int     i, imax, j, k;      /* BUG `imax' may be used uninitialized
 14                                  * in this function */
 15     double  big, dum, sum, temp;
 16     double *vv;
 17 
 18     vv = dvector_alloc(0, n);
 19     *d = 1.0;
 20     for (i = 0; i < n; i++)
 21     {
 22         big = 0.0;
 23         for (j = 0; j < n; j++)
 24             if ((temp = fabs(a[i][j])) > big)
 25                 big = temp;
 26         if (big == 0.0)
 27             big = TINY;
 28         vv[i] = 1.0 / big;
 29     }
 30     for (j = 0; j < n; j++)
 31     {
 32         for (i = 0; i < j; i++)
 33         {
 34             sum = a[i][j];
 35             for (k = 0; k < i; k++)
 36                 sum -= a[i][k] * a[k][j];
 37             a[i][j] = sum;
 38         }
 39         big = 0.0;
 40         for (i = j; i < n; i++)
 41         {
 42             sum = a[i][j];
 43             for (k = 0; k < j; k++)
 44                 sum -= a[i][k] * a[k][j];
 45             a[i][j] = sum;
 46             if ((dum = vv[i] * fabs(sum)) >= big)
 47             {
 48                 big = dum;
 49                 imax = i;
 50             }
 51         }
 52         if (j != imax)
 53         {
 54             for (k = 0; k < n; k++)
 55             {
 56                 dum = a[imax][k];
 57                 a[imax][k] = a[j][k];
 58                 a[j][k] = dum;
 59             }
 60             *d = -(*d);
 61             vv[imax] = vv[j];
 62         }
 63         indx[j] = imax;
 64         if (a[j][j] == 0.0)
 65             a[j][j] = TINY;
 66         if (j != n - 1)
 67         {
 68             dum = 1.0 / (a[j][j]);
 69             for (i = j + 1; i < n; i++)
 70                 a[i][j] *= dum;
 71         }
 72     }
 73     dvector_free((void *) vv, 0);
 74 }
 75 
 76 #undef TINY
 77 
 78 
 79 void    lubksb(double **a, int n, int *indx, double *b)
 80 {
 81     int     i, ii = -1, ip, j;
 82     double  sum;
 83 
 84     for (i = 0; i < n; i++)
 85     {
 86         ip = indx[i];
 87         sum = b[ip];
 88         b[ip] = b[i];
 89         if (ii >= 0)
 90             for (j = ii; j <= i - 1; j++)
 91                 sum -= a[i][j] * b[j];
 92         else if (sum)
 93             ii = i;
 94         b[i] = sum;
 95     }
 96     for (i = n - 1; i >= 0; i--)
 97     {
 98         sum = b[i];
 99         for (j = i + 1; j < n; j++)
100             sum -= a[i][j] * b[j];
101         b[i] = sum / a[i][i];
102     }
103 }
104 
105 Vector *matrix_solve(Matrix * mat, Vector * y)
106 {
107     Vector *x;
108     double **a, *b, d;
109     int    *index, n;
110 
111     if (mat == NULL || y == NULL)
112         return (NULL);
113     if (mat->m != y->n || mat->n != y->n)
114         return (NULL);
115 
116     n = mat->m;
117     a = darray_dcopy(mat->el.double_v, 0, 0, n, n);
118     x = vector_copy(y);
119     b = x->data;
120     index = ivector_alloc(0, n);
121 
122     ludcmp(a, n, index, &d);
123     lubksb(a, n, index, b);
124 
125     darray_free((char **) a, 0, 0, n, n);
126     ivector_free((void *) index, 0);
127 
128     return (x);
129 }
130 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.