Module | LUSolve
In: | bigdecimal/lib/bigdecimal/ludcmp.rb

ludcmp.rb

### Methods

ludecomp lusolve

### Public Instance methods

ludecomp(a,n,zero=0,one=1)

[Source]

    # File bigdecimal/lib/bigdecimal/ludcmp.rb, line 5
      def ludecomp(a,n,zero=0,one=1)
        prec = BigDecimal.limit(nil)
        ps     = []
        scales = []
        for i in 0...n do  # pick up largest(abs. val.) element in each row.
          ps <<= i
          nrmrow  = zero
          ixn = i*n
          for j in 0...n do
             biggst = a[ixn+j].abs
             nrmrow = biggst if biggst>nrmrow
          end
          if nrmrow>zero then
             scales <<= one.div(nrmrow,prec)
          else
             raise "Singular matrix"
          end
        end
        n1          = n - 1
        for k in 0...n1 do # Gaussian elimination with partial pivoting.
          biggst  = zero;
          for i in k...n do
             size = a[ps[i]*n+k].abs*scales[ps[i]]
             if size>biggst then
                biggst = size
                pividx  = i
             end
          end
          raise "Singular matrix" if biggst<=zero
          if pividx!=k then
            j = ps[k]
            ps[k] = ps[pividx]
            ps[pividx] = j
          end
          pivot   = a[ps[k]*n+k]
          for i in (k+1)...n do
            psin = ps[i]*n
            a[psin+k] = mult = a[psin+k].div(pivot,prec)
            if mult!=zero then
               pskn = ps[k]*n
               for j in (k+1)...n do
                 a[psin+j] -= mult.mult(a[pskn+j],prec)
               end
            end
          end
        end
        raise "Singular matrix" if a[ps[n1]*n+n1] == zero
        ps
      end

      def lusolve(a,b,ps,zero=0.0)
        prec = BigDecimal.limit(nil)
        n = ps.size
        x = []
        for i in 0...n do
          dot = zero
          psin = ps[i]*n
          for j in 0...i do
            dot = a[psin+j].mult(x[j],prec) + dot
          end
          x <<= b[ps[i]] - dot
        end
        (n-1).downto(0) do |i|
           dot = zero
           psin = ps[i]*n
           for j in (i+1)...n do
             dot = a[psin+j].mult(x[j],prec) + dot
           end
           x[i]  = (x[i]-dot).div(a[psin+i],prec)
        end
        x
      end
    end

lusolve(a,b,ps,zero=0.0)

[Source]

    # File bigdecimal/lib/bigdecimal/ludcmp.rb, line 55
      def lusolve(a,b,ps,zero=0.0)
        prec = BigDecimal.limit(nil)
        n = ps.size
        x = []
        for i in 0...n do
          dot = zero
          psin = ps[i]*n
          for j in 0...i do
            dot = a[psin+j].mult(x[j],prec) + dot
          end
          x <<= b[ps[i]] - dot
        end
        (n-1).downto(0) do |i|
           dot = zero
           psin = ps[i]*n
           for j in (i+1)...n do
             dot = a[psin+j].mult(x[j],prec) + dot
           end
           x[i]  = (x[i]-dot).div(a[psin+i],prec)
        end
        x
      end
    end

[Validate]