diff options
Diffstat (limited to 'lib/matrix')
-rw-r--r-- | lib/matrix/eigenvalue_decomposition.rb | 882 | ||||
-rw-r--r-- | lib/matrix/lup_decomposition.rb | 219 | ||||
-rw-r--r-- | lib/matrix/matrix.gemspec | 26 | ||||
-rw-r--r-- | lib/matrix/version.rb | 5 |
4 files changed, 0 insertions, 1132 deletions
diff --git a/lib/matrix/eigenvalue_decomposition.rb b/lib/matrix/eigenvalue_decomposition.rb deleted file mode 100644 index bf6637635a..0000000000 --- a/lib/matrix/eigenvalue_decomposition.rb +++ /dev/null @@ -1,882 +0,0 @@ -# frozen_string_literal: false -class Matrix - # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/ - - # Eigenvalues and eigenvectors of a real matrix. - # - # Computes the eigenvalues and eigenvectors of a matrix A. - # - # If A is diagonalizable, this provides matrices V and D - # such that A = V*D*V.inv, where D is the diagonal matrix with entries - # equal to the eigenvalues and V is formed by the eigenvectors. - # - # If A is symmetric, then V is orthogonal and thus A = V*D*V.t - - class EigenvalueDecomposition - - # Constructs the eigenvalue decomposition for a square matrix +A+ - # - def initialize(a) - # @d, @e: Arrays for internal storage of eigenvalues. - # @v: Array for internal storage of eigenvectors. - # @h: Array for internal storage of nonsymmetric Hessenberg form. - raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) - @size = a.row_count - @d = Array.new(@size, 0) - @e = Array.new(@size, 0) - - if (@symmetric = a.symmetric?) - @v = a.to_a - tridiagonalize - diagonalize - else - @v = Array.new(@size) { Array.new(@size, 0) } - @h = a.to_a - @ort = Array.new(@size, 0) - reduce_to_hessenberg - hessenberg_to_real_schur - end - end - - # Returns the eigenvector matrix +V+ - # - def eigenvector_matrix - Matrix.send(:new, build_eigenvectors.transpose) - end - alias_method :v, :eigenvector_matrix - - # Returns the inverse of the eigenvector matrix +V+ - # - def eigenvector_matrix_inv - r = Matrix.send(:new, build_eigenvectors) - r = r.transpose.inverse unless @symmetric - r - end - alias_method :v_inv, :eigenvector_matrix_inv - - # Returns the eigenvalues in an array - # - def eigenvalues - values = @d.dup - @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0} - values - end - - # Returns an array of the eigenvectors - # - def eigenvectors - build_eigenvectors.map{|ev| Vector.send(:new, ev)} - end - - # Returns the block diagonal eigenvalue matrix +D+ - # - def eigenvalue_matrix - Matrix.diagonal(*eigenvalues) - end - alias_method :d, :eigenvalue_matrix - - # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv] - # - def to_ary - [v, d, v_inv] - end - alias_method :to_a, :to_ary - - - private def build_eigenvectors - # JAMA stores complex eigenvectors in a strange way - # See http://web.archive.org/web/20111016032731/http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html - @e.each_with_index.map do |imag, i| - if imag == 0 - Array.new(@size){|j| @v[j][i]} - elsif imag > 0 - Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])} - else - Array.new(@size){|j| Complex(@v[j][i-1], -@v[j][i])} - end - end - end - - # Complex scalar division. - - private def cdiv(xr, xi, yr, yi) - if (yr.abs > yi.abs) - r = yi/yr - d = yr + r*yi - [(xr + r*xi)/d, (xi - r*xr)/d] - else - r = yr/yi - d = yi + r*yr - [(r*xr + xi)/d, (r*xi - xr)/d] - end - end - - - # Symmetric Householder reduction to tridiagonal form. - - private def tridiagonalize - - # This is derived from the Algol procedures tred2 by - # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - # Fortran subroutine in EISPACK. - - @size.times do |j| - @d[j] = @v[@size-1][j] - end - - # Householder reduction to tridiagonal form. - - (@size-1).downto(0+1) do |i| - - # Scale to avoid under/overflow. - - scale = 0.0 - h = 0.0 - i.times do |k| - scale = scale + @d[k].abs - end - if (scale == 0.0) - @e[i] = @d[i-1] - i.times do |j| - @d[j] = @v[i-1][j] - @v[i][j] = 0.0 - @v[j][i] = 0.0 - end - else - - # Generate Householder vector. - - i.times do |k| - @d[k] /= scale - h += @d[k] * @d[k] - end - f = @d[i-1] - g = Math.sqrt(h) - if (f > 0) - g = -g - end - @e[i] = scale * g - h -= f * g - @d[i-1] = f - g - i.times do |j| - @e[j] = 0.0 - end - - # Apply similarity transformation to remaining columns. - - i.times do |j| - f = @d[j] - @v[j][i] = f - g = @e[j] + @v[j][j] * f - (j+1).upto(i-1) do |k| - g += @v[k][j] * @d[k] - @e[k] += @v[k][j] * f - end - @e[j] = g - end - f = 0.0 - i.times do |j| - @e[j] /= h - f += @e[j] * @d[j] - end - hh = f / (h + h) - i.times do |j| - @e[j] -= hh * @d[j] - end - i.times do |j| - f = @d[j] - g = @e[j] - j.upto(i-1) do |k| - @v[k][j] -= (f * @e[k] + g * @d[k]) - end - @d[j] = @v[i-1][j] - @v[i][j] = 0.0 - end - end - @d[i] = h - end - - # Accumulate transformations. - - 0.upto(@size-1-1) do |i| - @v[@size-1][i] = @v[i][i] - @v[i][i] = 1.0 - h = @d[i+1] - if (h != 0.0) - 0.upto(i) do |k| - @d[k] = @v[k][i+1] / h - end - 0.upto(i) do |j| - g = 0.0 - 0.upto(i) do |k| - g += @v[k][i+1] * @v[k][j] - end - 0.upto(i) do |k| - @v[k][j] -= g * @d[k] - end - end - end - 0.upto(i) do |k| - @v[k][i+1] = 0.0 - end - end - @size.times do |j| - @d[j] = @v[@size-1][j] - @v[@size-1][j] = 0.0 - end - @v[@size-1][@size-1] = 1.0 - @e[0] = 0.0 - end - - - # Symmetric tridiagonal QL algorithm. - - private def diagonalize - # This is derived from the Algol procedures tql2, by - # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for - # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding - # Fortran subroutine in EISPACK. - - 1.upto(@size-1) do |i| - @e[i-1] = @e[i] - end - @e[@size-1] = 0.0 - - f = 0.0 - tst1 = 0.0 - eps = Float::EPSILON - @size.times do |l| - - # Find small subdiagonal element - - tst1 = [tst1, @d[l].abs + @e[l].abs].max - m = l - while (m < @size) do - if (@e[m].abs <= eps*tst1) - break - end - m+=1 - end - - # If m == l, @d[l] is an eigenvalue, - # otherwise, iterate. - - if (m > l) - iter = 0 - begin - iter = iter + 1 # (Could check iteration count here.) - - # Compute implicit shift - - g = @d[l] - p = (@d[l+1] - g) / (2.0 * @e[l]) - r = Math.hypot(p, 1.0) - if (p < 0) - r = -r - end - @d[l] = @e[l] / (p + r) - @d[l+1] = @e[l] * (p + r) - dl1 = @d[l+1] - h = g - @d[l] - (l+2).upto(@size-1) do |i| - @d[i] -= h - end - f += h - - # Implicit QL transformation. - - p = @d[m] - c = 1.0 - c2 = c - c3 = c - el1 = @e[l+1] - s = 0.0 - s2 = 0.0 - (m-1).downto(l) do |i| - c3 = c2 - c2 = c - s2 = s - g = c * @e[i] - h = c * p - r = Math.hypot(p, @e[i]) - @e[i+1] = s * r - s = @e[i] / r - c = p / r - p = c * @d[i] - s * g - @d[i+1] = h + s * (c * g + s * @d[i]) - - # Accumulate transformation. - - @size.times do |k| - h = @v[k][i+1] - @v[k][i+1] = s * @v[k][i] + c * h - @v[k][i] = c * @v[k][i] - s * h - end - end - p = -s * s2 * c3 * el1 * @e[l] / dl1 - @e[l] = s * p - @d[l] = c * p - - # Check for convergence. - - end while (@e[l].abs > eps*tst1) - end - @d[l] = @d[l] + f - @e[l] = 0.0 - end - - # Sort eigenvalues and corresponding vectors. - - 0.upto(@size-2) do |i| - k = i - p = @d[i] - (i+1).upto(@size-1) do |j| - if (@d[j] < p) - k = j - p = @d[j] - end - end - if (k != i) - @d[k] = @d[i] - @d[i] = p - @size.times do |j| - p = @v[j][i] - @v[j][i] = @v[j][k] - @v[j][k] = p - end - end - end - end - - # Nonsymmetric reduction to Hessenberg form. - - private def reduce_to_hessenberg - # This is derived from the Algol procedures orthes and ortran, - # by Martin and Wilkinson, Handbook for Auto. Comp., - # Vol.ii-Linear Algebra, and the corresponding - # Fortran subroutines in EISPACK. - - low = 0 - high = @size-1 - - (low+1).upto(high-1) do |m| - - # Scale column. - - scale = 0.0 - m.upto(high) do |i| - scale = scale + @h[i][m-1].abs - end - if (scale != 0.0) - - # Compute Householder transformation. - - h = 0.0 - high.downto(m) do |i| - @ort[i] = @h[i][m-1]/scale - h += @ort[i] * @ort[i] - end - g = Math.sqrt(h) - if (@ort[m] > 0) - g = -g - end - h -= @ort[m] * g - @ort[m] = @ort[m] - g - - # Apply Householder similarity transformation - # @h = (I-u*u'/h)*@h*(I-u*u')/h) - - m.upto(@size-1) do |j| - f = 0.0 - high.downto(m) do |i| - f += @ort[i]*@h[i][j] - end - f = f/h - m.upto(high) do |i| - @h[i][j] -= f*@ort[i] - end - end - - 0.upto(high) do |i| - f = 0.0 - high.downto(m) do |j| - f += @ort[j]*@h[i][j] - end - f = f/h - m.upto(high) do |j| - @h[i][j] -= f*@ort[j] - end - end - @ort[m] = scale*@ort[m] - @h[m][m-1] = scale*g - end - end - - # Accumulate transformations (Algol's ortran). - - @size.times do |i| - @size.times do |j| - @v[i][j] = (i == j ? 1.0 : 0.0) - end - end - - (high-1).downto(low+1) do |m| - if (@h[m][m-1] != 0.0) - (m+1).upto(high) do |i| - @ort[i] = @h[i][m-1] - end - m.upto(high) do |j| - g = 0.0 - m.upto(high) do |i| - g += @ort[i] * @v[i][j] - end - # Double division avoids possible underflow - g = (g / @ort[m]) / @h[m][m-1] - m.upto(high) do |i| - @v[i][j] += g * @ort[i] - end - end - end - end - end - - # Nonsymmetric reduction from Hessenberg to real Schur form. - - private def hessenberg_to_real_schur - - # This is derived from the Algol procedure hqr2, - # by Martin and Wilkinson, Handbook for Auto. Comp., - # Vol.ii-Linear Algebra, and the corresponding - # Fortran subroutine in EISPACK. - - # Initialize - - nn = @size - n = nn-1 - low = 0 - high = nn-1 - eps = Float::EPSILON - exshift = 0.0 - p = q = r = s = z = 0 - - # Store roots isolated by balanc and compute matrix norm - - norm = 0.0 - nn.times do |i| - if (i < low || i > high) - @d[i] = @h[i][i] - @e[i] = 0.0 - end - ([i-1, 0].max).upto(nn-1) do |j| - norm = norm + @h[i][j].abs - end - end - - # Outer loop over eigenvalue index - - iter = 0 - while (n >= low) do - - # Look for single small sub-diagonal element - - l = n - while (l > low) do - s = @h[l-1][l-1].abs + @h[l][l].abs - if (s == 0.0) - s = norm - end - if (@h[l][l-1].abs < eps * s) - break - end - l-=1 - end - - # Check for convergence - # One root found - - if (l == n) - @h[n][n] = @h[n][n] + exshift - @d[n] = @h[n][n] - @e[n] = 0.0 - n-=1 - iter = 0 - - # Two roots found - - elsif (l == n-1) - w = @h[n][n-1] * @h[n-1][n] - p = (@h[n-1][n-1] - @h[n][n]) / 2.0 - q = p * p + w - z = Math.sqrt(q.abs) - @h[n][n] = @h[n][n] + exshift - @h[n-1][n-1] = @h[n-1][n-1] + exshift - x = @h[n][n] - - # Real pair - - if (q >= 0) - if (p >= 0) - z = p + z - else - z = p - z - end - @d[n-1] = x + z - @d[n] = @d[n-1] - if (z != 0.0) - @d[n] = x - w / z - end - @e[n-1] = 0.0 - @e[n] = 0.0 - x = @h[n][n-1] - s = x.abs + z.abs - p = x / s - q = z / s - r = Math.sqrt(p * p+q * q) - p /= r - q /= r - - # Row modification - - (n-1).upto(nn-1) do |j| - z = @h[n-1][j] - @h[n-1][j] = q * z + p * @h[n][j] - @h[n][j] = q * @h[n][j] - p * z - end - - # Column modification - - 0.upto(n) do |i| - z = @h[i][n-1] - @h[i][n-1] = q * z + p * @h[i][n] - @h[i][n] = q * @h[i][n] - p * z - end - - # Accumulate transformations - - low.upto(high) do |i| - z = @v[i][n-1] - @v[i][n-1] = q * z + p * @v[i][n] - @v[i][n] = q * @v[i][n] - p * z - end - - # Complex pair - - else - @d[n-1] = x + p - @d[n] = x + p - @e[n-1] = z - @e[n] = -z - end - n -= 2 - iter = 0 - - # No convergence yet - - else - - # Form shift - - x = @h[n][n] - y = 0.0 - w = 0.0 - if (l < n) - y = @h[n-1][n-1] - w = @h[n][n-1] * @h[n-1][n] - end - - # Wilkinson's original ad hoc shift - - if (iter == 10) - exshift += x - low.upto(n) do |i| - @h[i][i] -= x - end - s = @h[n][n-1].abs + @h[n-1][n-2].abs - x = y = 0.75 * s - w = -0.4375 * s * s - end - - # MATLAB's new ad hoc shift - - if (iter == 30) - s = (y - x) / 2.0 - s *= s + w - if (s > 0) - s = Math.sqrt(s) - if (y < x) - s = -s - end - s = x - w / ((y - x) / 2.0 + s) - low.upto(n) do |i| - @h[i][i] -= s - end - exshift += s - x = y = w = 0.964 - end - end - - iter = iter + 1 # (Could check iteration count here.) - - # Look for two consecutive small sub-diagonal elements - - m = n-2 - while (m >= l) do - z = @h[m][m] - r = x - z - s = y - z - p = (r * s - w) / @h[m+1][m] + @h[m][m+1] - q = @h[m+1][m+1] - z - r - s - r = @h[m+2][m+1] - s = p.abs + q.abs + r.abs - p /= s - q /= s - r /= s - if (m == l) - break - end - if (@h[m][m-1].abs * (q.abs + r.abs) < - eps * (p.abs * (@h[m-1][m-1].abs + z.abs + - @h[m+1][m+1].abs))) - break - end - m-=1 - end - - (m+2).upto(n) do |i| - @h[i][i-2] = 0.0 - if (i > m+2) - @h[i][i-3] = 0.0 - end - end - - # Double QR step involving rows l:n and columns m:n - - m.upto(n-1) do |k| - notlast = (k != n-1) - if (k != m) - p = @h[k][k-1] - q = @h[k+1][k-1] - r = (notlast ? @h[k+2][k-1] : 0.0) - x = p.abs + q.abs + r.abs - next if x == 0 - p /= x - q /= x - r /= x - end - s = Math.sqrt(p * p + q * q + r * r) - if (p < 0) - s = -s - end - if (s != 0) - if (k != m) - @h[k][k-1] = -s * x - elsif (l != m) - @h[k][k-1] = -@h[k][k-1] - end - p += s - x = p / s - y = q / s - z = r / s - q /= p - r /= p - - # Row modification - - k.upto(nn-1) do |j| - p = @h[k][j] + q * @h[k+1][j] - if (notlast) - p += r * @h[k+2][j] - @h[k+2][j] = @h[k+2][j] - p * z - end - @h[k][j] = @h[k][j] - p * x - @h[k+1][j] = @h[k+1][j] - p * y - end - - # Column modification - - 0.upto([n, k+3].min) do |i| - p = x * @h[i][k] + y * @h[i][k+1] - if (notlast) - p += z * @h[i][k+2] - @h[i][k+2] = @h[i][k+2] - p * r - end - @h[i][k] = @h[i][k] - p - @h[i][k+1] = @h[i][k+1] - p * q - end - - # Accumulate transformations - - low.upto(high) do |i| - p = x * @v[i][k] + y * @v[i][k+1] - if (notlast) - p += z * @v[i][k+2] - @v[i][k+2] = @v[i][k+2] - p * r - end - @v[i][k] = @v[i][k] - p - @v[i][k+1] = @v[i][k+1] - p * q - end - end # (s != 0) - end # k loop - end # check convergence - end # while (n >= low) - - # Backsubstitute to find vectors of upper triangular form - - if (norm == 0.0) - return - end - - (nn-1).downto(0) do |k| - p = @d[k] - q = @e[k] - - # Real vector - - if (q == 0) - l = k - @h[k][k] = 1.0 - (k-1).downto(0) do |i| - w = @h[i][i] - p - r = 0.0 - l.upto(k) do |j| - r += @h[i][j] * @h[j][k] - end - if (@e[i] < 0.0) - z = w - s = r - else - l = i - if (@e[i] == 0.0) - if (w != 0.0) - @h[i][k] = -r / w - else - @h[i][k] = -r / (eps * norm) - end - - # Solve real equations - - else - x = @h[i][i+1] - y = @h[i+1][i] - q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - t = (x * s - z * r) / q - @h[i][k] = t - if (x.abs > z.abs) - @h[i+1][k] = (-r - w * t) / x - else - @h[i+1][k] = (-s - y * t) / z - end - end - - # Overflow control - - t = @h[i][k].abs - if ((eps * t) * t > 1) - i.upto(k) do |j| - @h[j][k] = @h[j][k] / t - end - end - end - end - - # Complex vector - - elsif (q < 0) - l = n-1 - - # Last vector component imaginary so matrix is triangular - - if (@h[n][n-1].abs > @h[n-1][n].abs) - @h[n-1][n-1] = q / @h[n][n-1] - @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1] - else - cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q) - @h[n-1][n-1] = cdivr - @h[n-1][n] = cdivi - end - @h[n][n-1] = 0.0 - @h[n][n] = 1.0 - (n-2).downto(0) do |i| - ra = 0.0 - sa = 0.0 - l.upto(n) do |j| - ra = ra + @h[i][j] * @h[j][n-1] - sa = sa + @h[i][j] * @h[j][n] - end - w = @h[i][i] - p - - if (@e[i] < 0.0) - z = w - r = ra - s = sa - else - l = i - if (@e[i] == 0) - cdivr, cdivi = cdiv(-ra, -sa, w, q) - @h[i][n-1] = cdivr - @h[i][n] = cdivi - else - - # Solve complex equations - - x = @h[i][i+1] - y = @h[i+1][i] - vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q - vi = (@d[i] - p) * 2.0 * q - if (vr == 0.0 && vi == 0.0) - vr = eps * norm * (w.abs + q.abs + - x.abs + y.abs + z.abs) - end - cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi) - @h[i][n-1] = cdivr - @h[i][n] = cdivi - if (x.abs > (z.abs + q.abs)) - @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x - @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x - else - cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q) - @h[i+1][n-1] = cdivr - @h[i+1][n] = cdivi - end - end - - # Overflow control - - t = [@h[i][n-1].abs, @h[i][n].abs].max - if ((eps * t) * t > 1) - i.upto(n) do |j| - @h[j][n-1] = @h[j][n-1] / t - @h[j][n] = @h[j][n] / t - end - end - end - end - end - end - - # Vectors of isolated roots - - nn.times do |i| - if (i < low || i > high) - i.upto(nn-1) do |j| - @v[i][j] = @h[i][j] - end - end - end - - # Back transformation to get eigenvectors of original matrix - - (nn-1).downto(low) do |j| - low.upto(high) do |i| - z = 0.0 - low.upto([j, high].min) do |k| - z += @v[i][k] * @h[k][j] - end - @v[i][j] = z - end - end - end - - end -end diff --git a/lib/matrix/lup_decomposition.rb b/lib/matrix/lup_decomposition.rb deleted file mode 100644 index e37def75f6..0000000000 --- a/lib/matrix/lup_decomposition.rb +++ /dev/null @@ -1,219 +0,0 @@ -# frozen_string_literal: false -class Matrix - # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/ - - # - # For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n - # unit lower triangular matrix L, an n-by-n upper triangular matrix U, - # and a m-by-m permutation matrix P so that L*U = P*A. - # If m < n, then L is m-by-m and U is m-by-n. - # - # The LUP decomposition with pivoting always exists, even if the matrix is - # singular, so the constructor will never fail. The primary use of the - # LU decomposition is in the solution of square systems of simultaneous - # linear equations. This will fail if singular? returns true. - # - - class LUPDecomposition - # Returns the lower triangular factor +L+ - - include Matrix::ConversionHelper - - def l - Matrix.build(@row_count, [@column_count, @row_count].min) do |i, j| - if (i > j) - @lu[i][j] - elsif (i == j) - 1 - else - 0 - end - end - end - - # Returns the upper triangular factor +U+ - - def u - Matrix.build([@column_count, @row_count].min, @column_count) do |i, j| - if (i <= j) - @lu[i][j] - else - 0 - end - end - end - - # Returns the permutation matrix +P+ - - def p - rows = Array.new(@row_count){Array.new(@row_count, 0)} - @pivots.each_with_index{|p, i| rows[i][p] = 1} - Matrix.send :new, rows, @row_count - end - - # Returns +L+, +U+, +P+ in an array - - def to_ary - [l, u, p] - end - alias_method :to_a, :to_ary - - # Returns the pivoting indices - - attr_reader :pivots - - # Returns +true+ if +U+, and hence +A+, is singular. - - def singular? - @column_count.times do |j| - if (@lu[j][j] == 0) - return true - end - end - false - end - - # Returns the determinant of +A+, calculated efficiently - # from the factorization. - - def det - if (@row_count != @column_count) - raise Matrix::ErrDimensionMismatch - end - d = @pivot_sign - @column_count.times do |j| - d *= @lu[j][j] - end - d - end - alias_method :determinant, :det - - # Returns +m+ so that <tt>A*m = b</tt>, - # or equivalently so that <tt>L*U*m = P*b</tt> - # +b+ can be a Matrix or a Vector - - def solve b - if (singular?) - raise Matrix::ErrNotRegular, "Matrix is singular." - end - if b.is_a? Matrix - if (b.row_count != @row_count) - raise Matrix::ErrDimensionMismatch - end - - # Copy right hand side with pivoting - nx = b.column_count - m = @pivots.map{|row| b.row(row).to_a} - - # Solve L*Y = P*b - @column_count.times do |k| - (k+1).upto(@column_count-1) do |i| - nx.times do |j| - m[i][j] -= m[k][j]*@lu[i][k] - end - end - end - # Solve U*m = Y - (@column_count-1).downto(0) do |k| - nx.times do |j| - m[k][j] = m[k][j].quo(@lu[k][k]) - end - k.times do |i| - nx.times do |j| - m[i][j] -= m[k][j]*@lu[i][k] - end - end - end - Matrix.send :new, m, nx - else # same algorithm, specialized for simpler case of a vector - b = convert_to_array(b) - if (b.size != @row_count) - raise Matrix::ErrDimensionMismatch - end - - # Copy right hand side with pivoting - m = b.values_at(*@pivots) - - # Solve L*Y = P*b - @column_count.times do |k| - (k+1).upto(@column_count-1) do |i| - m[i] -= m[k]*@lu[i][k] - end - end - # Solve U*m = Y - (@column_count-1).downto(0) do |k| - m[k] = m[k].quo(@lu[k][k]) - k.times do |i| - m[i] -= m[k]*@lu[i][k] - end - end - Vector.elements(m, false) - end - end - - def initialize a - raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) - # Use a "left-looking", dot-product, Crout/Doolittle algorithm. - @lu = a.to_a - @row_count = a.row_count - @column_count = a.column_count - @pivots = Array.new(@row_count) - @row_count.times do |i| - @pivots[i] = i - end - @pivot_sign = 1 - lu_col_j = Array.new(@row_count) - - # Outer loop. - - @column_count.times do |j| - - # Make a copy of the j-th column to localize references. - - @row_count.times do |i| - lu_col_j[i] = @lu[i][j] - end - - # Apply previous transformations. - - @row_count.times do |i| - lu_row_i = @lu[i] - - # Most of the time is spent in the following dot product. - - kmax = [i, j].min - s = 0 - kmax.times do |k| - s += lu_row_i[k]*lu_col_j[k] - end - - lu_row_i[j] = lu_col_j[i] -= s - end - - # Find pivot and exchange if necessary. - - p = j - (j+1).upto(@row_count-1) do |i| - if (lu_col_j[i].abs > lu_col_j[p].abs) - p = i - end - end - if (p != j) - @column_count.times do |k| - t = @lu[p][k]; @lu[p][k] = @lu[j][k]; @lu[j][k] = t - end - k = @pivots[p]; @pivots[p] = @pivots[j]; @pivots[j] = k - @pivot_sign = -@pivot_sign - end - - # Compute multipliers. - - if (j < @row_count && @lu[j][j] != 0) - (j+1).upto(@row_count-1) do |i| - @lu[i][j] = @lu[i][j].quo(@lu[j][j]) - end - end - end - end - end -end diff --git a/lib/matrix/matrix.gemspec b/lib/matrix/matrix.gemspec deleted file mode 100644 index 6f68cbd816..0000000000 --- a/lib/matrix/matrix.gemspec +++ /dev/null @@ -1,26 +0,0 @@ -# frozen_string_literal: true - -begin - require_relative "lib/matrix/version" -rescue LoadError - # for Ruby core repository - require_relative "version" -end - -Gem::Specification.new do |spec| - spec.name = "matrix" - spec.version = Matrix::VERSION - spec.authors = ["Marc-Andre Lafortune"] - spec.email = ["ruby-core@marc-andre.ca"] - - spec.summary = %q{An implementation of Matrix and Vector classes.} - spec.description = %q{An implementation of Matrix and Vector classes.} - spec.homepage = "https://github.com/ruby/matrix" - spec.licenses = ["Ruby", "BSD-2-Clause"] - spec.required_ruby_version = ">= 2.5.0" - - spec.files = [".gitignore", "Gemfile", "LICENSE.txt", "README.md", "Rakefile", "bin/console", "bin/setup", "lib/matrix.rb", "lib/matrix/eigenvalue_decomposition.rb", "lib/matrix/lup_decomposition.rb", "lib/matrix/version.rb", "matrix.gemspec"] - spec.bindir = "exe" - spec.executables = spec.files.grep(%r{^exe/}) { |f| File.basename(f) } - spec.require_paths = ["lib"] -end diff --git a/lib/matrix/version.rb b/lib/matrix/version.rb deleted file mode 100644 index 82b835f038..0000000000 --- a/lib/matrix/version.rb +++ /dev/null @@ -1,5 +0,0 @@ -# frozen_string_literal: true - -class Matrix - VERSION = "0.4.1" -end |