home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / math / ols / lustuff.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-28  |  8.9 KB  |  434 lines
  1. /*
  2.     The LU approach to numerical linear algebra.
  3.     Mostly off Numerical Recipes.
  4. */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "utils.h"
  9. #include "matrix.h"
  10.  
  11. void lubksb (float **a, int n, int *indx, float b[]);
  12. void d_lubksb (double **a, int n, int *indx, double b[]);
  13. int ludcmp (float **a, int n, int *indx, float *d);
  14. int d_ludcmp (double **a, int n, int *indx, double *d);
  15. int dnear (double x, double y);
  16. int fnear (float x, float y);
  17.  
  18. int d_LinSys (double **A, double *b, double *x, int K);
  19.  
  20. #ifdef TESTING
  21. int TestInversion ();
  22. int TestLinSys ();
  23.  
  24. int 
  25. main ()
  26. {
  27.   int returncode;
  28.  
  29.   returncode = 0;
  30.   if (1 == TestInversion ())
  31.     {
  32.       returncode = 1;
  33.       printf ("INVERSION TEST FAILED.\n");
  34.     }
  35.   if (1 == TestLinSys ())
  36.     {
  37.       returncode = 1;
  38.       printf ("LINEAR SYSTEMS TEST FAILED.\n");
  39.     }
  40.   return returncode;
  41. }
  42.  
  43. int 
  44. TestInversion ()
  45. {
  46.   double biggesterror, error, **A, **Acopy, **B, **product;
  47.   int trial, i, j, K, returncode, seed;
  48.  
  49.   seed = 0;
  50.   K = 20;            /* going to be playing with K x K matrices.  */
  51.   A = dmatrix (1, K, 1, K);
  52.   Acopy = dmatrix (1, K, 1, K);
  53.   B = dmatrix (1, K, 1, K);
  54.   product = dmatrix (1, K, 1, K);
  55.   printf ("TESTING MATRIX INVERSION.\n");
  56.   printf ("Today, we're going to play with %d x %d matrices.\n\n", K, K);
  57.  
  58.   for (trial = 1; trial <= 5; trial++)
  59.     {
  60.       printf ("Trial #%d:\n    generating matrix,\n", trial);
  61.       RandomSym (A, K, &seed);
  62.       CopyMatrix (A, Acopy, K, K);    /* because inversion clobbers A */
  63.       printf ("    inverting it,\n");
  64.       returncode = d_Inverse (A, B, K);
  65.       if (returncode != 0)
  66.     {
  67.       printf ("d_Inverse returned %d.\n", returncode);
  68.       return 1;
  69.     }
  70.       printf ("    post-multiplying,\n");
  71.       d_MMult (Acopy, B, product, K, K, K);
  72.       printf ("    finding biggest error = ");
  73.       biggesterror = -1.0e10;
  74.       for (i = 1; i <= K; i++)
  75.     for (j = 1; j <= K; j++)
  76.       {
  77.         if (i == j)
  78.           error = fabs (product[i][i] - 1.0);
  79.         else
  80.           error = fabs (product[i][j]);
  81.         if (error > biggesterror)
  82.           biggesterror = error;
  83.       }
  84.       printf ("%e.\n\n", biggesterror);
  85.     }
  86.   free_dmatrix (A, 1, K, 1, K);
  87.   free_dmatrix (B, 1, K, 1, K);
  88.   free_dmatrix (Acopy, 1, K, 1, K);
  89.   free_dmatrix (product, 1, K, 1, K);
  90.   return 0;
  91. }
  92.  
  93. int 
  94. TestLinSys ()
  95. {
  96.   double sum, biggesterror, error, **A, **Acopy;
  97.   double *x, *b, *test;
  98.   int row, trial, i, K, returncode, seed;
  99.  
  100.   seed = 0;
  101.   K = 20;            /* going to be playing with K x K matrices.  */
  102.   A = dmatrix (1, K, 1, K);
  103.   Acopy = dmatrix (1, K, 1, K);
  104.   x = dvector (1, K);
  105.   b = dvector (1, K);
  106.   test = dvector (1, K);
  107.   printf ("TESTING SOLUTION OF LINEAR SYSTEMS.\n");
  108.   printf ("Today, we're going to play with %d x %d systems.\n\n", K, K);
  109.  
  110.   for (trial = 1; trial <= 5; trial++)
  111.     {
  112.       printf ("Trial #%d:\n    generating system,\n", trial);
  113.       RandomSym (A, K, &seed);
  114.       CopyMatrix (A, Acopy, K, K);    /* because LinSys clobbers A. */
  115.       for (i = 1; i <= K; i++)
  116.     b[i] = ran1 (&seed);
  117.       printf ("    solving system,\n");
  118.       returncode = d_LinSys (A, b, x, K);
  119.       if (returncode != 0)
  120.     {
  121.       printf ("d_LinSys returned %d.\n", returncode);
  122.       return 1;
  123.     }
  124.       printf ("    post-multiplying,\n");
  125.       for (row = 1; row <= K; row++)
  126.     {
  127.       sum = 0.0;
  128.       for (i = 1; i <= K; i++)
  129.         {
  130.           sum += Acopy[row][i] * x[i];
  131.           test[row] = sum;
  132.         }
  133.     }
  134.       printf ("    finding biggest error = ");
  135.       biggesterror = -1.0e10;
  136.       for (i = 1; i <= K; i++)
  137.     {
  138.       error = fabs (test[i] - b[i]);
  139.       if (error > biggesterror)
  140.         biggesterror = error;
  141.     }
  142.       printf ("%e.\n\n", biggesterror);
  143.     }
  144.   free_dmatrix (Acopy, 1, K, 1, K);
  145.   free_dmatrix (A, 1, K, 1, K);
  146.   free_dvector (x, 1, K);
  147.   free_dvector (b, 1, K);
  148.   free_dvector (test, 1, K);
  149.   return 0;
  150. }
  151.  
  152. #endif
  153.  
  154. /* all routines are in matrix.h except low-level lu routines. */
  155.  
  156. #define TINY 1.0e-10;
  157. #define DTINY 1.0e-20;
  158.  
  159. int 
  160. fnear (float x, float y)
  161. {
  162.   return (fabs (x - y) < 1.0e-10) ? 1 : 0;
  163. }
  164.  
  165. int 
  166. dnear (double x, double y)
  167. {
  168.   return (fabs (x - y) < 1.0e-20) ? 1 : 0;
  169. }
  170.  
  171. void 
  172. lubksb (float **a, int n, int *indx, float b[])
  173. {
  174.   int i, ii = 0, ip, j;
  175.   float sum;
  176.  
  177.   for (i = 1; i <= n; i++)
  178.     {
  179.       ip = indx[i];
  180.       sum = b[ip];
  181.       b[ip] = b[i];
  182.       if (ii)
  183.     for (j = ii; j <= i - 1; j++)
  184.       sum -= a[i][j] * b[j];
  185.       else if (sum)
  186.     ii = i;
  187.       b[i] = sum;
  188.     }
  189.   for (i = n; i >= 1; i--)
  190.     {
  191.       sum = b[i];
  192.       for (j = i + 1; j <= n; j++)
  193.     sum -= a[i][j] * b[j];
  194.       b[i] = sum / a[i][i];
  195.     }
  196. }
  197.  
  198. void 
  199. d_lubksb (double **a, int n, int *indx, double b[])
  200. {
  201.   int i, ii = 0, ip, j;
  202.   double sum;
  203.  
  204.   for (i = 1; i <= n; i++)
  205.     {
  206.       ip = indx[i];
  207.       sum = b[ip];
  208.       b[ip] = b[i];
  209.       if (ii)
  210.     for (j = ii; j <= i - 1; j++)
  211.       sum -= a[i][j] * b[j];
  212.       else if (sum)
  213.     ii = i;
  214.       b[i] = sum;
  215.     }
  216.   for (i = n; i >= 1; i--)
  217.     {
  218.       sum = b[i];
  219.       for (j = i + 1; j <= n; j++)
  220.     sum -= a[i][j] * b[j];
  221.       b[i] = sum / a[i][i];
  222.     }
  223. }
  224.  
  225. int 
  226. ludcmp (float **a, int n, int *indx, float *d)
  227. {
  228.   int i, imax, j, k;
  229.   float big, dum, sum, temp;
  230.   float *vv;
  231.  
  232.   vv = vector (1, n);
  233.   *d = 1.0;
  234.   for (i = 1; i <= n; i++)
  235.     {
  236.       big = 0.0;
  237.       for (j = 1; j <= n; j++)
  238.     if ((temp = fabs (a[i][j])) > big)
  239.       big = temp;
  240.       if (fnear (big, 0.0))
  241.     return 1;
  242.       vv[i] = 1.0 / big;
  243.     }
  244.   for (j = 1; j <= n; j++)
  245.     {
  246.       for (i = 1; i < j; i++)
  247.     {
  248.       sum = a[i][j];
  249.       for (k = 1; k < i; k++)
  250.         sum -= a[i][k] * a[k][j];
  251.       a[i][j] = sum;
  252.     }
  253.       big = 0.0;
  254.       for (i = j; i <= n; i++)
  255.     {
  256.       sum = a[i][j];
  257.       for (k = 1; k < j; k++)
  258.         sum -= a[i][k] * a[k][j];
  259.       a[i][j] = sum;
  260.       if ((dum = vv[i] * fabs (sum)) >= big)
  261.         {
  262.           big = dum;
  263.           imax = i;
  264.         }
  265.     }
  266.       if (j != imax)
  267.     {
  268.       for (k = 1; k <= n; k++)
  269.         {
  270.           dum = a[imax][k];
  271.           a[imax][k] = a[j][k];
  272.           a[j][k] = dum;
  273.         }
  274.       *d = -(*d);
  275.       vv[imax] = vv[j];
  276.     }
  277.       indx[j] = imax;
  278.       if (a[j][j] == 0.0)
  279.     a[j][j] = TINY;
  280.       if (j != n)
  281.     {
  282.       dum = 1.0 / (a[j][j]);
  283.       for (i = j + 1; i <= n; i++)
  284.         a[i][j] *= dum;
  285.     }
  286.     }
  287.   free_vector (vv, 1, n);
  288.   return 0;
  289. }
  290.  
  291. int 
  292. d_ludcmp (double **a, int n, int *indx, double *d)
  293. {
  294.   int i, imax, j, k;
  295.   double big, dum, sum, temp;
  296.   double *vv;
  297.  
  298.   vv = dvector (1, n);
  299.   *d = 1.0;
  300.   for (i = 1; i <= n; i++)
  301.     {
  302.       big = 0.0;
  303.       for (j = 1; j <= n; j++)
  304.     if ((temp = fabs (a[i][j])) > big)
  305.       big = temp;
  306.       if (dnear (big, 0.0))
  307.     return 1;
  308.       vv[i] = 1.0 / big;
  309.     }
  310.   for (j = 1; j <= n; j++)
  311.     {
  312.       for (i = 1; i < j; i++)
  313.     {
  314.       sum = a[i][j];
  315.       for (k = 1; k < i; k++)
  316.         sum -= a[i][k] * a[k][j];
  317.       a[i][j] = sum;
  318.     }
  319.       big = 0.0;
  320.       for (i = j; i <= n; i++)
  321.     {
  322.       sum = a[i][j];
  323.       for (k = 1; k < j; k++)
  324.         sum -= a[i][k] * a[k][j];
  325.       a[i][j] = sum;
  326.       if ((dum = vv[i] * fabs (sum)) >= big)
  327.         {
  328.           big = dum;
  329.           imax = i;
  330.         }
  331.     }
  332.       if (j != imax)
  333.     {
  334.       for (k = 1; k <= n; k++)
  335.         {
  336.           dum = a[imax][k];
  337.           a[imax][k] = a[j][k];
  338.           a[j][k] = dum;
  339.         }
  340.       *d = -(*d);
  341.       vv[imax] = vv[j];
  342.     }
  343.       indx[j] = imax;
  344.       if (a[j][j] == 0.0)
  345.     a[j][j] = DTINY;
  346.       if (j != n)
  347.     {
  348.       dum = 1.0 / (a[j][j]);
  349.       for (i = j + 1; i <= n; i++)
  350.         a[i][j] *= dum;
  351.     }
  352.     }
  353.   free_dvector (vv, 1, n);
  354.   return 0;
  355. }
  356.  
  357. int 
  358. Inverse (float **A, float **Y, int n)
  359. {
  360.   int *indexes;
  361.   int i, j;
  362.   float d;
  363.   float *onecolumn;
  364.  
  365.   indexes = ivector (1, n);
  366.   if (ludcmp (A, n, indexes, &d) == 1)
  367.     return 1;
  368.   onecolumn = vector (1, n);
  369.  
  370.   for (i = 1; i <= n; i++)
  371.     {
  372.       for (j = 1; j <= n; j++)
  373.     onecolumn[j] = 0.0;
  374.       onecolumn[i] = 1.0;
  375.       lubksb (A, n, indexes, onecolumn);
  376.       for (j = 1; j <= n; j++)
  377.     Y[j][i] = onecolumn[j];
  378.     }
  379.   free_ivector (indexes, 1, n);
  380.   free_vector (onecolumn, 1, n);
  381.   return 0;
  382. }
  383.  
  384. int 
  385. d_Inverse (double **A, double **Y, int n)
  386. {
  387.   int *indexes;
  388.   int i, j;
  389.   double d;
  390.   double *onecolumn;
  391.  
  392.   indexes = ivector (1, n);
  393.   if (d_ludcmp (A, n, indexes, &d) == 1)
  394.     {
  395.       free_ivector (indexes, 1, n);
  396.       return 1;
  397.     }
  398.   onecolumn = dvector (1, n);
  399.  
  400.   for (i = 1; i <= n; i++)
  401.     {
  402.       for (j = 1; j <= n; j++)
  403.     onecolumn[j] = 0.0;
  404.       onecolumn[i] = 1.0;
  405.       d_lubksb (A, n, indexes, onecolumn);
  406.       for (j = 1; j <= n; j++)
  407.     Y[j][i] = onecolumn[j];
  408.     }
  409.   free_ivector (indexes, 1, n);
  410.   free_dvector (onecolumn, 1, n);
  411.   return 0;
  412. }
  413.  
  414. int 
  415. d_LinSys (double **A, double *b, double *x, int K)
  416. {
  417.   int *indexes;
  418.   double d;
  419.  
  420.   indexes = ivector (1, K);
  421.   if (1 == d_ludcmp (A, K, indexes, &d))
  422.     {
  423.       free_ivector (indexes, 1, K);
  424.       return 1;
  425.     }
  426.   CopyVector (b, x, K);        /* copy b into x */
  427.   d_lubksb (A, K, indexes, x);
  428.   free_ivector (indexes, 1, K);
  429.   return 0;
  430. }
  431.  
  432. #undef TINY
  433. #undef DTINY
  434.