home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / LINEARAL.PAK / GAUSSIAN.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  5.9 KB  |  192 lines
  1.  
  2. (* :Name: LinearAlgebra`GaussianElimination` *)
  3.  
  4. (* :Title: Gaussian Elimination and Back Substitution *)
  5.  
  6. (* :Author: Jerry B. Keiper *)
  7.  
  8. (* :Summary:
  9.     Implements LU decomposition with partial pivoting.  This package
  10.     is not intended for computational purposes, but rather for
  11.     educational purposes.  This package can be used in conjunction
  12.     with ComputerArithmetic.m or IntervalArithmetic.m.
  13. *)
  14.  
  15. (* :Context: LinearAlgebra`GaussianElimination` *) 
  16.  
  17. (* :Package Version: Mathematica 2.0 *)
  18.  
  19. (* :Copyright: Copyright 1990,  Wolfram Research, Inc.
  20.     Permission is hereby granted to modify and/or make copies of
  21.     this file for any purpose other than direct profit, or as part
  22.     of a commercial product, provided this copyright notice is left
  23.     intact.  Sale, other than for the cost of media, is prohibited.
  24.  
  25.     Permission is hereby granted to reproduce part or all of
  26.     this file, provided that the source is acknowledged.
  27. *)
  28.  
  29. (* :History:
  30.     Revised by Jerry B. Keiper, December, 1990.
  31.     Original version by Jerry B. Keiper, May, 1990.
  32. *)
  33.  
  34. (* :Keywords:
  35.     Gaussian elimination, back substitution, partial pivoting,
  36.     LU decomposition
  37. *)
  38.  
  39. (* :Source: Any elementary numerical analysis textbook. *)
  40.  
  41. (* :Mathematica Version: 2.0 *)
  42.  
  43. (* :Limitation: *)
  44.  
  45. BeginPackage["LinearAlgebra`GaussianElimination`"]
  46.  
  47. LUFactor::usage = 
  48. "LUFactor[mat] gives the LU decomposition along with a pivot list of the
  49. matrix mat.  The calculation is done using Gaussian elimination with
  50. partial pivoting.  The result is returned as a data object with a head of LU.
  51. LUFactor[ ] may be used in the context of ComputerArithmetic.m or
  52. IntervalArithmetic.m."
  53.  
  54. LUSolve::usage =
  55. "LUSolve[lu, b] solves the linear system represented by lu (the matrix) and
  56. right-hand side b.  The \"matrix\" lu must be the data object with a head of
  57. LU resulting from LUFactor[mat] and b must be an ordinary list of the
  58. appropriate size.  LUSolve[ ] may be used in the context of
  59. ComputerArithmetic.m or IntervalArithmetic.m."
  60.  
  61. LU::usage =
  62. "LU[a, pivots] is the data object returned by LUFactor[ ] and is to be given
  63. as the first argument to LUSolve[ ].  The matrix a represents the lower and
  64. upper triangular matrices in the LU decomposition, and pivots is a record
  65. of the pivots used in the decomposition."
  66.  
  67. Begin["LinearAlgebra`GaussianElimination`Private`"]
  68.  
  69. LUFactor[aa_?MatrixQ] :=
  70.     Module[{ans}, ans /; Head[ans = lufactor[aa]] == LU] /;
  71.     (Length[aa] == Length[aa[[1]]])
  72.  
  73. LUSolve[ap_LU, bb_List] :=
  74.     Module[{ans}, ans /; VectorQ[ans = lusolve[ap, bb]]] /;
  75.     (Length[ap[[1]]] == Length[bb])
  76.  
  77. divide[a_, b_] :=
  78.     If[FreeQ[{a, b}, NumericalMath`ComputerArithmetic`ComputerNumber],
  79.     a/b, NumericalMath`ComputerArithmetic`IdealDivide[a, b]];
  80.  
  81. lufactor[aa_] :=
  82.     Module[{a = aa, pivot, ii, iip, i, ip, j, k, mpiv, m, n=Length[aa], tmp},
  83.     pivot = Table[i, {i, n}];
  84.     For[ii=1, ii<=n, ii++,    (* for each row do ... *)
  85.         (* find a pivot *)
  86.         mpiv = Abs[a[[pivot[[ii]], ii]]];
  87.         k = ii;
  88.         For[i=ii+1, i<=n, i++,
  89.         tmp = Abs[a[[pivot[[i]], ii]]];
  90.         If[tmp > mpiv, mpiv = tmp; k = i]
  91.         ];
  92.         tmp = pivot[[ii]];
  93.         pivot[[ii]] = iip = pivot[[k]];
  94.         pivot[[k]] = tmp;
  95.         mpiv = a[[iip, ii]];
  96.         (* calculate and store the multipliers and reduce the other elements. *)
  97.         For[i=ii+1, i<=n, i++,
  98.         ip = pivot[[i]];
  99.         m = a[[ip, ii]] = divide[a[[ip, ii]], mpiv];
  100.         For[j=ii+1, j<=n, j++, a[[ip, j]] -= m*a[[iip, j]]]
  101.         ];
  102.         ];
  103.     LU[a, pivot]
  104.     ]
  105.             
  106. lusolve[ap_, bb_] :=
  107.     Module[{a, b = bb, pivot, ii, i, m, n, tmp},
  108.     a = ap[[1]];
  109.     pivot = ap[[2]];
  110.     n = Length[a];
  111.     (* forward elimination *)
  112.     For[ii=1, ii<=n, ii++,  (* for each row do ... *)
  113.         For[i=ii+1, i<=n, i++,
  114.         b[[pivot[[i]]]] -= a[[pivot[[i]], ii]]*b[[pivot[[ii]]]]
  115.         ];
  116.         ];
  117.     (* back substitution *)
  118.     Clear[i];
  119.     For[ii=n, ii>0, ii--,
  120.         tmp = Sum[a[[pivot[[ii]], i]]*b[[pivot[[i]]]], {i, ii+1, n}];
  121.         tmp = b[[pivot[[ii]]]] - tmp;
  122.         b[[pivot[[ii]]]] = divide[tmp, a[[pivot[[ii]], ii]]]
  123.         ];
  124.     tmp = Table[,{n}];
  125.     For[ii=1, ii<=n, ii++, tmp[[ii]] = b[[pivot[[ii]]]]];
  126.     tmp
  127.     ]
  128.  
  129. End[ ] (* "LinearAlgebra`GaussianElimination`Private`" *)
  130.  
  131. Protect[LUFactor, LUSolve, LU];
  132.  
  133. EndPackage[ ] (* "LinearAlgebra`GaussianElimination`" *)
  134.  
  135. (*
  136. test3 := Module[{a, x, b, r, i, j}, 
  137.     a = Table[N[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
  138.     b = Table[N[Random[Integer, {-10, 10}]], {i, 3}]; 
  139.     x = LUSolve[LUFactor[a], b];
  140.     r = a . x - b;
  141.     {a, b, x, r}
  142.     ]
  143.  
  144. test5 := Module[{a, b, i, j}, 
  145.     a = Table[N[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
  146.     b = Table[N[Random[Integer, {-10, 10}]], {i, 5}]; 
  147.     a . LUSolve[LUFactor[a], b] - b
  148.     ]
  149.  
  150. cntest3 := Module[{a, x, b, r, i, j}, 
  151.     a = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
  152.     b = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 3}]; 
  153.     x = LUSolve[LUFactor[a], b];
  154.     r = a . x - b;
  155.     {a, b, x, r}
  156.     ]
  157.  
  158. cntest5 := Module[{a, b, i, j}, 
  159.     a = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
  160.     b = Table[ComputerNumber[Random[Integer, {-10, 10}]], {i, 5}]; 
  161.     a . LUSolve[LUFactor[a], b] - b
  162.     ]
  163.  
  164. inttest3 := Module[{a, x, b, r, i, j}, 
  165.     a = Table[Interval[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]; 
  166.     b = Table[Interval[Random[Integer, {-10, 10}]], {i, 3}]; 
  167.     x = LUSolve[LUFactor[a], b];
  168.     r = a . x - b;
  169.     {a, b, x, r}
  170.     ]
  171.  
  172. inttest5 := Module[{a, b, i, j}, 
  173.     a = Table[Interval[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]; 
  174.     b = Table[Interval[Random[Integer, {-10, 10}]], {i, 5}]; 
  175.     a . LUSolve[LUFactor[a], b] - b
  176.     ]
  177.  
  178. ninttest3 := Module[{a, x, b, r, i, j}, 
  179.     a = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 3}, {j, 3}]];
  180.     b = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 3}]];
  181.     x = LUSolve[LUFactor[a], b];
  182.     r = a . x - b;
  183.     {a, b, x, r}
  184.     ]
  185.  
  186. ninttest5 := Module[{a, b, i, j}, 
  187.     a = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 5}, {j, 5}]];
  188.     b = N[Table[Interval[Random[Integer, {-10, 10}]], {i, 5}]];
  189.     a . LUSolve[LUFactor[a], b] - b
  190.     ]
  191. *)
  192.