| Conditions | 95 |
| Paths | > 20000 |
| Total Lines | 373 |
| Code Lines | 241 |
| Lines | 76 |
| Ratio | 20.38 % |
| 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 |
||
| 61 | public function __construct($Arg) { |
||
| 62 | |||
| 63 | // Initialize. |
||
| 64 | $A = $Arg->getArrayCopy(); |
||
| 65 | $this->m = $Arg->getRowDimension(); |
||
| 66 | $this->n = $Arg->getColumnDimension(); |
||
| 67 | $nu = min($this->m, $this->n); |
||
| 68 | $e = array(); |
||
| 69 | $work = array(); |
||
| 70 | $wantu = true; |
||
| 71 | $wantv = true; |
||
| 72 | $nct = min($this->m - 1, $this->n); |
||
| 73 | $nrt = max(0, min($this->n - 2, $this->m)); |
||
| 74 | |||
| 75 | // Reduce A to bidiagonal form, storing the diagonal elements |
||
| 76 | // in s and the super-diagonal elements in e. |
||
| 77 | for ($k = 0; $k < max($nct,$nrt); ++$k) { |
||
| 78 | |||
| 79 | if ($k < $nct) { |
||
| 80 | // Compute the transformation for the k-th column and |
||
| 81 | // place the k-th diagonal in s[$k]. |
||
| 82 | // Compute 2-norm of k-th column without under/overflow. |
||
| 83 | $this->s[$k] = 0; |
||
| 84 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 85 | $this->s[$k] = hypo($this->s[$k], $A[$i][$k]); |
||
| 86 | } |
||
| 87 | if ($this->s[$k] != 0.0) { |
||
| 88 | if ($A[$k][$k] < 0.0) { |
||
| 89 | $this->s[$k] = -$this->s[$k]; |
||
| 90 | } |
||
| 91 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 92 | $A[$i][$k] /= $this->s[$k]; |
||
| 93 | } |
||
| 94 | $A[$k][$k] += 1.0; |
||
| 95 | } |
||
| 96 | $this->s[$k] = -$this->s[$k]; |
||
| 97 | } |
||
| 98 | |||
| 99 | for ($j = $k + 1; $j < $this->n; ++$j) { |
||
| 100 | if (($k < $nct) & ($this->s[$k] != 0.0)) { |
||
| 101 | // Apply the transformation. |
||
| 102 | $t = 0; |
||
| 103 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 104 | $t += $A[$i][$k] * $A[$i][$j]; |
||
| 105 | } |
||
| 106 | $t = -$t / $A[$k][$k]; |
||
| 107 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 108 | $A[$i][$j] += $t * $A[$i][$k]; |
||
| 109 | } |
||
| 110 | // Place the k-th row of A into e for the |
||
| 111 | // subsequent calculation of the row transformation. |
||
| 112 | $e[$j] = $A[$k][$j]; |
||
| 113 | } |
||
| 114 | } |
||
| 115 | |||
| 116 | if ($wantu AND ($k < $nct)) { |
||
| 117 | // Place the transformation in U for subsequent back |
||
| 118 | // multiplication. |
||
| 119 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 120 | $this->U[$i][$k] = $A[$i][$k]; |
||
| 121 | } |
||
| 122 | } |
||
| 123 | |||
| 124 | if ($k < $nrt) { |
||
| 125 | // Compute the k-th row transformation and place the |
||
| 126 | // k-th super-diagonal in e[$k]. |
||
| 127 | // Compute 2-norm without under/overflow. |
||
| 128 | $e[$k] = 0; |
||
| 129 | for ($i = $k + 1; $i < $this->n; ++$i) { |
||
| 130 | $e[$k] = hypo($e[$k], $e[$i]); |
||
| 131 | } |
||
| 132 | if ($e[$k] != 0.0) { |
||
| 133 | if ($e[$k+1] < 0.0) { |
||
| 134 | $e[$k] = -$e[$k]; |
||
| 135 | } |
||
| 136 | for ($i = $k + 1; $i < $this->n; ++$i) { |
||
| 137 | $e[$i] /= $e[$k]; |
||
| 138 | } |
||
| 139 | $e[$k+1] += 1.0; |
||
| 140 | } |
||
| 141 | $e[$k] = -$e[$k]; |
||
| 142 | if (($k+1 < $this->m) AND ($e[$k] != 0.0)) { |
||
| 143 | // Apply the transformation. |
||
| 144 | View Code Duplication | for ($i = $k+1; $i < $this->m; ++$i) { |
|
| 145 | $work[$i] = 0.0; |
||
| 146 | } |
||
| 147 | for ($j = $k+1; $j < $this->n; ++$j) { |
||
| 148 | for ($i = $k+1; $i < $this->m; ++$i) { |
||
| 149 | $work[$i] += $e[$j] * $A[$i][$j]; |
||
| 150 | } |
||
| 151 | } |
||
| 152 | for ($j = $k + 1; $j < $this->n; ++$j) { |
||
| 153 | $t = -$e[$j] / $e[$k+1]; |
||
| 154 | for ($i = $k + 1; $i < $this->m; ++$i) { |
||
| 155 | $A[$i][$j] += $t * $work[$i]; |
||
| 156 | } |
||
| 157 | } |
||
| 158 | } |
||
| 159 | View Code Duplication | if ($wantv) { |
|
| 160 | // Place the transformation in V for subsequent |
||
| 161 | // back multiplication. |
||
| 162 | for ($i = $k + 1; $i < $this->n; ++$i) { |
||
| 163 | $this->V[$i][$k] = $e[$i]; |
||
| 164 | } |
||
| 165 | } |
||
| 166 | } |
||
| 167 | } |
||
| 168 | |||
| 169 | // Set up the final bidiagonal matrix or order p. |
||
| 170 | $p = min($this->n, $this->m + 1); |
||
| 171 | if ($nct < $this->n) { |
||
| 172 | $this->s[$nct] = $A[$nct][$nct]; |
||
| 173 | } |
||
| 174 | if ($this->m < $p) { |
||
| 175 | $this->s[$p-1] = 0.0; |
||
| 176 | } |
||
| 177 | if ($nrt + 1 < $p) { |
||
| 178 | $e[$nrt] = $A[$nrt][$p-1]; |
||
| 179 | } |
||
| 180 | $e[$p-1] = 0.0; |
||
| 181 | // If required, generate U. |
||
| 182 | if ($wantu) { |
||
| 183 | for ($j = $nct; $j < $nu; ++$j) { |
||
| 184 | for ($i = 0; $i < $this->m; ++$i) { |
||
| 185 | $this->U[$i][$j] = 0.0; |
||
| 186 | } |
||
| 187 | $this->U[$j][$j] = 1.0; |
||
| 188 | } |
||
| 189 | for ($k = $nct - 1; $k >= 0; --$k) { |
||
| 190 | if ($this->s[$k] != 0.0) { |
||
| 191 | View Code Duplication | for ($j = $k + 1; $j < $nu; ++$j) { |
|
| 192 | $t = 0; |
||
| 193 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 194 | $t += $this->U[$i][$k] * $this->U[$i][$j]; |
||
| 195 | } |
||
| 196 | $t = -$t / $this->U[$k][$k]; |
||
| 197 | for ($i = $k; $i < $this->m; ++$i) { |
||
| 198 | $this->U[$i][$j] += $t * $this->U[$i][$k]; |
||
| 199 | } |
||
| 200 | } |
||
| 201 | for ($i = $k; $i < $this->m; ++$i ) { |
||
| 202 | $this->U[$i][$k] = -$this->U[$i][$k]; |
||
| 203 | } |
||
| 204 | $this->U[$k][$k] = 1.0 + $this->U[$k][$k]; |
||
| 205 | View Code Duplication | for ($i = 0; $i < $k - 1; ++$i) { |
|
| 206 | $this->U[$i][$k] = 0.0; |
||
| 207 | } |
||
| 208 | } else { |
||
| 209 | View Code Duplication | for ($i = 0; $i < $this->m; ++$i) { |
|
| 210 | $this->U[$i][$k] = 0.0; |
||
| 211 | } |
||
| 212 | $this->U[$k][$k] = 1.0; |
||
| 213 | } |
||
| 214 | } |
||
| 215 | } |
||
| 216 | |||
| 217 | // If required, generate V. |
||
| 218 | if ($wantv) { |
||
| 219 | for ($k = $this->n - 1; $k >= 0; --$k) { |
||
| 220 | if (($k < $nrt) AND ($e[$k] != 0.0)) { |
||
| 221 | for ($j = $k + 1; $j < $nu; ++$j) { |
||
| 222 | $t = 0; |
||
| 223 | for ($i = $k + 1; $i < $this->n; ++$i) { |
||
| 224 | $t += $this->V[$i][$k]* $this->V[$i][$j]; |
||
| 225 | } |
||
| 226 | $t = -$t / $this->V[$k+1][$k]; |
||
| 227 | for ($i = $k + 1; $i < $this->n; ++$i) { |
||
| 228 | $this->V[$i][$j] += $t * $this->V[$i][$k]; |
||
| 229 | } |
||
| 230 | } |
||
| 231 | } |
||
| 232 | View Code Duplication | for ($i = 0; $i < $this->n; ++$i) { |
|
| 233 | $this->V[$i][$k] = 0.0; |
||
| 234 | } |
||
| 235 | $this->V[$k][$k] = 1.0; |
||
| 236 | } |
||
| 237 | } |
||
| 238 | |||
| 239 | // Main iteration loop for the singular values. |
||
| 240 | $pp = $p - 1; |
||
| 241 | $iter = 0; |
||
| 242 | $eps = pow(2.0, -52.0); |
||
| 243 | |||
| 244 | while ($p > 0) { |
||
| 245 | // Here is where a test for too many iterations would go. |
||
| 246 | // This section of the program inspects for negligible |
||
| 247 | // elements in the s and e arrays. On completion the |
||
| 248 | // variables kase and k are set as follows: |
||
| 249 | // kase = 1 if s(p) and e[k-1] are negligible and k<p |
||
| 250 | // kase = 2 if s(k) is negligible and k<p |
||
| 251 | // kase = 3 if e[k-1] is negligible, k<p, and |
||
| 252 | // s(k), ..., s(p) are not negligible (qr step). |
||
| 253 | // kase = 4 if e(p-1) is negligible (convergence). |
||
| 254 | for ($k = $p - 2; $k >= -1; --$k) { |
||
| 255 | if ($k == -1) { |
||
| 256 | break; |
||
| 257 | } |
||
| 258 | if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) { |
||
| 259 | $e[$k] = 0.0; |
||
| 260 | break; |
||
| 261 | } |
||
| 262 | } |
||
| 263 | if ($k == $p - 2) { |
||
| 264 | $kase = 4; |
||
| 265 | } else { |
||
| 266 | for ($ks = $p - 1; $ks >= $k; --$ks) { |
||
| 267 | if ($ks == $k) { |
||
| 268 | break; |
||
| 269 | } |
||
| 270 | $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.); |
||
| 271 | if (abs($this->s[$ks]) <= $eps * $t) { |
||
| 272 | $this->s[$ks] = 0.0; |
||
| 273 | break; |
||
| 274 | } |
||
| 275 | } |
||
| 276 | if ($ks == $k) { |
||
| 277 | $kase = 3; |
||
| 278 | } else if ($ks == $p-1) { |
||
| 279 | $kase = 1; |
||
| 280 | } else { |
||
| 281 | $kase = 2; |
||
| 282 | $k = $ks; |
||
| 283 | } |
||
| 284 | } |
||
| 285 | ++$k; |
||
| 286 | |||
| 287 | // Perform the task indicated by kase. |
||
| 288 | switch ($kase) { |
||
| 289 | // Deflate negligible s(p). |
||
| 290 | case 1: |
||
| 291 | $f = $e[$p-2]; |
||
| 292 | $e[$p-2] = 0.0; |
||
| 293 | for ($j = $p - 2; $j >= $k; --$j) { |
||
| 294 | $t = hypo($this->s[$j],$f); |
||
| 295 | $cs = $this->s[$j] / $t; |
||
| 296 | $sn = $f / $t; |
||
| 297 | $this->s[$j] = $t; |
||
| 298 | if ($j != $k) { |
||
| 299 | $f = -$sn * $e[$j-1]; |
||
| 300 | $e[$j-1] = $cs * $e[$j-1]; |
||
| 301 | } |
||
| 302 | View Code Duplication | if ($wantv) { |
|
| 303 | for ($i = 0; $i < $this->n; ++$i) { |
||
| 304 | $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1]; |
||
| 305 | $this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1]; |
||
| 306 | $this->V[$i][$j] = $t; |
||
| 307 | } |
||
| 308 | } |
||
| 309 | } |
||
| 310 | break; |
||
| 311 | // Split at negligible s(k). |
||
| 312 | case 2: |
||
| 313 | $f = $e[$k-1]; |
||
| 314 | $e[$k-1] = 0.0; |
||
| 315 | for ($j = $k; $j < $p; ++$j) { |
||
| 316 | $t = hypo($this->s[$j], $f); |
||
| 317 | $cs = $this->s[$j] / $t; |
||
| 318 | $sn = $f / $t; |
||
| 319 | $this->s[$j] = $t; |
||
| 320 | $f = -$sn * $e[$j]; |
||
| 321 | $e[$j] = $cs * $e[$j]; |
||
| 322 | View Code Duplication | if ($wantu) { |
|
| 323 | for ($i = 0; $i < $this->m; ++$i) { |
||
| 324 | $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1]; |
||
| 325 | $this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1]; |
||
| 326 | $this->U[$i][$j] = $t; |
||
| 327 | } |
||
| 328 | } |
||
| 329 | } |
||
| 330 | break; |
||
| 331 | // Perform one qr step. |
||
| 332 | case 3: |
||
| 333 | // Calculate the shift. |
||
| 334 | $scale = max(max(max(max( |
||
| 335 | abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])), |
||
| 336 | abs($this->s[$k])), abs($e[$k])); |
||
| 337 | $sp = $this->s[$p-1] / $scale; |
||
| 338 | $spm1 = $this->s[$p-2] / $scale; |
||
| 339 | $epm1 = $e[$p-2] / $scale; |
||
| 340 | $sk = $this->s[$k] / $scale; |
||
| 341 | $ek = $e[$k] / $scale; |
||
| 342 | $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0; |
||
| 343 | $c = ($sp * $epm1) * ($sp * $epm1); |
||
| 344 | $shift = 0.0; |
||
| 345 | if (($b != 0.0) || ($c != 0.0)) { |
||
| 346 | $shift = sqrt($b * $b + $c); |
||
| 347 | if ($b < 0.0) { |
||
| 348 | $shift = -$shift; |
||
| 349 | } |
||
| 350 | $shift = $c / ($b + $shift); |
||
| 351 | } |
||
| 352 | $f = ($sk + $sp) * ($sk - $sp) + $shift; |
||
| 353 | $g = $sk * $ek; |
||
| 354 | // Chase zeros. |
||
| 355 | for ($j = $k; $j < $p-1; ++$j) { |
||
| 356 | $t = hypo($f,$g); |
||
| 357 | $cs = $f/$t; |
||
| 358 | $sn = $g/$t; |
||
| 359 | if ($j != $k) { |
||
| 360 | $e[$j-1] = $t; |
||
| 361 | } |
||
| 362 | $f = $cs * $this->s[$j] + $sn * $e[$j]; |
||
| 363 | $e[$j] = $cs * $e[$j] - $sn * $this->s[$j]; |
||
| 364 | $g = $sn * $this->s[$j+1]; |
||
| 365 | $this->s[$j+1] = $cs * $this->s[$j+1]; |
||
| 366 | View Code Duplication | if ($wantv) { |
|
| 367 | for ($i = 0; $i < $this->n; ++$i) { |
||
| 368 | $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1]; |
||
| 369 | $this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1]; |
||
| 370 | $this->V[$i][$j] = $t; |
||
| 371 | } |
||
| 372 | } |
||
| 373 | $t = hypo($f,$g); |
||
| 374 | $cs = $f/$t; |
||
| 375 | $sn = $g/$t; |
||
| 376 | $this->s[$j] = $t; |
||
| 377 | $f = $cs * $e[$j] + $sn * $this->s[$j+1]; |
||
| 378 | $this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1]; |
||
| 379 | $g = $sn * $e[$j+1]; |
||
| 380 | $e[$j+1] = $cs * $e[$j+1]; |
||
| 381 | View Code Duplication | if ($wantu && ($j < $this->m - 1)) { |
|
| 382 | for ($i = 0; $i < $this->m; ++$i) { |
||
| 383 | $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1]; |
||
| 384 | $this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1]; |
||
| 385 | $this->U[$i][$j] = $t; |
||
| 386 | } |
||
| 387 | } |
||
| 388 | } |
||
| 389 | $e[$p-2] = $f; |
||
| 390 | $iter = $iter + 1; |
||
| 391 | break; |
||
| 392 | // Convergence. |
||
| 393 | case 4: |
||
| 394 | // Make the singular values positive. |
||
| 395 | if ($this->s[$k] <= 0.0) { |
||
| 396 | $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0); |
||
| 397 | View Code Duplication | if ($wantv) { |
|
| 398 | for ($i = 0; $i <= $pp; ++$i) { |
||
| 399 | $this->V[$i][$k] = -$this->V[$i][$k]; |
||
| 400 | } |
||
| 401 | } |
||
| 402 | } |
||
| 403 | // Order the singular values. |
||
| 404 | while ($k < $pp) { |
||
| 405 | if ($this->s[$k] >= $this->s[$k+1]) { |
||
| 406 | break; |
||
| 407 | } |
||
| 408 | $t = $this->s[$k]; |
||
| 409 | $this->s[$k] = $this->s[$k+1]; |
||
| 410 | $this->s[$k+1] = $t; |
||
| 411 | View Code Duplication | if ($wantv AND ($k < $this->n - 1)) { |
|
| 412 | for ($i = 0; $i < $this->n; ++$i) { |
||
| 413 | $t = $this->V[$i][$k+1]; |
||
| 414 | $this->V[$i][$k+1] = $this->V[$i][$k]; |
||
| 415 | $this->V[$i][$k] = $t; |
||
| 416 | } |
||
| 417 | } |
||
| 418 | View Code Duplication | if ($wantu AND ($k < $this->m-1)) { |
|
| 419 | for ($i = 0; $i < $this->m; ++$i) { |
||
| 420 | $t = $this->U[$i][$k+1]; |
||
| 421 | $this->U[$i][$k+1] = $this->U[$i][$k]; |
||
| 422 | $this->U[$i][$k] = $t; |
||
| 423 | } |
||
| 424 | } |
||
| 425 | ++$k; |
||
| 426 | } |
||
| 427 | $iter = 0; |
||
| 428 | --$p; |
||
| 429 | break; |
||
| 430 | } // end switch |
||
| 431 | } // end while |
||
| 432 | |||
| 433 | } // end constructor |
||
| 434 | |||
| 527 |