Conditions | 71 |
Paths | > 20000 |
Total Lines | 414 |
Code Lines | 270 |
Lines | 48 |
Ratio | 11.59 % |
Changes | 0 |
Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.
For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.
Commonly applied refactorings include:
If many parameters/temporary variables are present:
1 | <?php |
||
550 | private function hqr2(): void |
||
551 | { |
||
552 | // Initialize |
||
553 | $nn = $this->n; |
||
554 | $n = $nn - 1; |
||
555 | $low = 0; |
||
556 | $high = $nn - 1; |
||
557 | $eps = pow(2.0, -52.0); |
||
558 | $exshift = 0.0; |
||
559 | $p = $q = $r = $s = $z = 0; |
||
560 | // Store roots isolated by balanc and compute matrix norm |
||
561 | $norm = 0.0; |
||
562 | |||
563 | for ($i = 0; $i < $nn; ++$i) { |
||
564 | if (($i < $low) or ($i > $high)) { |
||
565 | $this->d[$i] = $this->H[$i][$i]; |
||
566 | $this->e[$i] = 0.0; |
||
567 | } |
||
568 | |||
569 | for ($j = max($i - 1, 0); $j < $nn; ++$j) { |
||
570 | $norm = $norm + abs($this->H[$i][$j]); |
||
571 | } |
||
572 | } |
||
573 | |||
574 | // Outer loop over eigenvalue index |
||
575 | $iter = 0; |
||
576 | while ($n >= $low) { |
||
577 | // Look for single small sub-diagonal element |
||
578 | $l = $n; |
||
579 | while ($l > $low) { |
||
580 | $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]); |
||
581 | if ($s == 0.0) { |
||
582 | $s = $norm; |
||
583 | } |
||
584 | |||
585 | if (abs($this->H[$l][$l - 1]) < $eps * $s) { |
||
586 | break; |
||
587 | } |
||
588 | |||
589 | --$l; |
||
590 | } |
||
591 | |||
592 | // Check for convergence |
||
593 | // One root found |
||
594 | if ($l == $n) { |
||
595 | $this->H[$n][$n] = $this->H[$n][$n] + $exshift; |
||
596 | $this->d[$n] = $this->H[$n][$n]; |
||
597 | $this->e[$n] = 0.0; |
||
598 | --$n; |
||
599 | $iter = 0; |
||
600 | // Two roots found |
||
601 | } elseif ($l == $n - 1) { |
||
602 | $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; |
||
603 | $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0; |
||
604 | $q = $p * $p + $w; |
||
605 | $z = sqrt(abs($q)); |
||
606 | $this->H[$n][$n] = $this->H[$n][$n] + $exshift; |
||
607 | $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift; |
||
608 | $x = $this->H[$n][$n]; |
||
609 | // Real pair |
||
610 | if ($q >= 0) { |
||
611 | if ($p >= 0) { |
||
612 | $z = $p + $z; |
||
613 | } else { |
||
614 | $z = $p - $z; |
||
615 | } |
||
616 | |||
617 | $this->d[$n - 1] = $x + $z; |
||
618 | $this->d[$n] = $this->d[$n - 1]; |
||
619 | if ($z != 0.0) { |
||
620 | $this->d[$n] = $x - $w / $z; |
||
621 | } |
||
622 | |||
623 | $this->e[$n - 1] = 0.0; |
||
624 | $this->e[$n] = 0.0; |
||
625 | $x = $this->H[$n][$n - 1]; |
||
626 | $s = abs($x) + abs($z); |
||
627 | $p = $x / $s; |
||
628 | $q = $z / $s; |
||
629 | $r = sqrt($p * $p + $q * $q); |
||
630 | $p = $p / $r; |
||
631 | $q = $q / $r; |
||
632 | // Row modification |
||
633 | View Code Duplication | for ($j = $n - 1; $j < $nn; ++$j) { |
|
634 | $z = $this->H[$n - 1][$j]; |
||
635 | $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j]; |
||
636 | $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; |
||
637 | } |
||
638 | |||
639 | // Column modification |
||
640 | View Code Duplication | for ($i = 0; $i <= $n; ++$i) { |
|
641 | $z = $this->H[$i][$n - 1]; |
||
642 | $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n]; |
||
643 | $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; |
||
644 | } |
||
645 | |||
646 | // Accumulate transformations |
||
647 | View Code Duplication | for ($i = $low; $i <= $high; ++$i) { |
|
648 | $z = $this->V[$i][$n - 1]; |
||
649 | $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n]; |
||
650 | $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; |
||
651 | } |
||
652 | |||
653 | // Complex pair |
||
654 | } else { |
||
655 | $this->d[$n - 1] = $x + $p; |
||
656 | $this->d[$n] = $x + $p; |
||
657 | $this->e[$n - 1] = $z; |
||
658 | $this->e[$n] = -$z; |
||
659 | } |
||
660 | |||
661 | $n = $n - 2; |
||
662 | $iter = 0; |
||
663 | // No convergence yet |
||
664 | } else { |
||
665 | // Form shift |
||
666 | $x = $this->H[$n][$n]; |
||
667 | $y = 0.0; |
||
668 | $w = 0.0; |
||
669 | if ($l < $n) { |
||
670 | $y = $this->H[$n - 1][$n - 1]; |
||
671 | $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n]; |
||
672 | } |
||
673 | |||
674 | // Wilkinson's original ad hoc shift |
||
675 | if ($iter == 10) { |
||
676 | $exshift += $x; |
||
677 | for ($i = $low; $i <= $n; ++$i) { |
||
678 | $this->H[$i][$i] -= $x; |
||
679 | } |
||
680 | |||
681 | $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]); |
||
682 | $x = $y = 0.75 * $s; |
||
683 | $w = -0.4375 * $s * $s; |
||
684 | } |
||
685 | |||
686 | // MATLAB's new ad hoc shift |
||
687 | if ($iter == 30) { |
||
688 | $s = ($y - $x) / 2.0; |
||
689 | $s = $s * $s + $w; |
||
690 | if ($s > 0) { |
||
691 | $s = sqrt($s); |
||
692 | if ($y < $x) { |
||
693 | $s = -$s; |
||
694 | } |
||
695 | |||
696 | $s = $x - $w / (($y - $x) / 2.0 + $s); |
||
697 | for ($i = $low; $i <= $n; ++$i) { |
||
698 | $this->H[$i][$i] -= $s; |
||
699 | } |
||
700 | |||
701 | $exshift += $s; |
||
702 | $x = $y = $w = 0.964; |
||
703 | } |
||
704 | } |
||
705 | |||
706 | // Could check iteration count here. |
||
707 | $iter = $iter + 1; |
||
708 | // Look for two consecutive small sub-diagonal elements |
||
709 | $m = $n - 2; |
||
710 | while ($m >= $l) { |
||
711 | $z = $this->H[$m][$m]; |
||
712 | $r = $x - $z; |
||
713 | $s = $y - $z; |
||
714 | $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1]; |
||
715 | $q = $this->H[$m + 1][$m + 1] - $z - $r - $s; |
||
716 | $r = $this->H[$m + 2][$m + 1]; |
||
717 | $s = abs($p) + abs($q) + abs($r); |
||
718 | $p = $p / $s; |
||
719 | $q = $q / $s; |
||
720 | $r = $r / $s; |
||
721 | if ($m == $l) { |
||
722 | break; |
||
723 | } |
||
724 | |||
725 | if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < |
||
726 | $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) { |
||
727 | break; |
||
728 | } |
||
729 | |||
730 | --$m; |
||
731 | } |
||
732 | |||
733 | for ($i = $m + 2; $i <= $n; ++$i) { |
||
734 | $this->H[$i][$i - 2] = 0.0; |
||
735 | if ($i > $m + 2) { |
||
736 | $this->H[$i][$i - 3] = 0.0; |
||
737 | } |
||
738 | } |
||
739 | |||
740 | // Double QR step involving rows l:n and columns m:n |
||
741 | for ($k = $m; $k <= $n - 1; ++$k) { |
||
742 | $notlast = ($k != $n - 1); |
||
743 | if ($k != $m) { |
||
744 | $p = $this->H[$k][$k - 1]; |
||
745 | $q = $this->H[$k + 1][$k - 1]; |
||
746 | $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); |
||
747 | $x = abs($p) + abs($q) + abs($r); |
||
748 | if ($x != 0.0) { |
||
749 | $p = $p / $x; |
||
750 | $q = $q / $x; |
||
751 | $r = $r / $x; |
||
752 | } |
||
753 | } |
||
754 | |||
755 | if ($x == 0.0) { |
||
756 | break; |
||
757 | } |
||
758 | |||
759 | $s = sqrt($p * $p + $q * $q + $r * $r); |
||
760 | if ($p < 0) { |
||
761 | $s = -$s; |
||
762 | } |
||
763 | |||
764 | if ($s != 0) { |
||
765 | if ($k != $m) { |
||
766 | $this->H[$k][$k - 1] = -$s * $x; |
||
767 | } elseif ($l != $m) { |
||
768 | $this->H[$k][$k - 1] = -$this->H[$k][$k - 1]; |
||
769 | } |
||
770 | |||
771 | $p = $p + $s; |
||
772 | $x = $p / $s; |
||
773 | $y = $q / $s; |
||
774 | $z = $r / $s; |
||
775 | $q = $q / $p; |
||
776 | $r = $r / $p; |
||
777 | // Row modification |
||
778 | View Code Duplication | for ($j = $k; $j < $nn; ++$j) { |
|
779 | $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j]; |
||
780 | if ($notlast) { |
||
781 | $p = $p + $r * $this->H[$k + 2][$j]; |
||
782 | $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z; |
||
783 | } |
||
784 | |||
785 | $this->H[$k][$j] = $this->H[$k][$j] - $p * $x; |
||
786 | $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y; |
||
787 | } |
||
788 | |||
789 | // Column modification |
||
790 | View Code Duplication | for ($i = 0; $i <= min($n, $k + 3); ++$i) { |
|
791 | $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1]; |
||
792 | if ($notlast) { |
||
793 | $p = $p + $z * $this->H[$i][$k + 2]; |
||
794 | $this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r; |
||
795 | } |
||
796 | |||
797 | $this->H[$i][$k] = $this->H[$i][$k] - $p; |
||
798 | $this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q; |
||
799 | } |
||
800 | |||
801 | // Accumulate transformations |
||
802 | View Code Duplication | for ($i = $low; $i <= $high; ++$i) { |
|
803 | $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1]; |
||
804 | if ($notlast) { |
||
805 | $p = $p + $z * $this->V[$i][$k + 2]; |
||
806 | $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r; |
||
807 | } |
||
808 | |||
809 | $this->V[$i][$k] = $this->V[$i][$k] - $p; |
||
810 | $this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q; |
||
811 | } |
||
812 | } // ($s != 0) |
||
813 | } // k loop |
||
814 | } // check convergence |
||
815 | } // while ($n >= $low) |
||
816 | |||
817 | // Backsubstitute to find vectors of upper triangular form |
||
818 | if ($norm == 0.0) { |
||
819 | return; |
||
820 | } |
||
821 | |||
822 | for ($n = $nn - 1; $n >= 0; --$n) { |
||
823 | $p = $this->d[$n]; |
||
824 | $q = $this->e[$n]; |
||
825 | // Real vector |
||
826 | if ($q == 0) { |
||
827 | $l = $n; |
||
828 | $this->H[$n][$n] = 1.0; |
||
829 | for ($i = $n - 1; $i >= 0; --$i) { |
||
830 | $w = $this->H[$i][$i] - $p; |
||
831 | $r = 0.0; |
||
832 | View Code Duplication | for ($j = $l; $j <= $n; ++$j) { |
|
833 | $r = $r + $this->H[$i][$j] * $this->H[$j][$n]; |
||
834 | } |
||
835 | |||
836 | if ($this->e[$i] < 0.0) { |
||
837 | $z = $w; |
||
838 | $s = $r; |
||
839 | } else { |
||
840 | $l = $i; |
||
841 | if ($this->e[$i] == 0.0) { |
||
842 | if ($w != 0.0) { |
||
843 | $this->H[$i][$n] = -$r / $w; |
||
844 | } else { |
||
845 | $this->H[$i][$n] = -$r / ($eps * $norm); |
||
846 | } |
||
847 | |||
848 | // Solve real equations |
||
849 | } else { |
||
850 | $x = $this->H[$i][$i + 1]; |
||
851 | $y = $this->H[$i + 1][$i]; |
||
852 | $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i]; |
||
853 | $t = ($x * $s - $z * $r) / $q; |
||
854 | $this->H[$i][$n] = $t; |
||
855 | if (abs($x) > abs($z)) { |
||
856 | $this->H[$i + 1][$n] = (-$r - $w * $t) / $x; |
||
857 | } else { |
||
858 | $this->H[$i + 1][$n] = (-$s - $y * $t) / $z; |
||
859 | } |
||
860 | } |
||
861 | |||
862 | // Overflow control |
||
863 | $t = abs($this->H[$i][$n]); |
||
864 | if (($eps * $t) * $t > 1) { |
||
865 | View Code Duplication | for ($j = $i; $j <= $n; ++$j) { |
|
866 | $this->H[$j][$n] = $this->H[$j][$n] / $t; |
||
867 | } |
||
868 | } |
||
869 | } |
||
870 | } |
||
871 | |||
872 | // Complex vector |
||
873 | } elseif ($q < 0) { |
||
874 | $l = $n - 1; |
||
875 | // Last vector component imaginary so matrix is triangular |
||
876 | if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) { |
||
877 | $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1]; |
||
878 | $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1]; |
||
879 | } else { |
||
880 | $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q); |
||
881 | $this->H[$n - 1][$n - 1] = $this->cdivr; |
||
882 | $this->H[$n - 1][$n] = $this->cdivi; |
||
883 | } |
||
884 | |||
885 | $this->H[$n][$n - 1] = 0.0; |
||
886 | $this->H[$n][$n] = 1.0; |
||
887 | for ($i = $n - 2; $i >= 0; --$i) { |
||
888 | // double ra,sa,vr,vi; |
||
889 | $ra = 0.0; |
||
890 | $sa = 0.0; |
||
891 | for ($j = $l; $j <= $n; ++$j) { |
||
892 | $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n - 1]; |
||
893 | $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n]; |
||
894 | } |
||
895 | |||
896 | $w = $this->H[$i][$i] - $p; |
||
897 | if ($this->e[$i] < 0.0) { |
||
898 | $z = $w; |
||
899 | $r = $ra; |
||
900 | $s = $sa; |
||
901 | } else { |
||
902 | $l = $i; |
||
903 | if ($this->e[$i] == 0) { |
||
904 | $this->cdiv(-$ra, -$sa, $w, $q); |
||
905 | $this->H[$i][$n - 1] = $this->cdivr; |
||
906 | $this->H[$i][$n] = $this->cdivi; |
||
907 | } else { |
||
908 | // Solve complex equations |
||
909 | $x = $this->H[$i][$i + 1]; |
||
910 | $y = $this->H[$i + 1][$i]; |
||
911 | $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q; |
||
912 | $vi = ($this->d[$i] - $p) * 2.0 * $q; |
||
913 | if ($vr == 0.0 & $vi == 0.0) { |
||
914 | $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z)); |
||
915 | } |
||
916 | |||
917 | $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi); |
||
918 | $this->H[$i][$n - 1] = $this->cdivr; |
||
919 | $this->H[$i][$n] = $this->cdivi; |
||
920 | if (abs($x) > (abs($z) + abs($q))) { |
||
921 | $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x; |
||
922 | $this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x; |
||
923 | } else { |
||
924 | $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q); |
||
925 | $this->H[$i + 1][$n - 1] = $this->cdivr; |
||
926 | $this->H[$i + 1][$n] = $this->cdivi; |
||
927 | } |
||
928 | } |
||
929 | |||
930 | // Overflow control |
||
931 | $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n])); |
||
932 | if (($eps * $t) * $t > 1) { |
||
933 | for ($j = $i; $j <= $n; ++$j) { |
||
934 | $this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t; |
||
935 | $this->H[$j][$n] = $this->H[$j][$n] / $t; |
||
936 | } |
||
937 | } |
||
938 | } // end else |
||
939 | } // end for |
||
940 | } // end else for complex case |
||
941 | } // end for |
||
942 | |||
943 | // Vectors of isolated roots |
||
944 | for ($i = 0; $i < $nn; ++$i) { |
||
945 | if ($i < $low | $i > $high) { |
||
946 | for ($j = $i; $j < $nn; ++$j) { |
||
947 | $this->V[$i][$j] = $this->H[$i][$j]; |
||
948 | } |
||
949 | } |
||
950 | } |
||
951 | |||
952 | // Back transformation to get eigenvectors of original matrix |
||
953 | for ($j = $nn - 1; $j >= $low; --$j) { |
||
954 | for ($i = $low; $i <= $high; ++$i) { |
||
955 | $z = 0.0; |
||
956 | for ($k = $low; $k <= min($j, $high); ++$k) { |
||
957 | $z = $z + $this->V[$i][$k] * $this->H[$k][$j]; |
||
958 | } |
||
959 | |||
960 | $this->V[$i][$j] = $z; |
||
961 | } |
||
962 | } |
||
963 | } |
||
964 | } |
||
965 |
If the size of the collection does not change during the iteration, it is generally a good practice to compute it beforehand, and not on each iteration: