Passed
03:28
created

#### Math/LinearAlgebra/EigenvalueDecomposition.php (3 issues)

Check for code that has been commented out.

Comprehensibility Unused Code Minor

#### Upgrade to new PHP Analysis Engine

These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more

 1 `A = \$Arg;` 99 ` \$this->n = count(\$Arg[0]);` 100 ` \$this->symmetric = true;` 101 102 ` for (\$j = 0; (\$j < \$this->n) && \$this->symmetric; ++\$j) {` 103 ` for (\$i = 0; (\$i < \$this->n) & \$this->symmetric; ++\$i) {` 104 ` \$this->symmetric = (\$this->A[\$i][\$j] == \$this->A[\$j][\$i]);` 105 ` }` 106 ` }` 107 108 ` if (\$this->symmetric) {` 109 ` \$this->V = \$this->A;` 110 ` // Tridiagonalize.` 111 ` \$this->tred2();` 112 ` // Diagonalize.` 113 ` \$this->tql2();` 114 ` } else {` 115 ` \$this->H = \$this->A;` 116 ` \$this->ort = [];` 117 ` // Reduce to Hessenberg form.` 118 ` \$this->orthes();` 119 ` // Reduce Hessenberg to real Schur form.` 120 ` \$this->hqr2();` 121 ` }` 122 ` }` 123 124 ` /**` 125 ` * Return the eigenvector matrix` 126 ` */` 127 ` public function getEigenvectors(): array` 128 ` {` 129 ` \$vectors = \$this->V;` 130 131 ` // Always return the eigenvectors of length 1.0` 132 ` \$vectors = new Matrix(\$vectors);` 133 ` \$vectors = array_map(function (\$vect) {` 134 ` \$sum = 0;` 135 View Code Duplication ` for (\$i = 0; \$i < count(\$vect); ++\$i) {` 136 ` \$sum += \$vect[\$i] ** 2;` 137 ` }` 138 139 ` \$sum = sqrt(\$sum);` 140 View Code Duplication ` for (\$i = 0; \$i < count(\$vect); ++\$i) {` 141 ` \$vect[\$i] /= \$sum;` 142 ` }` 143 144 ` return \$vect;` 145 ` }, \$vectors->transpose()->toArray());` 146 147 ` return \$vectors;` 148 ` }` 149 150 ` /**` 151 ` * Return the real parts of the eigenvalues
` 152 ` * d = real(diag(D));` 153 ` */` 154 ` public function getRealEigenvalues(): array` 155 ` {` 156 ` return \$this->d;` 157 ` }` 158 159 ` /**` 160 ` * Return the imaginary parts of the eigenvalues
` 161 ` * d = imag(diag(D))` 162 ` */` 163 ` public function getImagEigenvalues(): array` 164 ` {` 165 ` return \$this->e;` 166 ` }` 167 168 ` /**` 169 ` * Return the block diagonal eigenvalue matrix` 170 ` */` 171 ` public function getDiagonalEigenvalues(): array` 172 ` {` 173 ` \$D = [];` 174 175 ` for (\$i = 0; \$i < \$this->n; ++\$i) {` 176 ` \$D[\$i] = array_fill(0, \$this->n, 0.0);` 177 ` \$D[\$i][\$i] = \$this->d[\$i];` 178 ` if (\$this->e[\$i] == 0) {` 179 ` continue;` 180 ` }` 181 182 ` \$o = (\$this->e[\$i] > 0) ? \$i + 1 : \$i - 1;` 183 ` \$D[\$i][\$o] = \$this->e[\$i];` 184 ` }` 185 186 ` return \$D;` 187 ` }` 188 189 ` /**` 190 ` * Symmetric Householder reduction to tridiagonal form.` 191 ` */` 192 ` private function tred2(): void` 193 ` {` 194 ` // This is derived from the Algol procedures tred2 by` 195 ` // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for` 196 ` // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding` 197 ` // Fortran subroutine in EISPACK.` 198 ` \$this->d = \$this->V[\$this->n - 1];` 199 ` // Householder reduction to tridiagonal form.` 200 ` for (\$i = \$this->n - 1; \$i > 0; --\$i) {` 201 ` \$i_ = \$i - 1;` 202 ` // Scale to avoid under/overflow.` 203 ` \$h = \$scale = 0.0;` 204 ` \$scale += array_sum(array_map('abs', \$this->d));` 205 ` if (\$scale == 0.0) {` 206 ` \$this->e[\$i] = \$this->d[\$i_];` 207 ` \$this->d = array_slice(\$this->V[\$i_], 0, \$i_);` 208 ` for (\$j = 0; \$j < \$i; ++\$j) {` 209 ` \$this->V[\$j][\$i] = \$this->V[\$i][\$j] = 0.0;` 210 ` }` 211 ` } else {` 212 ` // Generate Householder vector.` 213 ` for (\$k = 0; \$k < \$i; ++\$k) {` 214 ` \$this->d[\$k] /= \$scale;` 215 ` \$h += pow(\$this->d[\$k], 2);` 216 ` }` 217 218 ` \$f = \$this->d[\$i_];` 219 ` \$g = sqrt(\$h);` 220 ` if (\$f > 0) {` 221 ` \$g = -\$g;` 222 ` }` 223 224 ` \$this->e[\$i] = \$scale * \$g;` 225 ` \$h = \$h - \$f * \$g;` 226 ` \$this->d[\$i_] = \$f - \$g;` 227 228 ` for (\$j = 0; \$j < \$i; ++\$j) {` 229 ` \$this->e[\$j] = 0.0;` 230 ` }` 231 232 ` // Apply similarity transformation to remaining columns.` 233 ` for (\$j = 0; \$j < \$i; ++\$j) {` 234 ` \$f = \$this->d[\$j];` 235 ` \$this->V[\$j][\$i] = \$f;` 236 ` \$g = \$this->e[\$j] + \$this->V[\$j][\$j] * \$f;` 237 238 ` for (\$k = \$j + 1; \$k <= \$i_; ++\$k) {` 239 ` \$g += \$this->V[\$k][\$j] * \$this->d[\$k];` 240 ` \$this->e[\$k] += \$this->V[\$k][\$j] * \$f;` 241 ` }` 242 243 ` \$this->e[\$j] = \$g;` 244 ` }` 245 246 ` \$f = 0.0;` 247 ` if (\$h === 0 || \$h < 1e-32) {` 248 ` \$h = 1e-32;` 249 ` }` 250 251 ` for (\$j = 0; \$j < \$i; ++\$j) {` 252 ` \$this->e[\$j] /= \$h;` 253 ` \$f += \$this->e[\$j] * \$this->d[\$j];` 254 ` }` 255 256 ` \$hh = \$f / (2 * \$h);` 257 ` for (\$j = 0; \$j < \$i; ++\$j) {` 258 ` \$this->e[\$j] -= \$hh * \$this->d[\$j];` 259 ` }` 260 261 ` for (\$j = 0; \$j < \$i; ++\$j) {` 262 ` \$f = \$this->d[\$j];` 263 ` \$g = \$this->e[\$j];` 264 ` for (\$k = \$j; \$k <= \$i_; ++\$k) {` 265 ` \$this->V[\$k][\$j] -= (\$f * \$this->e[\$k] + \$g * \$this->d[\$k]);` 266 ` }` 267 268 ` \$this->d[\$j] = \$this->V[\$i - 1][\$j];` 269 ` \$this->V[\$i][\$j] = 0.0;` 270 ` }` 271 ` }` 272 273 ` \$this->d[\$i] = \$h;` 274 ` }` 275 276 ` // Accumulate transformations.` 277 ` for (\$i = 0; \$i < \$this->n - 1; ++\$i) {` 278 ` \$this->V[\$this->n - 1][\$i] = \$this->V[\$i][\$i];` 279 ` \$this->V[\$i][\$i] = 1.0;` 280 ` \$h = \$this->d[\$i + 1];` 281 ` if (\$h != 0.0) {` 282 ` for (\$k = 0; \$k <= \$i; ++\$k) {` 283 ` \$this->d[\$k] = \$this->V[\$k][\$i + 1] / \$h;` 284 ` }` 285 286 ` for (\$j = 0; \$j <= \$i; ++\$j) {` 287 ` \$g = 0.0;` 288 ` for (\$k = 0; \$k <= \$i; ++\$k) {` 289 ` \$g += \$this->V[\$k][\$i + 1] * \$this->V[\$k][\$j];` 290 ` }` 291 292 ` for (\$k = 0; \$k <= \$i; ++\$k) {` 293 ` \$this->V[\$k][\$j] -= \$g * \$this->d[\$k];` 294 ` }` 295 ` }` 296 ` }` 297 298 ` for (\$k = 0; \$k <= \$i; ++\$k) {` 299 ` \$this->V[\$k][\$i + 1] = 0.0;` 300 ` }` 301 ` }` 302 303 ` \$this->d = \$this->V[\$this->n - 1];` 304 ` \$this->V[\$this->n - 1] = array_fill(0, \$j, 0.0);` 305 ` \$this->V[\$this->n - 1][\$this->n - 1] = 1.0;` 306 ` \$this->e[0] = 0.0;` 307 ` }` 308 309 ` /**` 310 ` * Symmetric tridiagonal QL algorithm.` 311 ` *` 312 ` * This is derived from the Algol procedures tql2, by` 313 ` * Bowdler, Martin, Reinsch, and Wilkinson, Handbook for` 314 ` * Auto. Comp., Vol.ii-Linear Algebra, and the corresponding` 315 ` * Fortran subroutine in EISPACK.` 316 ` */` 317 ` private function tql2(): void` 318 ` {` 319 ` for (\$i = 1; \$i < \$this->n; ++\$i) {` 320 ` \$this->e[\$i - 1] = \$this->e[\$i];` 321 ` }` 322 323 ` \$this->e[\$this->n - 1] = 0.0;` 324 ` \$f = 0.0;` 325 ` \$tst1 = 0.0;` 326 ` \$eps = pow(2.0, -52.0);` 327 328 ` for (\$l = 0; \$l < \$this->n; ++\$l) {` 329 ` // Find small subdiagonal element` 330 ` \$tst1 = max(\$tst1, abs(\$this->d[\$l]) + abs(\$this->e[\$l]));` 331 ` \$m = \$l;` 332 ` while (\$m < \$this->n) {` 333 ` if (abs(\$this->e[\$m]) <= \$eps * \$tst1) {` 334 ` break;` 335 ` }` 336 337 ` ++\$m;` 338 ` }` 339 340 ` // If m == l, \$this->d[l] is an eigenvalue,` 341 ` // otherwise, iterate.` 342 ` if (\$m > \$l) {` 343 ` \$iter = 0;` 344 ` do {` 345 ` // Could check iteration count here.` 346 ` \$iter += 1;` 347 ` // Compute implicit shift` 348 ` \$g = \$this->d[\$l];` 349 ` \$p = (\$this->d[\$l + 1] - \$g) / (2.0 * \$this->e[\$l]);` 350 ` \$r = hypot(\$p, 1.0);` 351 ` if (\$p < 0) {` 352 ` \$r *= -1;` 353 ` }` 354 355 ` \$this->d[\$l] = \$this->e[\$l] / (\$p + \$r);` 356 ` \$this->d[\$l + 1] = \$this->e[\$l] * (\$p + \$r);` 357 ` \$dl1 = \$this->d[\$l + 1];` 358 ` \$h = \$g - \$this->d[\$l];` 359 ` for (\$i = \$l + 2; \$i < \$this->n; ++\$i) {` 360 ` \$this->d[\$i] -= \$h;` 361 ` }` 362 363 ` \$f += \$h;` 364 ` // Implicit QL transformation.` 365 ` \$p = \$this->d[\$m];` 366 ` \$c = 1.0;` 367 ` \$c2 = \$c3 = \$c;` 368 ` \$el1 = \$this->e[\$l + 1];` 369 ` \$s = \$s2 = 0.0;` 370 ` for (\$i = \$m - 1; \$i >= \$l; --\$i) {` 371 ` \$c3 = \$c2;` 372 ` \$c2 = \$c;` 373 ` \$s2 = \$s;` 374 ` \$g = \$c * \$this->e[\$i];` 375 ` \$h = \$c * \$p;` 376 ` \$r = hypot(\$p, \$this->e[\$i]);` 377 ` \$this->e[\$i + 1] = \$s * \$r;` 378 ` \$s = \$this->e[\$i] / \$r;` 379 ` \$c = \$p / \$r;` 380 ` \$p = \$c * \$this->d[\$i] - \$s * \$g;` 381 ` \$this->d[\$i + 1] = \$h + \$s * (\$c * \$g + \$s * \$this->d[\$i]);` 382 ` // Accumulate transformation.` 383 ` for (\$k = 0; \$k < \$this->n; ++\$k) {` 384 ` \$h = \$this->V[\$k][\$i + 1];` 385 ` \$this->V[\$k][\$i + 1] = \$s * \$this->V[\$k][\$i] + \$c * \$h;` 386 ` \$this->V[\$k][\$i] = \$c * \$this->V[\$k][\$i] - \$s * \$h;` 387 ` }` 388 ` }` 389 390 ` \$p = -\$s * \$s2 * \$c3 * \$el1 * \$this->e[\$l] / \$dl1;` 391 ` \$this->e[\$l] = \$s * \$p;` 392 ` \$this->d[\$l] = \$c * \$p;` 393 ` // Check for convergence.` 394 ` } while (abs(\$this->e[\$l]) > \$eps * \$tst1);` 395 ` }` 396 397 ` \$this->d[\$l] = \$this->d[\$l] + \$f;` 398 ` \$this->e[\$l] = 0.0;` 399 ` }` 400 401 ` // Sort eigenvalues and corresponding vectors.` 402 ` for (\$i = 0; \$i < \$this->n - 1; ++\$i) {` 403 ` \$k = \$i;` 404 ` \$p = \$this->d[\$i];` 405 ` for (\$j = \$i + 1; \$j < \$this->n; ++\$j) {` 406 ` if (\$this->d[\$j] < \$p) {` 407 ` \$k = \$j;` 408 ` \$p = \$this->d[\$j];` 409 ` }` 410 ` }` 411 412 ` if (\$k != \$i) {` 413 ` \$this->d[\$k] = \$this->d[\$i];` 414 ` \$this->d[\$i] = \$p;` 415 View Code Duplication ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 416 ` \$p = \$this->V[\$j][\$i];` 417 ` \$this->V[\$j][\$i] = \$this->V[\$j][\$k];` 418 ` \$this->V[\$j][\$k] = \$p;` 419 ` }` 420 ` }` 421 ` }` 422 ` }` 423 424 ` /**` 425 ` * Nonsymmetric reduction to Hessenberg form.` 426 ` *` 427 ` * This is derived from the Algol procedures orthes and ortran,` 428 ` * by Martin and Wilkinson, Handbook for Auto. Comp.,` 429 ` * Vol.ii-Linear Algebra, and the corresponding` 430 ` * Fortran subroutines in EISPACK.` 431 ` */` 432 ` private function orthes(): void` 433 ` {` 434 ` \$low = 0;` 435 ` \$high = \$this->n - 1;` 436 437 ` for (\$m = \$low + 1; \$m <= \$high - 1; ++\$m) {` 438 ` // Scale column.` 439 ` \$scale = 0.0;` 440 View Code Duplication ` for (\$i = \$m; \$i <= \$high; ++\$i) {` 441 ` \$scale = \$scale + abs(\$this->H[\$i][\$m - 1]);` 442 ` }` 443 444 ` if (\$scale != 0.0) {` 445 ` // Compute Householder transformation.` 446 ` \$h = 0.0;` 447 ` for (\$i = \$high; \$i >= \$m; --\$i) {` 448 ` \$this->ort[\$i] = \$this->H[\$i][\$m - 1] / \$scale;` 449 ` \$h += \$this->ort[\$i] * \$this->ort[\$i];` 450 ` }` 451 452 ` \$g = sqrt(\$h);` 453 ` if (\$this->ort[\$m] > 0) {` 454 ` \$g *= -1;` 455 ` }` 456 457 ` \$h -= \$this->ort[\$m] * \$g;` 458 ` \$this->ort[\$m] -= \$g;` 459 ` // Apply Householder similarity transformation` 460 ` // H = (I -u * u' / h) * H * (I -u * u') / h)` 461 View Code Duplication ` for (\$j = \$m; \$j < \$this->n; ++\$j) {` 462 ` \$f = 0.0;` 463 ` for (\$i = \$high; \$i >= \$m; --\$i) {` 464 ` \$f += \$this->ort[\$i] * \$this->H[\$i][\$j];` 465 ` }` 466 467 ` \$f /= \$h;` 468 ` for (\$i = \$m; \$i <= \$high; ++\$i) {` 469 ` \$this->H[\$i][\$j] -= \$f * \$this->ort[\$i];` 470 ` }` 471 ` }` 472 473 View Code Duplication ` for (\$i = 0; \$i <= \$high; ++\$i) {` 474 ` \$f = 0.0;` 475 ` for (\$j = \$high; \$j >= \$m; --\$j) {` 476 ` \$f += \$this->ort[\$j] * \$this->H[\$i][\$j];` 477 ` }` 478 479 ` \$f = \$f / \$h;` 480 ` for (\$j = \$m; \$j <= \$high; ++\$j) {` 481 ` \$this->H[\$i][\$j] -= \$f * \$this->ort[\$j];` 482 ` }` 483 ` }` 484 485 ` \$this->ort[\$m] = \$scale * \$this->ort[\$m];` 486 ` \$this->H[\$m][\$m - 1] = \$scale * \$g;` 487 ` }` 488 ` }` 489 490 ` // Accumulate transformations (Algol's ortran).` 491 ` for (\$i = 0; \$i < \$this->n; ++\$i) {` 492 ` for (\$j = 0; \$j < \$this->n; ++\$j) {` 493 ` \$this->V[\$i][\$j] = (\$i == \$j ? 1.0 : 0.0);` 494 ` }` 495 ` }` 496 497 ` for (\$m = \$high - 1; \$m >= \$low + 1; --\$m) {` 498 ` if (\$this->H[\$m][\$m - 1] != 0.0) {` 499 View Code Duplication ` for (\$i = \$m + 1; \$i <= \$high; ++\$i) {` 500 ` \$this->ort[\$i] = \$this->H[\$i][\$m - 1];` 501 ` }` 502 503 ` for (\$j = \$m; \$j <= \$high; ++\$j) {` 504 ` \$g = 0.0;` 505 ` for (\$i = \$m; \$i <= \$high; ++\$i) {` 506 ` \$g += \$this->ort[\$i] * \$this->V[\$i][\$j];` 507 ` }` 508 509 ` // Double division avoids possible underflow` 510 ` \$g = (\$g / \$this->ort[\$m]) / \$this->H[\$m][\$m - 1];` 511 ` for (\$i = \$m; \$i <= \$high; ++\$i) {` 512 ` \$this->V[\$i][\$j] += \$g * \$this->ort[\$i];` 513 ` }` 514 ` }` 515 ` }` 516 ` }` 517 ` }` 518 519 ` /**` 520 ` * Performs complex division.` 521 ` *` 522 ` * @param int|float \$xr` 523 ` * @param int|float \$xi` 524 ` * @param int|float \$yr` 525 ` * @param int|float \$yi` 526 ` */` 527 ` private function cdiv(\$xr, \$xi, \$yr, \$yi): void` 528 ` {` 529 ` if (abs(\$yr) > abs(\$yi)) {` 530 ` \$r = \$yi / \$yr;` 531 ` \$d = \$yr + \$r * \$yi;` 532 ` \$this->cdivr = (\$xr + \$r * \$xi) / \$d;` 533 ` \$this->cdivi = (\$xi - \$r * \$xr) / \$d;` 534 ` } else {` 535 ` \$r = \$yr / \$yi;` 536 ` \$d = \$yi + \$r * \$yr;` 537 ` \$this->cdivr = (\$r * \$xr + \$xi) / \$d;` 538 ` \$this->cdivi = (\$r * \$xi - \$xr) / \$d;` 539 ` }` 540 ` }` 541 542 ` /**` 543 ` * Nonsymmetric reduction from Hessenberg to real Schur form.` 544 ` *` 545 ` * Code is derived from the Algol procedure hqr2,` 546 ` * by Martin and Wilkinson, Handbook for Auto. Comp.,` 547 ` * Vol.ii-Linear Algebra, and the corresponding` 548 ` * Fortran subroutine in EISPACK.` 549 ` */` 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)` 0 ignored issues – show Unused Code Comprehensibility introduced 2017-04-19 11:50 UTC by `56%` of this comment could be valid code. Did you maybe forget this after debugging? Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it. The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production. This check looks for comments that seem to be mostly valid code and reports them. Loading history... 813 ` } // k loop` 814 ` } // check convergence` 815 ` } // while (\$n >= \$low)` 0 ignored issues – show Unused Code Comprehensibility introduced 2017-04-19 11:50 UTC by `46%` of this comment could be valid code. Did you maybe forget this after debugging? Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it. The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production. This check looks for comments that seem to be mostly valid code and reports them. Loading history... 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;` 0 ignored issues – show Unused Code Comprehensibility introduced 2017-04-19 11:50 UTC by `37%` of this comment could be valid code. Did you maybe forget this after debugging? Sometimes obsolete code just ends up commented out instead of removed. In this case it is better to remove the code once you have checked you do not need it. The code might also have been commented out for debugging purposes. In this case it is vital that someone uncomments it again or your project may behave in very unexpected ways in production. This check looks for comments that seem to be mostly valid code and reports them. Loading history... 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