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 QZVEC { 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 QZVEC() {} 020 021 public static boolean VERBOSE = false; 022 023 /** 024 * <p> 025 * This subroutine is the optional fourth step of the qz algorithm for solving generalized matrix eigenvalue 026 * problems, Siam J. Numer. Anal. 10, 241-256 (1973) by Moler and Stewart. This description has been adapted from 027 * the original version (dated August 1983). 028 * </p> 029 * <p> 030 * This subroutine accepts a pair of real matrices, one of them in quasi-triangular form (in which each 2-by-2 block 031 * corresponds to a pair of complex eigenvalues) and the other in upper triangular form. It computes the 032 * eigenvectors of the triangular problem and transforms the results back to the original coordinate system. It is 033 * usually preceded by qzhes, qzit, and qzval. 034 * </p> 035 * <p> 036 * On output: 037 * </p> 038 * <ul> 039 * <li><strong>a</strong> is unaltered. Its subdiagonal elements provide information about the storage of the 040 * complex eigenvectors.</li> 041 * <li><strong>b</strong> has been destroyed.</li> 042 * <li><strong>alfr</strong>, <strong>alfi</strong> and <strong>beta</strong> are unaltered.</li> 043 * <li><strong>z</strong> contains the real and imaginary parts of the eigenvectors. If alfi(i) = 0.0, the i-th 044 * eigenvalue is real and the i-th column of z contains its eigenvector. if alfi(i) ≠ 0.0, the i-th eigenvalue 045 * is complex. if alfi(i) > 0.0, the eigenvalue is the first of a complex pair and the i-th and (i+1)-th columns 046 * of z contain its eigenvector. If alfi(i) < 0.0, the eigenvalue is the second of a complex pair and the 047 * (i-1)-th and i-th columns of z contain the conjugate of its eigenvector. Each eigenvector is normalized so that 048 * the modulus of its largest component is 1.0.</li> 049 * </ul> 050 * 051 * @param a contains a real upper quasi-triangular matrix. 052 * @param b contains a real upper triangular matrix. In addition, location b[n-1][0] contains the tolerance quantity 053 * (epsb) computed and saved in qzit. 054 * @param alfr , <strong>alfi</strong>, and <strong>beta</strong> are vectors with components whose ratios ((alfr + 055 * i * alfi) / beta) are the generalized eigenvalues. They are usually obtained from qzval. 056 * @param alfi see <strong>alfr</strong> 057 * @param beta see <strong>alfr</strong> 058 * @param z contains the transformation matrix produced in the reductions by qzhes, qzit, and qzval, if performed. 059 * If the eigenvectors of the triangular problem are desired, z must contain the identity matrix. 060 */ 061 public static void qzvec(double[][] a, double[][] b, double[] alfr, double[] alfi, double[] beta, double[][] z) { 062 063 final int n = a.length; 064 065 int i = 0, j, k, m; 066 int na = 0, ii, en = 0, jj, nn, enm2; 067 double d = 0, q; 068 double r = 0, s = 0, t = 0, w = 0, x = 0, y = 0, t1 = 0, t2 = 0, w1 = 0, x1 = 0, z1 = 0, di = 0; 069 double ra = 0, dr = 0, sa = 0; 070 double ti = 0, rr, tr = 0, zz = 0; 071 double alfm, almi, betm, almr; 072 073 double epsb = b[n - 1][0]; 074 int isw = 1; 075 076 if (VERBOSE) System.out.println("QZVEC: starting Loop A ***********************"); 077 LoopA: for (nn = 0; nn < n; nn++) { 078 079 StateA stateA = StateA.Ainit; 080 StateLoopA: while (stateA != StateA.Afinal) { 081 082 switch(stateA) { 083 case Ainit : 084 en = n - 1 - nn; // TODO: check!! 085 na = en - 1; 086 if (isw == 2) { 087 stateA = StateA.A795; 088 break; 089 } 090 if (alfi[en] != 0.0) { 091 stateA = StateA.A710; 092 break; 093 } 094 095 // Real vector 096 m = en; 097 b[en][en] = 1.0; 098 if (na == -1) { 099 stateA = StateA.A800; 100 break; 101 } 102 alfm = alfr[m]; 103 betm = beta[m]; 104 105 if (VERBOSE) System.out.println("QZVEC: starting Loop B ***********************"); 106 107 // for i=en-1 step -1 until 1 do -- 108 LoopB: for (ii = 0; ii <= na; ii++) { 109 110 StateB stateB = StateB.Binit; 111 StateLoopB: while (stateB != StateB.Bfinal) { 112 switch (stateB) { 113 case Binit : 114 i = en - ii - 1; 115 w = betm * a[i][i] - alfm * b[i][i]; 116 r = 0.0; 117 118 for (j = m; j <= en; j++) { 119 r = r +(betm * a[i][j] - alfm * b[i][j]) * b[j][en]; 120 } // L610 121 122 if (i == 0 || isw == 2) { 123 stateB = StateB.B630; 124 break; 125 } 126 if (betm * a[i][i - 1] == 0.0) { 127 stateB = StateB.B630; 128 break; 129 } 130 zz = w; 131 s = r; 132 stateB = StateB.B690; 133 break; 134 135 case B630 : 136 m = i; 137 if (isw == 2) { 138 stateB = StateB.B640; 139 break; 140 } 141 // real 1-by-1 block 142 t = w; 143 if (w == 0.0) { 144 t = epsb; 145 } 146 b[i][en] = -r / t; 147 stateB = StateB.B700; 148 break; 149 150 case B640 : 151 // real 2-by-2 block 152 x = betm * a[i][i + 1] - alfm * b[i][i + 1]; 153 y = betm * a[i + 1][i]; 154 q = w * zz - x * y; 155 t = (x * s - zz * r) / q; 156 b[i][en] = t; 157 if (Math.abs(x) <= Math.abs(zz)) { 158 stateB = StateB.B650; 159 break; 160 } 161 b[i + 1][en] = (-r - w * t) / x; 162 stateB = StateB.B690; 163 break; 164 165 case B650 : 166 b[i + 1][en] = (-s - y * t) / zz; 167 stateB = StateB.B690; 168 break; 169 170 case B690 : 171 isw = 3 - isw; 172 stateB = StateB.B700; 173 break; 174 175 case B700 : 176 stateB = StateB.Bfinal; 177 break; 178 179 case Bfinal : 180 throw new RuntimeException("this should never happen!"); 181 182 } // end switch (stateB) 183 } // end StateLoopB 184 } // L700: end of LoopB 185 if (VERBOSE) System.out.println("QZVEC: finished Loop B ***********************"); 186 // end real vector 187 stateA = StateA.A800; 188 break; 189 190 case A710: 191 // complex vector 192 m = na; 193 almr = alfr[m]; 194 almi = alfi[m]; 195 betm = beta[m]; 196 // last vector component chosen imaginary so that eigenvector matrix is triangular 197 y = betm * a[en][na]; 198 b[na][na] = -almi * b[en][en] / y; 199 b[na][en] = (almr * b[en][en] - betm * a[en][en]) / y; 200 b[en][na] = 0.0; 201 b[en][en] = 1.0; 202 enm2 = na; // TODO: check! 203 if (enm2 == 0) { 204 stateA = StateA.A795; 205 break; 206 } 207 208 if (VERBOSE) System.out.println("QZVEC: starting Loop C, enm2 = " + enm2); 209 LoopC: for (ii = 0; ii < enm2; ii++) { 210 if (VERBOSE) System.out.println(" ii = " + ii); 211 StateC stateC = StateC.Cinit; 212 StateLoopC: while (stateC != StateC.Cfinal) { 213 if (VERBOSE) System.out.println(" StateLoopC: " + stateC); 214 switch (stateC) { 215 216 case Cinit : 217 i = na - ii - 1; 218 w = betm * a[i][i] - almr * b[i][i]; 219 w1 = -almi * b[i][i]; 220 ra = 0.0; 221 sa = 0.0; 222 223 for (j = m; j <= en; j++) 224 { 225 x = betm * a[i][j] - almr * b[i][j]; 226 x1 = -almi * b[i][j]; 227 ra = ra + x * b[j][na] - x1 * b[j][en]; 228 sa = sa + x * b[j][en] + x1 * b[j][na]; 229 } // L760 230 231 if (i == 0 || isw == 2) { 232 stateC = StateC.C770; 233 break; 234 } 235 if (betm * a[i][i - 1] == 0.0) { 236 stateC = StateC.C770; 237 break; 238 } 239 zz = w; 240 z1 = w1; 241 r = ra; 242 s = sa; 243 isw = 2; 244 stateC = StateC.C790; 245 break; 246 247 case C770 : 248 m = i; 249 if (isw == 2) { 250 stateC = StateC.C780; 251 break; 252 } 253 // complex 1-by-1 block 254 tr = -ra; 255 ti = -sa; 256 stateC = StateC.C773; 257 break; 258 259 case C773: 260 dr = w; 261 di = w1; 262 stateC = StateC.C775; 263 break; 264 265 case C775: 266 // complex divide (t1,t2) = (tr,ti) / (dr,di) 267 if (Math.abs(di) > Math.abs(dr)) { 268 stateC = StateC.C777; 269 break; 270 } 271 rr = di / dr; 272 d = dr + di * rr; 273 t1 = (tr + ti * rr) / d; 274 t2 = (ti - tr * rr) / d; 275 276 switch (isw) { // go to (787,782), isw 277 case 1: 278 stateC = StateC.C787; 279 continue StateLoopC; 280 //break; 281 case 2: 282 stateC = StateC.C782; 283 continue StateLoopC; 284 //break; 285 default: 286 throw new RuntimeException("illegal state isw = " + isw); 287 } 288 //stateC = StateC.C777; 289 //break; 290 291 case C777: 292 rr = dr / di; 293 d = dr * rr + di; 294 t1 = (tr * rr + ti) / d; 295 t2 = (ti * rr - tr) / d; 296 switch (isw) { // go to (787,782), isw 297 case 1: 298 stateC = StateC.C787; 299 continue StateLoopC; 300 //break; 301 case 2: 302 stateC = StateC.C782; 303 continue StateLoopC; 304 //break; 305 default: 306 throw new RuntimeException("illegal state isw = " + isw); 307 } 308 //stateC = StateC.C780; 309 //break; 310 311 case C780: 312 // complex 2-by-2 block 313 x = betm * a[i][i + 1] - almr * b[i][i + 1]; 314 x1 = -almi * b[i][i + 1]; 315 y = betm * a[i + 1][i]; 316 tr = y * ra - w * r + w1 * s; 317 ti = y * sa - w * s - w1 * r; 318 dr = w * zz - w1 * z1 - x * y; 319 di = w * z1 + w1 * zz - x1 * y; 320 if (dr == 0.0 && di == 0.0) { 321 dr = epsb; 322 } 323 stateC = StateC.C775; 324 break; 325 326 case C782: 327 b[i + 1][na] = t1; 328 b[i + 1][en] = t2; 329 isw = 1; 330 if (Math.abs(y) > Math.abs(w) + Math.abs(w1)) { 331 stateC = StateC.C785; 332 break; 333 } 334 tr = -ra - x * b[(i + 1)][na] + x1 * b[(i + 1)][en]; 335 ti = -sa - x * b[(i + 1)][en] - x1 * b[(i + 1)][na]; 336 stateC = StateC.C773; 337 break; 338 339 case C785: 340 t1 = (-r - zz * b[(i + 1)][na] + z1 * b[(i + 1)][en]) / y; 341 t2 = (-s - zz * b[(i + 1)][en] - z1 * b[(i + 1)][na]) / y; 342 stateC = StateC.C787; 343 break; 344 345 case C787: 346 b[i][na] = t1; 347 b[i][en] = t2; 348 stateC = StateC.C790; 349 break; 350 351 case C790: 352 stateC = StateC.Cfinal; 353 break; 354 355 case Cfinal: 356 throw new RuntimeException("this should never happen!"); 357 358 } // end switch(stateC) 359 } // end StateLoopC 360 361 } // L790: end of LoopC: 362 363 if (VERBOSE) System.out.println("QZVEC: finished Loop C ***********************"); 364 stateA = StateA.A795; 365 break; 366 // end complex vector 367 368 case A795: 369 isw = 3 - isw; 370 stateA = StateA.A800; 371 break; 372 373 case A800: 374 stateA = StateA.Afinal; 375 break; 376 377 case Afinal: 378 throw new RuntimeException("this should never happen!"); 379 } // end switch(stateA) 380 381 } // end StateLoopA 382 } // L800: end LoopA 383 if (VERBOSE) System.out.println("QZVEC: finished Loop A ***********************"); 384 385 // end back substitution. transform to original coordinate system. 386 for (jj = 0; jj < n; jj++) { 387 j = n - jj - 1; 388 389 for (i = 0; i < n; i++) { 390 zz = 0.0; 391 for (k = 0; k <= j; k++) { 392 zz = zz + z[i][k] * b[k][j]; 393 } // L860 394 z[i][j] = zz; 395 } 396 } // L880: end for (jj ... 397 398 if (VERBOSE) System.out.println("QZVEC: starting Loop D ***********************"); 399 // normalize so that modulus of largest component of each vector is 1. 400 // (isw is 1 initially from before) 401 LoopD: for (j = 0; j < n; j++) { 402 403 StateD stateD = StateD.Dinit; 404 StateLoopD: while (stateD != StateD.Dfinal) { 405 406 switch (stateD) { 407 case Dinit : 408 d = 0.0; 409 if (isw == 2) { 410 stateD = StateD.D920; 411 break; 412 } 413 if (alfi[j] != 0.0) { 414 stateD = StateD.D945; 415 break; 416 } 417 418 for (i = 0; i < n; i++) { 419 // if ((Math.abs(z[i][j])) > d) d = (Math.abs(z[i][j])); 420 d = Math.max(d, Math.abs(z[i][j])); 421 } // L890 422 423 for (i = 0; i < n; i++) { 424 z[i][j] = z[i][j] / d; 425 } // L900 426 427 stateD = StateD.D950; 428 break; 429 430 case D920: 431 for (i = 0; i < n; i++) { 432 r = Math.abs(z[i][j - 1]) + Math.abs(z[i][j]); 433 if (r != 0.0) { 434 // r = r * dsqrt((z(i,j-1)/r)**2 + (z(i,j)/r)**2) 435 r = r * Math.hypot(z[i][j - 1] / r, z[i][j] / r); 436 } 437 if (r > d) { 438 d = r; 439 } 440 } // L930 441 442 for (i = 0; i < n; i++) { 443 z[i][j - 1] = z[i][j - 1] / d; 444 z[i][j] = z[i][j] / d; 445 } // L940 446 447 stateD = StateD.D945; 448 break; 449 450 case D945: 451 isw = 3 - isw; 452 stateD = StateD.D950; 453 break; 454 455 case D950: 456 stateD = StateD.Dfinal; 457 break; 458 459 case Dfinal: 460 throw new RuntimeException("this should never happen!"); 461 462 } // end switch(stateD) 463 } // end StateLoopD 464 465 } // L950: end LoopD 466 if (VERBOSE) System.out.println("QZVEC: finished Loop D ***********************"); 467 468 } // end of qzvec() 469 470 471 private enum StateA { 472 Ainit, Afinal, A795, A710, A800 473 } 474 475 private enum StateB { 476 Binit, Bfinal, B630, B690, B640, B700, B650 477 } 478 479 private enum StateC { 480 Cinit, Cfinal, C770, C790, C780, C773, C775, C777, C787, C782, C785 481 } 482 483 private enum StateD { 484 Dinit, Dfinal, D920, D945, D950 485 } 486 487}