= -1; --$k) { if ($k == -1) { break; } if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) { $e[$k] = 0.0; break; } } if ($k == $p - 2) { $kase = 4; } else { for ($ks = $p - 1; $ks >= $k; --$ks) { if ($ks == $k) { break; } $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.); if (abs($this->s[$ks]) <= $eps * $t) { $this->s[$ks] = 0.0; break; } } if ($ks == $k) { $kase = 3; } elseif ($ks == $p-1) { $kase = 1; } else { $kase = 2; $k = $ks; } } ++$k; // Perform the task indicated by kase. switch ($kase) { // Deflate negligible s(p). case 1: $f = $e[$p-2]; $e[$p-2] = 0.0; for ($j = $p - 2; $j >= $k; --$j) { $t = hypo($this->s[$j], $f); $cs = $this->s[$j] / $t; $sn = $f / $t; $this->s[$j] = $t; if ($j != $k) { $f = -$sn * $e[$j-1]; $e[$j-1] = $cs * $e[$j-1]; } if ($wantv) { for ($i = 0; $i < $this->n; ++$i) { $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1]; $this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1]; $this->V[$i][$j] = $t; } } } break; // Split at negligible s(k). case 2: $f = $e[$k-1]; $e[$k-1] = 0.0; for ($j = $k; $j < $p; ++$j) { $t = hypo($this->s[$j], $f); $cs = $this->s[$j] / $t; $sn = $f / $t; $this->s[$j] = $t; $f = -$sn * $e[$j]; $e[$j] = $cs * $e[$j]; if ($wantu) { for ($i = 0; $i < $this->m; ++$i) { $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1]; $this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1]; $this->U[$i][$j] = $t; } } } break; // Perform one qr step. case 3: // Calculate the shift. $scale = max(max(max(max(abs($this->s[$p-1]), abs($this->s[$p-2])), abs($e[$p-2])), abs($this->s[$k])), abs($e[$k])); $sp = $this->s[$p-1] / $scale; $spm1 = $this->s[$p-2] / $scale; $epm1 = $e[$p-2] / $scale; $sk = $this->s[$k] / $scale; $ek = $e[$k] / $scale; $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0; $c = ($sp * $epm1) * ($sp * $epm1); $shift = 0.0; if (($b != 0.0) || ($c != 0.0)) { $shift = sqrt($b * $b + $c); if ($b < 0.0) { $shift = -$shift; } $shift = $c / ($b + $shift); } $f = ($sk + $sp) * ($sk - $sp) + $shift; $g = $sk * $ek; // Chase zeros. for ($j = $k; $j < $p-1; ++$j) { $t = hypo($f, $g); $cs = $f/$t; $sn = $g/$t; if ($j != $k) { $e[$j-1] = $t; } $f = $cs * $this->s[$j] + $sn * $e[$j]; $e[$j] = $cs * $e[$j] - $sn * $this->s[$j]; $g = $sn * $this->s[$j+1]; $this->s[$j+1] = $cs * $this->s[$j+1]; if ($wantv) { for ($i = 0; $i < $this->n; ++$i) { $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1]; $this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1]; $this->V[$i][$j] = $t; } } $t = hypo($f, $g); $cs = $f/$t; $sn = $g/$t; $this->s[$j] = $t; $f = $cs * $e[$j] + $sn * $this->s[$j+1]; $this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1]; $g = $sn * $e[$j+1]; $e[$j+1] = $cs * $e[$j+1]; if ($wantu && ($j < $this->m - 1)) { for ($i = 0; $i < $this->m; ++$i) { $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1]; $this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1]; $this->U[$i][$j] = $t; } } } $e[$p-2] = $f; $iter = $iter + 1; break; // Convergence. case 4: // Make the singular values positive. if ($this->s[$k] <= 0.0) { $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0); if ($wantv) { for ($i = 0; $i <= $pp; ++$i) { $this->V[$i][$k] = -$this->V[$i][$k]; } } } // Order the singular values. while ($k < $pp) { if ($this->s[$k] >= $this->s[$k+1]) { break; } $t = $this->s[$k]; $this->s[$k] = $this->s[$k+1]; $this->s[$k+1] = $t; if ($wantv and ($k < $this->n - 1)) { for ($i = 0; $i < $this->n; ++$i) { $t = $this->V[$i][$k+1]; $this->V[$i][$k+1] = $this->V[$i][$k]; $this->V[$i][$k] = $t; } } if ($wantu and ($k < $this->m-1)) { for ($i = 0; $i < $this->m; ++$i) { $t = $this->U[$i][$k+1]; $this->U[$i][$k+1] = $this->U[$i][$k]; $this->U[$i][$k] = $t; } } ++$k; } $iter = 0; --$p; break; } // end switch } // end while } // end constructor /** * Return the left singular vectors * * @access public * @return U */ public function getU() { return new Matrix($this->U, $this->m, min($this->m + 1, $this->n)); } /** * Return the right singular vectors * * @access public * @return V */ public function getV() { return new Matrix($this->V); } /** * Return the one-dimensional array of singular values * * @access public * @return diagonal of S. */ public function getSingularValues() { return $this->s; } /** * Return the diagonal matrix of singular values * * @access public * @return S */ public function getS() { for ($i = 0; $i < $this->n; ++$i) { for ($j = 0; $j < $this->n; ++$j) { $S[$i][$j] = 0.0; } $S[$i][$i] = $this->s[$i]; } return new Matrix($S); } /** * Two norm * * @access public * @return max(S) */ public function norm2() { return $this->s[0]; } /** * Two norm condition number * * @access public * @return max(S)/min(S) */ public function cond() { return $this->s[0] / $this->s[min($this->m, $this->n) - 1]; } /** * Effective numerical matrix rank * * @access public * @return Number of nonnegligible singular values. */ public function rank() { $eps = pow(2.0, -52.0); $tol = max($this->m, $this->n) * $this->s[0] * $eps; $r = 0; for ($i = 0; $i < count($this->s); ++$i) { if ($this->s[$i] > $tol) { ++$r; } } return $r; } }