1class Matrix 2 # Adapted from JAMA: http://math.nist.gov/javanumerics/jama/ 3 4 # Eigenvalues and eigenvectors of a real matrix. 5 # 6 # Computes the eigenvalues and eigenvectors of a matrix A. 7 # 8 # If A is diagonalizable, this provides matrices V and D 9 # such that A = V*D*V.inv, where D is the diagonal matrix with entries 10 # equal to the eigenvalues and V is formed by the eigenvectors. 11 # 12 # If A is symmetric, then V is orthogonal and thus A = V*D*V.t 13 14 class EigenvalueDecomposition 15 16 # Constructs the eigenvalue decomposition for a square matrix +A+ 17 # 18 def initialize(a) 19 # @d, @e: Arrays for internal storage of eigenvalues. 20 # @v: Array for internal storage of eigenvectors. 21 # @h: Array for internal storage of nonsymmetric Hessenberg form. 22 raise TypeError, "Expected Matrix but got #{a.class}" unless a.is_a?(Matrix) 23 @size = a.row_count 24 @d = Array.new(@size, 0) 25 @e = Array.new(@size, 0) 26 27 if (@symmetric = a.symmetric?) 28 @v = a.to_a 29 tridiagonalize 30 diagonalize 31 else 32 @v = Array.new(@size) { Array.new(@size, 0) } 33 @h = a.to_a 34 @ort = Array.new(@size, 0) 35 reduce_to_hessenberg 36 hessenberg_to_real_schur 37 end 38 end 39 40 # Returns the eigenvector matrix +V+ 41 # 42 def eigenvector_matrix 43 Matrix.send :new, build_eigenvectors.transpose 44 end 45 alias v eigenvector_matrix 46 47 # Returns the inverse of the eigenvector matrix +V+ 48 # 49 def eigenvector_matrix_inv 50 r = Matrix.send :new, build_eigenvectors 51 r = r.transpose.inverse unless @symmetric 52 r 53 end 54 alias v_inv eigenvector_matrix_inv 55 56 # Returns the eigenvalues in an array 57 # 58 def eigenvalues 59 values = @d.dup 60 @e.each_with_index{|imag, i| values[i] = Complex(values[i], imag) unless imag == 0} 61 values 62 end 63 64 # Returns an array of the eigenvectors 65 # 66 def eigenvectors 67 build_eigenvectors.map{|ev| Vector.send :new, ev} 68 end 69 70 # Returns the block diagonal eigenvalue matrix +D+ 71 # 72 def eigenvalue_matrix 73 Matrix.diagonal(*eigenvalues) 74 end 75 alias d eigenvalue_matrix 76 77 # Returns [eigenvector_matrix, eigenvalue_matrix, eigenvector_matrix_inv] 78 # 79 def to_ary 80 [v, d, v_inv] 81 end 82 alias_method :to_a, :to_ary 83 84 private 85 def build_eigenvectors 86 # JAMA stores complex eigenvectors in a strange way 87 # See http://web.archive.org/web/20111016032731/http://cio.nist.gov/esd/emaildir/lists/jama/msg01021.html 88 @e.each_with_index.map do |imag, i| 89 if imag == 0 90 Array.new(@size){|j| @v[j][i]} 91 elsif imag > 0 92 Array.new(@size){|j| Complex(@v[j][i], @v[j][i+1])} 93 else 94 Array.new(@size){|j| Complex(@v[j][i-1], -@v[j][i])} 95 end 96 end 97 end 98 # Complex scalar division. 99 100 def cdiv(xr, xi, yr, yi) 101 if (yr.abs > yi.abs) 102 r = yi/yr 103 d = yr + r*yi 104 [(xr + r*xi)/d, (xi - r*xr)/d] 105 else 106 r = yr/yi 107 d = yi + r*yr 108 [(r*xr + xi)/d, (r*xi - xr)/d] 109 end 110 end 111 112 113 # Symmetric Householder reduction to tridiagonal form. 114 115 def tridiagonalize 116 117 # This is derived from the Algol procedures tred2 by 118 # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 119 # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 120 # Fortran subroutine in EISPACK. 121 122 @size.times do |j| 123 @d[j] = @v[@size-1][j] 124 end 125 126 # Householder reduction to tridiagonal form. 127 128 (@size-1).downto(0+1) do |i| 129 130 # Scale to avoid under/overflow. 131 132 scale = 0.0 133 h = 0.0 134 i.times do |k| 135 scale = scale + @d[k].abs 136 end 137 if (scale == 0.0) 138 @e[i] = @d[i-1] 139 i.times do |j| 140 @d[j] = @v[i-1][j] 141 @v[i][j] = 0.0 142 @v[j][i] = 0.0 143 end 144 else 145 146 # Generate Householder vector. 147 148 i.times do |k| 149 @d[k] /= scale 150 h += @d[k] * @d[k] 151 end 152 f = @d[i-1] 153 g = Math.sqrt(h) 154 if (f > 0) 155 g = -g 156 end 157 @e[i] = scale * g 158 h -= f * g 159 @d[i-1] = f - g 160 i.times do |j| 161 @e[j] = 0.0 162 end 163 164 # Apply similarity transformation to remaining columns. 165 166 i.times do |j| 167 f = @d[j] 168 @v[j][i] = f 169 g = @e[j] + @v[j][j] * f 170 (j+1).upto(i-1) do |k| 171 g += @v[k][j] * @d[k] 172 @e[k] += @v[k][j] * f 173 end 174 @e[j] = g 175 end 176 f = 0.0 177 i.times do |j| 178 @e[j] /= h 179 f += @e[j] * @d[j] 180 end 181 hh = f / (h + h) 182 i.times do |j| 183 @e[j] -= hh * @d[j] 184 end 185 i.times do |j| 186 f = @d[j] 187 g = @e[j] 188 j.upto(i-1) do |k| 189 @v[k][j] -= (f * @e[k] + g * @d[k]) 190 end 191 @d[j] = @v[i-1][j] 192 @v[i][j] = 0.0 193 end 194 end 195 @d[i] = h 196 end 197 198 # Accumulate transformations. 199 200 0.upto(@size-1-1) do |i| 201 @v[@size-1][i] = @v[i][i] 202 @v[i][i] = 1.0 203 h = @d[i+1] 204 if (h != 0.0) 205 0.upto(i) do |k| 206 @d[k] = @v[k][i+1] / h 207 end 208 0.upto(i) do |j| 209 g = 0.0 210 0.upto(i) do |k| 211 g += @v[k][i+1] * @v[k][j] 212 end 213 0.upto(i) do |k| 214 @v[k][j] -= g * @d[k] 215 end 216 end 217 end 218 0.upto(i) do |k| 219 @v[k][i+1] = 0.0 220 end 221 end 222 @size.times do |j| 223 @d[j] = @v[@size-1][j] 224 @v[@size-1][j] = 0.0 225 end 226 @v[@size-1][@size-1] = 1.0 227 @e[0] = 0.0 228 end 229 230 231 # Symmetric tridiagonal QL algorithm. 232 233 def diagonalize 234 # This is derived from the Algol procedures tql2, by 235 # Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 236 # Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 237 # Fortran subroutine in EISPACK. 238 239 1.upto(@size-1) do |i| 240 @e[i-1] = @e[i] 241 end 242 @e[@size-1] = 0.0 243 244 f = 0.0 245 tst1 = 0.0 246 eps = Float::EPSILON 247 @size.times do |l| 248 249 # Find small subdiagonal element 250 251 tst1 = [tst1, @d[l].abs + @e[l].abs].max 252 m = l 253 while (m < @size) do 254 if (@e[m].abs <= eps*tst1) 255 break 256 end 257 m+=1 258 end 259 260 # If m == l, @d[l] is an eigenvalue, 261 # otherwise, iterate. 262 263 if (m > l) 264 iter = 0 265 begin 266 iter = iter + 1 # (Could check iteration count here.) 267 268 # Compute implicit shift 269 270 g = @d[l] 271 p = (@d[l+1] - g) / (2.0 * @e[l]) 272 r = Math.hypot(p, 1.0) 273 if (p < 0) 274 r = -r 275 end 276 @d[l] = @e[l] / (p + r) 277 @d[l+1] = @e[l] * (p + r) 278 dl1 = @d[l+1] 279 h = g - @d[l] 280 (l+2).upto(@size-1) do |i| 281 @d[i] -= h 282 end 283 f += h 284 285 # Implicit QL transformation. 286 287 p = @d[m] 288 c = 1.0 289 c2 = c 290 c3 = c 291 el1 = @e[l+1] 292 s = 0.0 293 s2 = 0.0 294 (m-1).downto(l) do |i| 295 c3 = c2 296 c2 = c 297 s2 = s 298 g = c * @e[i] 299 h = c * p 300 r = Math.hypot(p, @e[i]) 301 @e[i+1] = s * r 302 s = @e[i] / r 303 c = p / r 304 p = c * @d[i] - s * g 305 @d[i+1] = h + s * (c * g + s * @d[i]) 306 307 # Accumulate transformation. 308 309 @size.times do |k| 310 h = @v[k][i+1] 311 @v[k][i+1] = s * @v[k][i] + c * h 312 @v[k][i] = c * @v[k][i] - s * h 313 end 314 end 315 p = -s * s2 * c3 * el1 * @e[l] / dl1 316 @e[l] = s * p 317 @d[l] = c * p 318 319 # Check for convergence. 320 321 end while (@e[l].abs > eps*tst1) 322 end 323 @d[l] = @d[l] + f 324 @e[l] = 0.0 325 end 326 327 # Sort eigenvalues and corresponding vectors. 328 329 0.upto(@size-2) do |i| 330 k = i 331 p = @d[i] 332 (i+1).upto(@size-1) do |j| 333 if (@d[j] < p) 334 k = j 335 p = @d[j] 336 end 337 end 338 if (k != i) 339 @d[k] = @d[i] 340 @d[i] = p 341 @size.times do |j| 342 p = @v[j][i] 343 @v[j][i] = @v[j][k] 344 @v[j][k] = p 345 end 346 end 347 end 348 end 349 350 # Nonsymmetric reduction to Hessenberg form. 351 352 def reduce_to_hessenberg 353 # This is derived from the Algol procedures orthes and ortran, 354 # by Martin and Wilkinson, Handbook for Auto. Comp., 355 # Vol.ii-Linear Algebra, and the corresponding 356 # Fortran subroutines in EISPACK. 357 358 low = 0 359 high = @size-1 360 361 (low+1).upto(high-1) do |m| 362 363 # Scale column. 364 365 scale = 0.0 366 m.upto(high) do |i| 367 scale = scale + @h[i][m-1].abs 368 end 369 if (scale != 0.0) 370 371 # Compute Householder transformation. 372 373 h = 0.0 374 high.downto(m) do |i| 375 @ort[i] = @h[i][m-1]/scale 376 h += @ort[i] * @ort[i] 377 end 378 g = Math.sqrt(h) 379 if (@ort[m] > 0) 380 g = -g 381 end 382 h -= @ort[m] * g 383 @ort[m] = @ort[m] - g 384 385 # Apply Householder similarity transformation 386 # @h = (I-u*u'/h)*@h*(I-u*u')/h) 387 388 m.upto(@size-1) do |j| 389 f = 0.0 390 high.downto(m) do |i| 391 f += @ort[i]*@h[i][j] 392 end 393 f = f/h 394 m.upto(high) do |i| 395 @h[i][j] -= f*@ort[i] 396 end 397 end 398 399 0.upto(high) do |i| 400 f = 0.0 401 high.downto(m) do |j| 402 f += @ort[j]*@h[i][j] 403 end 404 f = f/h 405 m.upto(high) do |j| 406 @h[i][j] -= f*@ort[j] 407 end 408 end 409 @ort[m] = scale*@ort[m] 410 @h[m][m-1] = scale*g 411 end 412 end 413 414 # Accumulate transformations (Algol's ortran). 415 416 @size.times do |i| 417 @size.times do |j| 418 @v[i][j] = (i == j ? 1.0 : 0.0) 419 end 420 end 421 422 (high-1).downto(low+1) do |m| 423 if (@h[m][m-1] != 0.0) 424 (m+1).upto(high) do |i| 425 @ort[i] = @h[i][m-1] 426 end 427 m.upto(high) do |j| 428 g = 0.0 429 m.upto(high) do |i| 430 g += @ort[i] * @v[i][j] 431 end 432 # Double division avoids possible underflow 433 g = (g / @ort[m]) / @h[m][m-1] 434 m.upto(high) do |i| 435 @v[i][j] += g * @ort[i] 436 end 437 end 438 end 439 end 440 end 441 442 443 444 # Nonsymmetric reduction from Hessenberg to real Schur form. 445 446 def hessenberg_to_real_schur 447 448 # This is derived from the Algol procedure hqr2, 449 # by Martin and Wilkinson, Handbook for Auto. Comp., 450 # Vol.ii-Linear Algebra, and the corresponding 451 # Fortran subroutine in EISPACK. 452 453 # Initialize 454 455 nn = @size 456 n = nn-1 457 low = 0 458 high = nn-1 459 eps = Float::EPSILON 460 exshift = 0.0 461 p=q=r=s=z=0 462 463 # Store roots isolated by balanc and compute matrix norm 464 465 norm = 0.0 466 nn.times do |i| 467 if (i < low || i > high) 468 @d[i] = @h[i][i] 469 @e[i] = 0.0 470 end 471 ([i-1, 0].max).upto(nn-1) do |j| 472 norm = norm + @h[i][j].abs 473 end 474 end 475 476 # Outer loop over eigenvalue index 477 478 iter = 0 479 while (n >= low) do 480 481 # Look for single small sub-diagonal element 482 483 l = n 484 while (l > low) do 485 s = @h[l-1][l-1].abs + @h[l][l].abs 486 if (s == 0.0) 487 s = norm 488 end 489 if (@h[l][l-1].abs < eps * s) 490 break 491 end 492 l-=1 493 end 494 495 # Check for convergence 496 # One root found 497 498 if (l == n) 499 @h[n][n] = @h[n][n] + exshift 500 @d[n] = @h[n][n] 501 @e[n] = 0.0 502 n-=1 503 iter = 0 504 505 # Two roots found 506 507 elsif (l == n-1) 508 w = @h[n][n-1] * @h[n-1][n] 509 p = (@h[n-1][n-1] - @h[n][n]) / 2.0 510 q = p * p + w 511 z = Math.sqrt(q.abs) 512 @h[n][n] = @h[n][n] + exshift 513 @h[n-1][n-1] = @h[n-1][n-1] + exshift 514 x = @h[n][n] 515 516 # Real pair 517 518 if (q >= 0) 519 if (p >= 0) 520 z = p + z 521 else 522 z = p - z 523 end 524 @d[n-1] = x + z 525 @d[n] = @d[n-1] 526 if (z != 0.0) 527 @d[n] = x - w / z 528 end 529 @e[n-1] = 0.0 530 @e[n] = 0.0 531 x = @h[n][n-1] 532 s = x.abs + z.abs 533 p = x / s 534 q = z / s 535 r = Math.sqrt(p * p+q * q) 536 p /= r 537 q /= r 538 539 # Row modification 540 541 (n-1).upto(nn-1) do |j| 542 z = @h[n-1][j] 543 @h[n-1][j] = q * z + p * @h[n][j] 544 @h[n][j] = q * @h[n][j] - p * z 545 end 546 547 # Column modification 548 549 0.upto(n) do |i| 550 z = @h[i][n-1] 551 @h[i][n-1] = q * z + p * @h[i][n] 552 @h[i][n] = q * @h[i][n] - p * z 553 end 554 555 # Accumulate transformations 556 557 low.upto(high) do |i| 558 z = @v[i][n-1] 559 @v[i][n-1] = q * z + p * @v[i][n] 560 @v[i][n] = q * @v[i][n] - p * z 561 end 562 563 # Complex pair 564 565 else 566 @d[n-1] = x + p 567 @d[n] = x + p 568 @e[n-1] = z 569 @e[n] = -z 570 end 571 n -= 2 572 iter = 0 573 574 # No convergence yet 575 576 else 577 578 # Form shift 579 580 x = @h[n][n] 581 y = 0.0 582 w = 0.0 583 if (l < n) 584 y = @h[n-1][n-1] 585 w = @h[n][n-1] * @h[n-1][n] 586 end 587 588 # Wilkinson's original ad hoc shift 589 590 if (iter == 10) 591 exshift += x 592 low.upto(n) do |i| 593 @h[i][i] -= x 594 end 595 s = @h[n][n-1].abs + @h[n-1][n-2].abs 596 x = y = 0.75 * s 597 w = -0.4375 * s * s 598 end 599 600 # MATLAB's new ad hoc shift 601 602 if (iter == 30) 603 s = (y - x) / 2.0 604 s *= s + w 605 if (s > 0) 606 s = Math.sqrt(s) 607 if (y < x) 608 s = -s 609 end 610 s = x - w / ((y - x) / 2.0 + s) 611 low.upto(n) do |i| 612 @h[i][i] -= s 613 end 614 exshift += s 615 x = y = w = 0.964 616 end 617 end 618 619 iter = iter + 1 # (Could check iteration count here.) 620 621 # Look for two consecutive small sub-diagonal elements 622 623 m = n-2 624 while (m >= l) do 625 z = @h[m][m] 626 r = x - z 627 s = y - z 628 p = (r * s - w) / @h[m+1][m] + @h[m][m+1] 629 q = @h[m+1][m+1] - z - r - s 630 r = @h[m+2][m+1] 631 s = p.abs + q.abs + r.abs 632 p /= s 633 q /= s 634 r /= s 635 if (m == l) 636 break 637 end 638 if (@h[m][m-1].abs * (q.abs + r.abs) < 639 eps * (p.abs * (@h[m-1][m-1].abs + z.abs + 640 @h[m+1][m+1].abs))) 641 break 642 end 643 m-=1 644 end 645 646 (m+2).upto(n) do |i| 647 @h[i][i-2] = 0.0 648 if (i > m+2) 649 @h[i][i-3] = 0.0 650 end 651 end 652 653 # Double QR step involving rows l:n and columns m:n 654 655 m.upto(n-1) do |k| 656 notlast = (k != n-1) 657 if (k != m) 658 p = @h[k][k-1] 659 q = @h[k+1][k-1] 660 r = (notlast ? @h[k+2][k-1] : 0.0) 661 x = p.abs + q.abs + r.abs 662 next if x == 0 663 p /= x 664 q /= x 665 r /= x 666 end 667 s = Math.sqrt(p * p + q * q + r * r) 668 if (p < 0) 669 s = -s 670 end 671 if (s != 0) 672 if (k != m) 673 @h[k][k-1] = -s * x 674 elsif (l != m) 675 @h[k][k-1] = -@h[k][k-1] 676 end 677 p += s 678 x = p / s 679 y = q / s 680 z = r / s 681 q /= p 682 r /= p 683 684 # Row modification 685 686 k.upto(nn-1) do |j| 687 p = @h[k][j] + q * @h[k+1][j] 688 if (notlast) 689 p += r * @h[k+2][j] 690 @h[k+2][j] = @h[k+2][j] - p * z 691 end 692 @h[k][j] = @h[k][j] - p * x 693 @h[k+1][j] = @h[k+1][j] - p * y 694 end 695 696 # Column modification 697 698 0.upto([n, k+3].min) do |i| 699 p = x * @h[i][k] + y * @h[i][k+1] 700 if (notlast) 701 p += z * @h[i][k+2] 702 @h[i][k+2] = @h[i][k+2] - p * r 703 end 704 @h[i][k] = @h[i][k] - p 705 @h[i][k+1] = @h[i][k+1] - p * q 706 end 707 708 # Accumulate transformations 709 710 low.upto(high) do |i| 711 p = x * @v[i][k] + y * @v[i][k+1] 712 if (notlast) 713 p += z * @v[i][k+2] 714 @v[i][k+2] = @v[i][k+2] - p * r 715 end 716 @v[i][k] = @v[i][k] - p 717 @v[i][k+1] = @v[i][k+1] - p * q 718 end 719 end # (s != 0) 720 end # k loop 721 end # check convergence 722 end # while (n >= low) 723 724 # Backsubstitute to find vectors of upper triangular form 725 726 if (norm == 0.0) 727 return 728 end 729 730 (nn-1).downto(0) do |n| 731 p = @d[n] 732 q = @e[n] 733 734 # Real vector 735 736 if (q == 0) 737 l = n 738 @h[n][n] = 1.0 739 (n-1).downto(0) do |i| 740 w = @h[i][i] - p 741 r = 0.0 742 l.upto(n) do |j| 743 r += @h[i][j] * @h[j][n] 744 end 745 if (@e[i] < 0.0) 746 z = w 747 s = r 748 else 749 l = i 750 if (@e[i] == 0.0) 751 if (w != 0.0) 752 @h[i][n] = -r / w 753 else 754 @h[i][n] = -r / (eps * norm) 755 end 756 757 # Solve real equations 758 759 else 760 x = @h[i][i+1] 761 y = @h[i+1][i] 762 q = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] 763 t = (x * s - z * r) / q 764 @h[i][n] = t 765 if (x.abs > z.abs) 766 @h[i+1][n] = (-r - w * t) / x 767 else 768 @h[i+1][n] = (-s - y * t) / z 769 end 770 end 771 772 # Overflow control 773 774 t = @h[i][n].abs 775 if ((eps * t) * t > 1) 776 i.upto(n) do |j| 777 @h[j][n] = @h[j][n] / t 778 end 779 end 780 end 781 end 782 783 # Complex vector 784 785 elsif (q < 0) 786 l = n-1 787 788 # Last vector component imaginary so matrix is triangular 789 790 if (@h[n][n-1].abs > @h[n-1][n].abs) 791 @h[n-1][n-1] = q / @h[n][n-1] 792 @h[n-1][n] = -(@h[n][n] - p) / @h[n][n-1] 793 else 794 cdivr, cdivi = cdiv(0.0, -@h[n-1][n], @h[n-1][n-1]-p, q) 795 @h[n-1][n-1] = cdivr 796 @h[n-1][n] = cdivi 797 end 798 @h[n][n-1] = 0.0 799 @h[n][n] = 1.0 800 (n-2).downto(0) do |i| 801 ra = 0.0 802 sa = 0.0 803 l.upto(n) do |j| 804 ra = ra + @h[i][j] * @h[j][n-1] 805 sa = sa + @h[i][j] * @h[j][n] 806 end 807 w = @h[i][i] - p 808 809 if (@e[i] < 0.0) 810 z = w 811 r = ra 812 s = sa 813 else 814 l = i 815 if (@e[i] == 0) 816 cdivr, cdivi = cdiv(-ra, -sa, w, q) 817 @h[i][n-1] = cdivr 818 @h[i][n] = cdivi 819 else 820 821 # Solve complex equations 822 823 x = @h[i][i+1] 824 y = @h[i+1][i] 825 vr = (@d[i] - p) * (@d[i] - p) + @e[i] * @e[i] - q * q 826 vi = (@d[i] - p) * 2.0 * q 827 if (vr == 0.0 && vi == 0.0) 828 vr = eps * norm * (w.abs + q.abs + 829 x.abs + y.abs + z.abs) 830 end 831 cdivr, cdivi = cdiv(x*r-z*ra+q*sa, x*s-z*sa-q*ra, vr, vi) 832 @h[i][n-1] = cdivr 833 @h[i][n] = cdivi 834 if (x.abs > (z.abs + q.abs)) 835 @h[i+1][n-1] = (-ra - w * @h[i][n-1] + q * @h[i][n]) / x 836 @h[i+1][n] = (-sa - w * @h[i][n] - q * @h[i][n-1]) / x 837 else 838 cdivr, cdivi = cdiv(-r-y*@h[i][n-1], -s-y*@h[i][n], z, q) 839 @h[i+1][n-1] = cdivr 840 @h[i+1][n] = cdivi 841 end 842 end 843 844 # Overflow control 845 846 t = [@h[i][n-1].abs, @h[i][n].abs].max 847 if ((eps * t) * t > 1) 848 i.upto(n) do |j| 849 @h[j][n-1] = @h[j][n-1] / t 850 @h[j][n] = @h[j][n] / t 851 end 852 end 853 end 854 end 855 end 856 end 857 858 # Vectors of isolated roots 859 860 nn.times do |i| 861 if (i < low || i > high) 862 i.upto(nn-1) do |j| 863 @v[i][j] = @h[i][j] 864 end 865 end 866 end 867 868 # Back transformation to get eigenvectors of original matrix 869 870 (nn-1).downto(low) do |j| 871 low.upto(high) do |i| 872 z = 0.0 873 low.upto([j, high].min) do |k| 874 z += @v[i][k] * @h[k][j] 875 end 876 @v[i][j] = z 877 end 878 end 879 end 880 881 end 882end 883