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 QZVAL { 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 QZVAL() {} 020 021 /** 022 * <p> 023 * This subroutine is the third step of the qz algorithm for solving generalized matrix eigenvalue problems, Siam J. 024 * Numer. Anal. 10, 241-256 (1973) by Moler and Stewart. This description has been adapted from the original version 025 * (dated August 1983). 026 * </p> 027 * <p> 028 * This subroutine accepts a pair of real matrices, one of them in quasi-triangular form and the other in upper 029 * triangular form. it reduces the quasi-triangular matrix further, so that any remaining 2-by-2 blocks correspond 030 * to pairs of complex Eigenvalues, and returns quantities whose ratios give the generalized eigenvalues. It is 031 * usually preceded by qzhes and qzit and may be followed by qzvec. 032 * </p> 033 * <p> 034 * On output: 035 * </p> 036 * <ul> 037 * <li><strong>a</strong> has been reduced further to a quasi-triangular matrix 038 * in which all nonzero subdiagonal elements correspond to pairs of complex 039 * eigenvalues.</li> 040 * <li><strong>b</strong> is still in upper triangular form, although its 041 * elements have been altered. b[n-1][0] is unaltered.</li> 042 * <li><strong>alfr</strong> and <strong>alfi</strong> contain the real and imaginary parts of the diagonal elements 043 * of the triangular matrix that would be obtained if a were reduced completely to triangular form by unitary 044 * transformations. Non-zero values of <strong>alfi</strong> occur in pairs, the first member positive and the 045 * second negative.</li> 046 * <li><strong>beta</strong> contains the diagonal elements of the corresponding <strong>b</strong>, normalized to 047 * be real and non-negative. The generalized eigenvalues are then the ratios ((alfr + i * alfi) / beta).</li> 048 * <li><strong>z</strong> contains the product of the right hand transformations (for all three steps) if 049 * <strong>matz</strong> has been set to true.</li> 050 * </ul> 051 * 052 * @param a contains a real upper quasi-triangular matrix. 053 * @param b contains a real upper triangular matrix. In addition, location b[n-1][0] contains the tolerance quantity 054 * (epsb) computed and saved in qzit. 055 * @param alfr and alfi contain the real and imaginary parts of the diagonal elements of the triangular matrix that 056 * would be obtained if a were reduced completely to triangular form by unitary transformations. Non-zero values of 057 * <strong>alfi</strong> occur in pairs, the first member positive and the second negative. 058 * @param alfi see description for <strong>alfr</strong>. 059 * @param beta contains the diagonal elements of the corresponding 060 * <strong>b</strong>, normalized to be real and non-negative. The 061 * generalized eigenvalues are then the ratios ((alfr + i * alfi) / beta). 062 * @param matz should be set to true. if the right hand transformations are to be accumulated for later use in 063 * computing eigenvectors, and to false otherwise. 064 * @param z contains, if <strong>matz</strong> has been set to true, the transformation matrix produced in the 065 * reductions by qzhes and qzit, if performed, or else the identity matrix. If 066 * <strong>matz</strong> has been set to false, z is not referenced. 067 */ 068 public static void qzval(double[][] a, double[][] b, double[] alfr, double[] alfi, double[] beta, boolean matz, 069 double[][] z) { 070 071 final int n = a.length; 072 073 int i = 0, j; 074 int na = 0, en = 0, nn; 075 double c = 0, d = 0, e = 0; 076 double r, s, t; 077 double a1 = 0, a2 = 0, u1, u2, v1, v2; 078 double a11 = 0, a12 = 0, a21 = 0, a22 = 0; 079 double b11 = 0, b12 = 0, b22 = 0; 080 double di = 0, ei = 0; 081 double an = 0, bn = 0; 082 double cq = 0, dr = 0; 083 double cz = 0, ti = 0, tr = 0; 084 double a1i = 0, a2i = 0, a11i, a12i, a22i, a11r, a12r, a22r; 085 double sqi = 0, ssi = 0, sqr = 0, szi = 0, ssr = 0, szr = 0; 086 087 double epsb = b[n - 1][0]; 088 int isw = 1; 089 090 // find eigenvalues of quasi-triangular matrices. 091 for (nn = 0; nn < n; nn++) { 092 093 State state = State.Linit; 094 StateLoop: while (state != State.Lfinal) { 095 096 switch (state) { 097 case Linit: 098 en = n - 1 - nn; // TODO: check!! 099 na = en - 1; 100 if (isw == 2) { 101 state = State.L505; 102 break; 103 } 104 if (en == 0) { 105 state = State.L410; 106 break; 107 } 108 if (a[en][na] != 0.0) { 109 state = State.L420; 110 break; 111 } 112 state = State.L410; 113 break; 114 115 case L410: 116 // 1-by-1 block, one real root 117 alfr[en] = a[en][en]; 118 if (b[en][en] < 0.0) { 119 alfr[en] = -alfr[en]; 120 } 121 beta[en] = (Math.abs(b[en][en])); 122 alfi[en] = 0.0; 123 state = State.L510; 124 break; 125 126 case L420: 127 // 2-by-2 block 128 if (Math.abs(b[na][na]) <= epsb) { 129 state = State.L455; 130 break; 131 } 132 if (Math.abs(b[en][en]) > epsb) { 133 state = State.L430; 134 break; 135 } 136 a1 = a[en][en]; 137 a2 = a[en][na]; 138 bn = 0.0; 139 state = State.L435; 140 break; 141 142 case L430: 143 an = Math.abs(a[na][na]) + Math.abs(a[na][en]) + Math.abs(a[en][na]) + Math.abs(a[en][en]); 144 bn = Math.abs(b[na][na]) + Math.abs(b[na][en]) + Math.abs(b[en][en]); 145 a11 = a[na][na] / an; 146 a12 = a[na][en] / an; 147 a21 = a[en][na] / an; 148 a22 = a[en][en] / an; 149 b11 = b[na][na] / bn; 150 b12 = b[na][en] / bn; 151 b22 = b[en][en] / bn; 152 e = a11 / b11; 153 ei = a22 / b22; 154 s = a21 / (b11 * b22); 155 t = (a22 - e * b22) / b22; 156 157 if (Math.abs(e) > Math.abs(ei)) { 158 e = ei; 159 t = (a11 - e * b11) / b11; 160 } 161 162 c = 0.5 * (t - s * b12); 163 d = c * c + s * (a12 - e * b12); 164 165 if (d < 0.0) { 166 state = State.L480; 167 break; 168 } 169 170 // two real roots. zero both a(en,na) and b(en,na) 171 e = e + c + Math.copySign(Math.sqrt(d), c); 172 a11 -= e * b11; 173 a12 -= e * b12; 174 a22 -= e * b22; 175 176 if (Math.abs(a11) + Math.abs(a12) < Math.abs(a21) + Math.abs(a22)) { 177 a1 = a22; 178 a2 = a21; 179 } else { 180 a1 = a12; 181 a2 = a11; 182 } 183 state = State.L435; 184 break; 185 186 case L435: 187 // choose and apply real z 188 s = Math.abs(a1) + Math.abs(a2); 189 u1 = a1 / s; 190 u2 = a2 / s; 191 r = Math.copySign(Math.hypot(u1, u2), u1); 192 v1 = -(u1 + r) / r; 193 v2 = -u2 / r; 194 u2 = v2 / v1; 195 196 for (i = 0; i <= en; i++) { // TODO: check! 197 t = a[i][en] + u2 * a[i][na]; 198 a[i][en] = a[i][en] + t * v1; 199 a[i][na] = a[i][na] + t * v2; 200 201 t = b[i][en] + u2 * b[i][na]; 202 b[i][en] = b[i][en] + t * v1; 203 b[i][na] = b[i][na] + t * v2; 204 } // L440 205 206 if (matz) { 207 for (i = 0; i < n; i++) { 208 t = z[i][en] + u2 * z[i][na]; 209 z[i][en] = z[i][en] + t * v1; 210 z[i][na] = z[i][na] + t * v2; 211 } // L445 212 } 213 214 // L450 215 if (bn == 0.0) { 216 state = State.L475; 217 break; 218 } 219 if (an < Math.abs(e) * bn) { 220 state = State.L455; 221 break; 222 } 223 a1 = b[na][na]; 224 a2 = b[en][na]; 225 state = State.L460; 226 break; 227 228 case L455: 229 a1 = a[na][na]; 230 a2 = a[en][na]; 231 state = State.L475; 232 break; 233 234 case L460: 235 // choose and apply real q 236 s = Math.abs(a1) + Math.abs(a2); 237 if (s == 0.0) { 238 state = State.L475; 239 break; 240 } 241 u1 = a1 / s; 242 u2 = a2 / s; 243 r = Math.copySign(Math.hypot(u1, u2), u1); 244 v1 = -(u1 + r) / r; 245 v2 = -u2 / r; 246 u2 = v2 / v1; 247 248 for (j = na; j < n; j++) { 249 t = a[na][j] + u2 * a[en][j]; 250 a[na][j] = a[na][j] + t * v1; 251 a[en][j] = a[en][j] + t * v2; 252 253 t = b[na][j] + u2 * b[en][j]; 254 b[na][j] = b[na][j] + t * v1; 255 b[en][j] = b[en][j] + t * v2; 256 } // L470 257 258 state = State.L475; 259 break; 260 261 case L475: 262 a[en][na] = 0.0; 263 b[en][na] = 0.0; 264 alfr[na] = a[na][na]; 265 alfr[en] = a[en][en]; 266 267 if (b[na][na] < 0.0) 268 alfr[na] = -alfr[na]; 269 270 if (b[en][en] < 0.0) 271 alfr[en] = -alfr[en]; 272 273 beta[na] = (Math.abs(b[na][na])); 274 beta[en] = (Math.abs(b[en][en])); 275 alfi[en] = 0.0; 276 alfi[na] = 0.0; 277 state = State.L505; 278 break; 279 280 case L480: 281 // two complex roots 282 e = e + c; 283 ei = Math.sqrt(-d); 284 a11r = a11 - e * b11; 285 a11i = ei * b11; 286 a12r = a12 - e * b12; 287 a12i = ei * b12; 288 a22r = a22 - e * b22; 289 a22i = ei * b22; 290 291 if (Math.abs(a11r) + Math.abs(a11i) + Math.abs(a12r) + Math.abs(a12i) < 292 Math.abs(a21) + Math.abs(a22r) + Math.abs(a22i)) { 293 // L482; 294 a1 = a22r; 295 a1i = a22i; 296 a2 = -a21; 297 a2i = 0.0; 298 } else { 299 a1 = a12r; 300 a1i = a12i; 301 a2 = -a11r; 302 a2i = -a11i; 303 } 304 // goto L485; 305 state = State.L485; 306 break; 307 308 case L485: 309 // choose complex z 310 cz = Math.hypot(a1, a1i); 311 if (cz == 0.0) { 312 // L487 313 szr = 1.0; 314 szi = 0.0; 315 } else { 316 szr = (a1 * a2 + a1i * a2i) / cz; 317 szi = (a1 * a2i - a1i * a2) / cz; 318 r = Math.sqrt(cz * cz + szr * szr + szi * szi); 319 cz = cz / r; 320 szr = szr / r; 321 szi = szi / r; 322 } 323 324 // L490: 325 if (an < (Math.abs(e) + ei) * bn) { 326 // L492 327 a1 = cz * a11 + szr * a12; 328 a1i = szi * a12; 329 a2 = cz * a21 + szr * a22; 330 a2i = szi * a22; 331 } else { 332 a1 = cz * b11 + szr * b12; 333 a1i = szi * b12; 334 a2 = szr * b22; 335 a2i = szi * b22; 336 } 337 338 // L495: choose complex q 339 cq = Math.hypot(a1, a1i); 340 if (cq == 0.0) { 341 sqr = 1.0; 342 sqi = 0.0; 343 } else { 344 sqr = (a1 * a2 + a1i * a2i) / cq; 345 sqi = (a1 * a2i - a1i * a2) / cq; 346 r = Math.sqrt(cq * cq + sqr * sqr + sqi * sqi); 347 cq = cq / r; 348 sqr = sqr / r; 349 sqi = sqi / r; 350 } 351 352 // L500: compute diagonal elements that would result if transformations were applied 353 ssr = sqr * szr + sqi * szi; 354 ssi = sqr * szi - sqi * szr; 355 i = 0; 356 tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22; 357 ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22; 358 dr = cq * cz * b11 + cq * szr * b12 + ssr * b22; 359 di = cq * szi * b12 + ssi * b22; 360 state = State.L503; 361 break; 362 363 case L502: 364 i = 1; 365 tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22; 366 ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21; 367 dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22; 368 di = -ssi * b11 - sqi * cz * b12; 369 state = State.L503; 370 break; 371 372 case L503: 373 t = ti * dr - tr * di; 374 j = na; 375 if (t < 0.0) { 376 j = en; 377 } 378 r = Math.hypot(dr, di); 379 beta[j] = bn * r; 380 alfr[j] = an * (tr * dr + ti * di) / r; 381 alfi[j] = an * t / r; 382 if (i == 0) { 383 // goto L502; 384 state = State.L502; 385 break; 386 } 387 state = State.L505; 388 break; 389 390 case L505: 391 isw = 3 - isw; 392 state = State.L510; 393 break; 394 395 case L510: 396 state = State.Lfinal; 397 break; 398 399 case Lfinal: 400 throw new RuntimeException("this should never happen!"); 401 // break StateLoop; 402 } 403 } // end StateLoop: while (state 404 405 } // L510 end for (nn .. 406 407 b[n - 1][0] = epsb; 408 409 } // end of qzval() 410 411 private enum State { 412 Linit, Lfinal, L505, L410, L420, L455, L430, L435, L480, L475, L485, L503, L502, L460, L510 413 } 414 415}