home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / rlib / lu.r < prev    next >
Encoding:
Text File  |  1994-09-23  |  2.1 KB  |  90 lines  |  [TEXT/RLAB]
  1. //-------------------------------------------------------------------//
  2. //  Syntax:    lu ( A )
  3.  
  4. //  Description:
  5.  
  6. //  Lu computes the "LU" factors of the input MATRIX. The input MATRIX
  7. //  must be square. lu() returns l, u, and pvt in a LIST. The
  8. //  factorization has the form: 
  9.  
  10. //  A = p*l*u            where A is the input matrix
  11.  
  12. //  The function lu is a user-function that utilizes the built-in
  13. //  function factor. Factor does the factorization using the LAPACK
  14. //  functions DGETRF/ZGETRF, and DGECON/ZGECON.
  15.  
  16. //  Example:
  17.  
  18. //  > a = [1,2,3;4,5,6;7,8,0]
  19. //   a =
  20. //   matrix columns 1 thru 3
  21. //             1           2           3
  22. //             4           5           6
  23. //             7           8           0
  24. //  > f = lu(a)
  25. //    l            pvt            u            
  26. //  > f.p
  27. //   matrix columns 1 thru 3
  28. //             0           0           1
  29. //             1           0           0
  30. //             0           1           0
  31. //  > f.l
  32. //   matrix columns 1 thru 3
  33. //             1           0           0
  34. //         0.143           1           0
  35. //         0.571         0.5           1
  36. //  > f.u
  37. //   matrix columns 1 thru 3
  38. //             7           8           0
  39. //             0       0.857           3
  40. //             0           0         4.5
  41. //  > f.p*f.l*f.u
  42. //   matrix columns 1 thru 3
  43. //             1           2           3
  44. //             4           5           6
  45. //             7           8           0
  46.  
  47. //  See Also: backsub, factor, inv, solve
  48. //-------------------------------------------------------------------//
  49.  
  50. static (swap);
  51.  
  52. lu = function ( A )
  53. {
  54.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  55.  
  56.   x = factor (A, "g");    // Do the factorization
  57.  
  58.   //
  59.   // Now create l, u, and pvt from lu and pvt.
  60.   //
  61.  
  62.   l = tril (x.lu, -1) + eye (size (x.lu));
  63.   u = triu (x.lu);
  64.   pvt = eye (size (x.lu));
  65.  
  66.   //
  67.   // Now re-arange the columns of pvt
  68.   //
  69.  
  70.   for (i in 1:max (size (x.lu)))
  71.   {
  72.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  73.   }
  74.   return << l = l; u = u; pvt = pvt >>;
  75. };
  76.  
  77. //
  78. //  In vector V, swap elements I, J
  79. //
  80.  
  81. swap = function ( V, I, J )
  82. {
  83.   local (v, tmp);
  84.   v = V;
  85.   tmp = v[I];
  86.   v[I] = v[J];
  87.   v[J] = tmp;
  88.   return v;
  89. };
  90.