Passed
03:28
created

Math/LinearAlgebra/EigenvalueDecomposition.php (15 issues)

duplicate/similar code.

Duplication Informational

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); 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 136 \$sum += \$vect[\$i] ** 2; 137 } 138 139 \$sum = sqrt(\$sum); 140 View Code Duplication for (\$i = 0; \$i < count(\$vect); ++\$i) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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; 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) { 0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-12-01 10:32 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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) { 0 ignored issues – show Duplication introduced 2017-04-19 11:50 UTC by This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history... 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