001/******************************************************************************* 002 * This software is provided as a supplement to the authors' textbooks on digital 003 * image processing published by Springer-Verlag in various languages and editions. 004 * Permission to use and distribute this software is granted under the BSD 2-Clause 005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 006 * Copyright (c) 2006-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.common.math.eigen.eispack; 010 011public abstract class QZIT { 012 013 /* 014 EISPACK Routines, see http://www.netlib.org/eispack/ for the original FORTRAN code. 015 Untangled to goto-free Java by W. Burger using a sequential state machine concept, inspired by D. E. Knuth, 016 "Structured Programming with Goto Statements", Computing Surveys, Vol. 6, No. 4 (1974). 017 */ 018 019 private QZIT() {} 020 021 /** 022 * <p> 023 * This subroutine is the second step of the qz algorithm for solving generalized matrix eigenvalue problems, Siam 024 * J. Numer. anal. 10, 241-256(1973) by Moler and Stewart, as modified in technical note NASA TN D-7305(1973) by 025 * Ward. This description has been adapted from the original version (dated August 1983). 026 * </p> 027 * <p> 028 * This subroutine accepts a pair of real matrices, one of them in upper Hessenberg form and the other in upper 029 * triangular form. It reduces the Hessenberg matrix to quasi-triangular form using orthogonal transformations while 030 * maintaining the triangular form of the other matrix. It is usually preceded by qzhes and followed by qzval and, 031 * possibly, qzvec. 032 * </p> 033 * <p> 034 * On output: 035 * </p> 036 * <ul> 037 * <li><strong>a</strong> has been reduced to quasi-triangular form. The elements below the first subdiagonal are 038 * still zero and no two consecutive subdiagonal elements are nonzero.</li> 039 * <li><strong>b</strong> is still in upper triangular form, although its elements have been altered. The location 040 * b(n,1) is used to store eps1 times the norm of b for later use by qzval and qzvec. 041 * <li><strong>z</strong> contains the product of the right hand transformations (for both steps) if matz has been 042 * set to true.</li> 043 * </ul> 044 * 045 * @param a contains a real upper hessenberg matrix. 046 * @param b contains a real upper triangular matrix. 047 * @param eps1 is a tolerance used to determine negligible elements. eps1 = 0.0 (or negative) may be input, in which 048 * case an element will be neglected only if it is less than roundoff error times the norm of its matrix. If the 049 * input eps1 is positive, then an element will be considered negligible if it is less than eps1 times the norm of 050 * its matrix. A positive value of eps1 may result in faster execution, but less accurate results. 051 * @param matz should be set to true if the right hand transformations are to be accumulated for later use in 052 * computing eigenvectors, and to false otherwise. 053 * @param z contains, if matz has been set to true, the transformation matrix produced in the reduction by qzhes, if 054 * performed, or else the identity matrix. If matz has been set to false, z is not referenced. 055 * @return -1 for normal return, j if the limit of 30*n iterations is exhausted while the j-th eigenvalue is being 056 * is being sought. 057 */ 058 public static int qzit(double[][] a, double[][] b, double eps1, boolean matz, double[][] z) { 059 060 final int n = a.length; 061 062 int i, j, k, l = 0; 063 double r, s, t, a1 = 0, a2 = 0, a3 = 0; 064 int k1, k2, l1 = 0, ll; 065 double u1, u2, u3; 066 double v1, v2, v3; 067 double a11 = 0, a12, a21 = 0, a22, a33 = 0, a34 = 0, a43 = 0, a44 = 0; 068 double b11 = 0, b12, b22 = 0, b33, b34 = 0, b44; 069 int na = 0, en, ld = 0; 070 double ep; 071 double sh = 0; 072 int km1, lm1 = 0; 073 double ani, bni; 074 int ish = 0, itn, its = 0, enm2 = 0, lor1; 075 double epsa, epsb, anorm, bnorm; 076 int enorn; 077 boolean notlas; 078 079 int ierr = -1; // no error 080 081 // compute epsa and epsb 082 anorm = 0.0; 083 bnorm = 0.0; 084 085 for (i = 0; i < n; i++) { 086 ani = 0.0; 087 if (i != 0) { 088 ani = (Math.abs(a[i][(i - 1)])); 089 } 090 bni = 0.0; 091 092 for (j = i; j < n; j++) { 093 ani += Math.abs(a[i][j]); 094 bni += Math.abs(b[i][j]); 095 } 096 097 if (ani > anorm) { 098 anorm = ani; 099 } 100 if (bni > bnorm) { 101 bnorm = bni; 102 } 103 } 104 105 if (anorm == 0.0) { 106 anorm = 1.0; 107 } 108 if (bnorm == 0.0) { 109 bnorm = 1.0; 110 } 111 112// ep = eps1; 113// if (ep == 0.0) { 114// // use round-off level if eps1 is zero 115// ep = epsilon(1.0); 116// } 117 ep = (eps1 > 0) ? eps1 : epslon(1.0); 118 epsa = ep * anorm; 119 epsb = ep * bnorm; 120 121 // rReduce a to quasi-triangular form, while keeping b triangular 122 lor1 = 0; 123 enorn = n; 124 en = n - 1; 125 itn = n * 30; 126 127 // ------------------------------------------ 128 129 State state = State.L60; 130 StateLoop: while (state != State.Final) { 131 switch (state) { 132 133 case L60: 134 // Begin QZ step 135 if (en <= 1) { 136 state = State.L1001; 137 break; 138 } 139 if (!matz) { 140 enorn = en + 1; 141 } 142 its = 0; 143 na = en - 1; 144 enm2 = na; 145 state = State.L70; 146 break; 147 148 case L70: 149 ish = 2; 150 // check for convergence or reducibility. 151 for (ll = 0; ll <= en; ++ll) { 152 lm1 = en - ll - 1; 153 l = lm1 + 1; 154 if (l == 0) { 155 state = State.L95; 156 continue StateLoop; // check again!! 157 } 158 if (Math.abs(a[l][lm1]) <= epsa) { 159 state = State.L90; 160 break; 161 } 162 } 163 state = State.L90; 164 break; 165 166 case L90: 167 a[l][lm1] = 0.0; 168 if (l < na) { 169 state = State.L95; 170 break; 171 } 172 // 1-by-1 or 2-by-2 block isolated 173 en = lm1; 174 state = State.L60; 175 break; 176 177 case L95: 178 // check for small top of b 179 ld = l; 180 state = State.L100; 181 break; 182 183 case L100: 184 l1 = l + 1; 185 b11 = b[l][l]; 186 187 if (Math.abs(b11) > epsb) { 188 state = State.L120; 189 break; 190 } 191 192 b[l][l] = 0.0; 193 s = Math.abs(a[l][l]) + Math.abs(a[l1][l]); 194 u1 = a[l][l] / s; 195 u2 = a[l1][l] / s; 196// r = Special.Sign(Math.sqrt(u1 * u1 + u2 * u2), u1); 197 r = Math.copySign(Math.hypot(u1, u2), u1); 198 v1 = -(u1 + r) / r; 199 v2 = -u2 / r; 200 u2 = v2 / v1; 201 202 for (j = l; j < enorn; j++) { 203 t = a[l][j] + u2 * a[l1][j]; 204 a[l][j] = a[l][j] + t * v1; 205 a[l1][j] = a[l1][j] + t * v2; 206 207 t = b[l][j] + u2 * b[l1][j]; 208 b[l][j] = b[l][j] + t * v1; 209 b[l1][j] = b[l1][j] + t * v2; 210 } 211 212 if (l != 0) { 213 a[l][lm1] = -a[l][lm1]; 214 } 215 lm1 = l; 216 l = l1; 217 state = State.L90; 218 break; 219 220 case L120: 221 a11 = a[l][l] / b11; 222 a21 = a[l1][l] / b11; 223 if (ish == 1) { 224 state = State.L140; 225 break; 226 } 227 // iteration strategy 228 if (itn == 0) { 229 state = State.L1000; 230 break; 231 } 232 if (its == 10) { 233 state = State.L155; 234 break; 235 } 236 237 // determine type of shift 238 b22 = b[l1][l1]; 239 if (Math.abs(b22) < epsb) { 240 b22 = epsb; 241 } 242 b33 = b[na][na]; 243 if (Math.abs(b33) < epsb) { 244 b33 = epsb; 245 } 246 b44 = b[en][en]; 247 if (Math.abs(b44) < epsb) 248 b44 = epsb; 249 a33 = a[na][na] / b33; 250 a34 = a[na][en] / b44; 251 a43 = a[en][na] / b33; 252 a44 = a[en][en] / b44; 253 b34 = b[na][en] / b44; 254 t = 0.5 * (a43 * b34 - a33 - a44); 255 r = t * t + a34 * a43 - a33 * a44; 256 if (r < 0.0) { 257 state = State.L150; 258 break; 259 } 260 261 // determine single shift zero-th column of a 262 ish = 1; 263 r = Math.sqrt(r); 264 sh = -t + r; 265 s = -t - r; 266 if (Math.abs(s - a44) < Math.abs(sh - a44)) { 267 sh = s; 268 } 269 270 // look for two consecutive small sub-diagonal elements of a. 271 for (ll = ld; ll < enm2; ll++) { 272 l = enm2 + ld - ll - 1; 273 if (l == ld) { 274 state = State.L140; 275 //break loop130; 276 continue StateLoop; 277 } 278 lm1 = l - 1; 279 l1 = l + 1; 280// t = a[l + 1][l + 1]; // CHECK THIS!! 281 t = a[l][l]; 282 283 if (Math.abs(b[l][l]) > epsb) { 284 t = t - sh * b[l][l]; 285 } 286 287 if (Math.abs(a[l][lm1]) <= Math.abs(t / a[l1][l]) * epsa) { 288 state = State.L100; 289 continue StateLoop; 290 } 291 } 292 state = State.L140; 293 break; 294 295 case L140: 296 a1 = a11 - sh; 297 a2 = a21; 298 if (l != ld) { 299 a[l][lm1] = -a[l][lm1]; 300 } 301 state = State.L160; 302 break; 303 304 case L150: 305 // determine double shift zero-th column of a 306 a12 = a[l][l1] / b22; 307 a22 = a[l1][l1] / b22; 308 b12 = b[l][l1] / b22; 309 a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 + a12 - a11 * b12; 310 a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34; 311 a3 = a[l1 + 1][l1] / b22; 312 state = State.L160; 313 break; 314 315 case L155: 316 // ad hoc shift 317 a1 = 0.0; 318 a2 = 1.0; 319 a3 = 1.1605; // magic! 320 state = State.L160; 321 break; 322 323 case L160: 324 its = its + 1; 325 itn = itn - 1; 326 if (!matz) { 327 lor1 = ld; 328 } 329 330 mainloop: for (k = l; k <= na; k++) { 331 notlas = (k != na && ish == 2); 332 k1 = k + 1; 333 k2 = k + 2; 334 km1 = Math.max(k, l + 1) - 1; 335 ll = Math.min(en, k1 + ish); 336 337 if (!notlas) { 338 // zero a(k+1,k-1) 339 if (k != l) { 340 a1 = a[k][km1]; 341 a2 = a[k1][km1]; 342 } 343 344 s = Math.abs(a1) + Math.abs(a2); 345 if (s == 0.0) { 346 state = State.L70; 347 continue StateLoop; 348 } 349 u1 = a1 / s; 350 u2 = a2 / s; 351 r = Math.copySign(Math.hypot(u1, u2), u1); 352 v1 = -(u1 + r) / r; 353 v2 = -u2 / r; 354 u2 = v2 / v1; 355 356 for (j = km1; j < enorn; j++) { 357 t = a[k][j] + u2 * a[k1][j]; 358 a[k][j] = a[k][j] + t * v1; 359 a[k1][j] = a[k1][j] + t * v2; 360 361 t = b[k][j] + u2 * b[k1][j]; 362 b[k][j] = b[k][j] + t * v1; 363 b[k1][j] = b[k1][j] + t * v2; 364 } 365 366 if (k != l) { 367 a[k1][km1] = 0.0; 368 } 369 // goto L240; 370 } 371 372 else { 373 // zero a(k+1,k-1) and a(k+2,k-1) 374 //L190: 375 if (k != l) { 376 a1 = a[k][km1]; 377 a2 = a[k1][km1]; 378 a3 = a[k2][km1]; 379 } 380 //L200: 381 s = Math.abs(a1) + Math.abs(a2) + Math.abs(a3); 382 if (s == 0.0) { 383 continue mainloop; // goto L260; 384 } 385 u1 = a1 / s; 386 u2 = a2 / s; 387 u3 = a3 / s; 388 r = Math.copySign(Math.sqrt(u1 * u1 + u2 * u2 + u3 * u3), u1); 389 v1 = -(u1 + r) / r; 390 v2 = -u2 / r; 391 v3 = -u3 / r; 392 u2 = v2 / v1; 393 u3 = v3 / v1; 394 395 for (j = km1; j < enorn; j++) { 396 t = a[k][j] + u2 * a[k1][j] + u3 * a[k2][j]; 397 a[k][j] = a[k][j] + t * v1; 398 a[k1][j] = a[k1][j] + t * v2; 399 a[k2][j] = a[k2][j] + t * v3; 400 401 t = b[k][j] + u2 * b[k1][j] + u3 * b[k2][j]; 402 b[k][j] = b[k][j] + t * v1; 403 b[k1][j] = b[k1][j] + t * v2; 404 b[k2][j] = b[k2][j] + t * v3; 405 } 406 407 if (k != l) { 408 a[k1][km1] = 0.0; 409 a[k2][km1] = 0.0; 410 } 411 412 // zero b(k+2,k+1) and b(k+2,k) 413 s = Math.abs(b[k2][k2]) + Math.abs(b[k2][k1]) + Math.abs(b[k2][k]); 414 if (s != 0.0) { 415 u1 = b[k2][k2] / s; 416 u2 = b[k2][k1] / s; 417 u3 = b[k2][k] / s; 418 r = Math.copySign(Math.sqrt(u1 * u1 + u2 * u2 + u3 * u3), u1); 419 v1 = -(u1 + r) / r; 420 v2 = -u2 / r; 421 v3 = -u3 / r; 422 u2 = v2 / v1; 423 u3 = v3 / v1; 424 425 for (i = lor1; i < ll + 1; i++) { 426 t = a[i][k2] + u2 * a[i][k1] + u3 * a[i][k]; 427 a[i][k2] = a[i][k2] + t * v1; 428 a[i][k1] = a[i][k1] + t * v2; 429 a[i][k] = a[i][k] + t * v3; 430 431 t = b[i][k2] + u2 * b[i][k1] + u3 * b[i][k]; 432 b[i][k2] = b[i][k2] + t * v1; 433 b[i][k1] = b[i][k1] + t * v2; 434 b[i][k] = b[i][k] + t * v3; 435 } 436 437 b[k2][k] = 0.0; 438 b[k2][k1] = 0.0; 439 440 if (matz) { 441 for (i = 0; i < n; i++) { 442 t = z[i][k2] + u2 * z[i][k1] + u3 * z[i][k]; 443 z[i][k2] = z[i][k2] + t * v1; 444 z[i][k1] = z[i][k1] + t * v2; 445 z[i][k] = z[i][k] + t * v3; 446 } 447 } 448 } 449 } 450 451 // L240: 452 // zero b(k+1,k) 453 s = (Math.abs(b[k1][k1])) + (Math.abs(b[k1][k])); 454 if (s == 0.0) { 455 continue; // goto L260; 456 } 457 u1 = b[k1][k1] / s; 458 u2 = b[k1][k] / s; 459 r = Math.copySign(Math.hypot(u1, u2), u1); 460 v1 = -(u1 + r) / r; 461 v2 = -u2 / r; 462 u2 = v2 / v1; 463 464 for (i = lor1; i <= ll; i++) { 465 t = a[i][k1] + u2 * a[i][k]; 466 a[i][k1] = a[i][k1] + t * v1; 467 a[i][k] = a[i][k] + t * v2; 468 469 t = b[i][k1] + u2 * b[i][k]; 470 b[i][k1] = b[i][k1] + t * v1; 471 b[i][k] = b[i][k] + t * v2; 472 } 473 474 b[k1][k] = 0.0; 475 476 if (matz) { 477 for (i = 0; i < n; i++) { 478 t = z[i][k1] + u2 * z[i][k]; 479 z[i][k1] = z[i][k1] + t * v1; 480 z[i][k] = z[i][k] + t * v2; 481 } 482 } 483 // L260: ; 484 } 485 486 state = State.L70; // end qz step 487 break; 488 489 case L1000: 490 // set error -- all eigenvalues have not converged after 30*n iterations 491 ierr = en + 1; 492 state = State.L1001; 493 break; 494 495 case L1001: 496 // save epsb for use by qzval and qzvec 497 if (n > 1) { 498 b[n - 1][0] = epsb; 499 } 500 state = State.Final; 501 break; 502 503 case Final: 504 throw new RuntimeException("this should never happen!"); 505 } 506 } 507 508 return ierr; // return eigenvalue index if iteration count exceeded 509 } // end of qzit() 510 511 512 private enum State { 513 L60, L70, L90, L95, L100, L120, L140, L150, L155, L160, L1000, L1001, Final; 514 } 515 516 private static double EpsDouble; 517 static { 518 double eps; 519 double a = 4.0 / 3.0; 520 do { 521 double b = a - 1.0; 522 double c = b + b + b; 523 eps = Math.abs(c - 1.0); 524 } while (eps == 0); 525 EpsDouble = eps; 526 } 527 528 private static float EpsFloat; 529 static { 530 float eps; 531 float a = 4.0f / 3.0f; 532 do { 533 float b = a - 1.0f; 534 float c = b + b + b; 535 eps = Math.abs(c - 1.0f); 536 } while (eps == 0); 537 EpsFloat = eps; 538 } 539 540 /** 541 * Estimates unit round-off in quantities of size x. This is a port of the epsilon function from EISPACK. See also 542 * https://www.researchgate.net/publication/2860253_A_Comment_on_the_Eispack_Machine_Epsilon_Routine 543 * 544 * @param x some quantity 545 * @return returns the smallest number e of the same kind as x such that 1 + e > 1. 546 */ 547 private static double epslon(double x) { 548// double eps; 549// double a = 4.0 / 3.0; 550// do { 551// double b = a - 1.0; 552// double c = b + b + b; 553// eps = Math.abs(c - 1.0); 554// } while (eps == 0); 555 return EpsDouble * Math.abs(x); 556 } 557 558 private static float epslon(float x) { 559 return EpsFloat * Math.abs(x); 560 561 } 562 563/* original FORTRAN code: 564 a = 4.0d0/3.0d0 565 10 b = a - 1.0d0 566 c = b + b + b 567 eps = dabs(c-1.0d0) 568 if (eps .eq. 0.0d0) go to 10 569 epslon = eps*dabs(x) 570 return 571 end 572*/ 573 574/* 575program test_epsilon 576 real :: x = 3.143 577 real(8) :: y = 2.33 578 print *, epsilon(x) 579 print *, epsilon(y) 580end program test_epsilon 581 */ 582 583 // public static void main(String[] args) { 584 // double x = 3.143; 585 // System.out.println(epslon(x)); 586 // float y = 2.33f; 587 // System.out.println(epslon(y)); 588 // 589 // } 590 591}