/tools/Ruby/lib/ruby/1.8/bigdecimal/ludcmp.rb

http://github.com/agross/netopenspace · Ruby · 84 lines · 73 code · 1 blank · 10 comment · 18 complexity · 9ed34f0d05b908c9554ce9484e182ac3 MD5 · raw file

  1. #
  2. # Solves a*x = b for x, using LU decomposition.
  3. #
  4. module LUSolve
  5. # Performs LU decomposition of the n by n matrix a.
  6. def ludecomp(a,n,zero=0,one=1)
  7. prec = BigDecimal.limit(nil)
  8. ps = []
  9. scales = []
  10. for i in 0...n do # pick up largest(abs. val.) element in each row.
  11. ps <<= i
  12. nrmrow = zero
  13. ixn = i*n
  14. for j in 0...n do
  15. biggst = a[ixn+j].abs
  16. nrmrow = biggst if biggst>nrmrow
  17. end
  18. if nrmrow>zero then
  19. scales <<= one.div(nrmrow,prec)
  20. else
  21. raise "Singular matrix"
  22. end
  23. end
  24. n1 = n - 1
  25. for k in 0...n1 do # Gaussian elimination with partial pivoting.
  26. biggst = zero;
  27. for i in k...n do
  28. size = a[ps[i]*n+k].abs*scales[ps[i]]
  29. if size>biggst then
  30. biggst = size
  31. pividx = i
  32. end
  33. end
  34. raise "Singular matrix" if biggst<=zero
  35. if pividx!=k then
  36. j = ps[k]
  37. ps[k] = ps[pividx]
  38. ps[pividx] = j
  39. end
  40. pivot = a[ps[k]*n+k]
  41. for i in (k+1)...n do
  42. psin = ps[i]*n
  43. a[psin+k] = mult = a[psin+k].div(pivot,prec)
  44. if mult!=zero then
  45. pskn = ps[k]*n
  46. for j in (k+1)...n do
  47. a[psin+j] -= mult.mult(a[pskn+j],prec)
  48. end
  49. end
  50. end
  51. end
  52. raise "Singular matrix" if a[ps[n1]*n+n1] == zero
  53. ps
  54. end
  55. # Solves a*x = b for x, using LU decomposition.
  56. #
  57. # a is a matrix, b is a constant vector, x is the solution vector.
  58. #
  59. # ps is the pivot, a vector which indicates the permutation of rows performed
  60. # during LU decomposition.
  61. def lusolve(a,b,ps,zero=0.0)
  62. prec = BigDecimal.limit(nil)
  63. n = ps.size
  64. x = []
  65. for i in 0...n do
  66. dot = zero
  67. psin = ps[i]*n
  68. for j in 0...i do
  69. dot = a[psin+j].mult(x[j],prec) + dot
  70. end
  71. x <<= b[ps[i]] - dot
  72. end
  73. (n-1).downto(0) do |i|
  74. dot = zero
  75. psin = ps[i]*n
  76. for j in (i+1)...n do
  77. dot = a[psin+j].mult(x[j],prec) + dot
  78. end
  79. x[i] = (x[i]-dot).div(a[psin+i],prec)
  80. end
  81. x
  82. end
  83. end