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 QZHES { 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 QZHES() {} 020 021 /** 022 * <p> 023 * This subroutine is the first 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 general matrices and reduces one of them to upper Hessenberg form and the 029 * other to upper triangular form using orthogonal transformations. It is usually followed by qzit, qzval and, 030 * possibly, qzvec. 031 * </p> 032 * <p> 033 * On output: 034 * </p> 035 * <ul> 036 * <li><strong>a</strong> has been reduced to upper hessenberg form. The elements below the first subdiagonal have 037 * been set to zero.</li> 038 * <li><strong>b</strong> has been reduced to upper triangular form. The elements below the main diagonal have been 039 * set to zero.</li> 040 * <li><strong>z</strong> contains the product of the right hand transformations if matz has been set to true. 041 * Otherwise, <strong>z</strong> is not referenced.</li> 042 * </ul> 043 * 044 * @param a contains a real general matrix. 045 * @param b contains a real general matrix. 046 * @param matz should be set to true if the right hand transformations are to be accumulated for later use in 047 * computing eigenvectors, and to false otherwise. 048 * @param z on output, contains the product of the right hand transformations if matz has been set to true. 049 * Otherwise, z is not referenced. 050 */ 051 public static void qzhes(double[][] a, double[][] b, boolean matz, double[][] z) { 052 053 final int n = a.length; 054 055 int i, j, k, l, nm1, nm2; 056 double r, s, t; 057 int l1; 058 double u1, u2, v1, v2; 059 int lb, nk1; 060 double rho; 061 062 // initialize z (with identity matrix) 063 if (matz) { 064 for (j = 0; j < n; j++) { 065 for (i = 0; i < n; i++) { 066 z[i][j] = 0.0; 067 } // L2 068 z[j][j] = 1.0; 069 } // L3 070 } 071 072 // L10: reduce b to upper triangular form 073 if (n <= 1) { 074 return; 075 } 076 077 nm1 = n - 1; // TODO: check!! 078 079 for (l = 0; l < nm1; l++) { 080 l1 = l + 1; 081 s = 0.0; 082 083 for (i = l1; i < n; i++) { 084 s = s + Math.abs(b[i][l]); 085 } // L20 086 087 if (s != 0.0) { 088 s = s + Math.abs(b[l][l]); 089 r = 0.0; 090 for (i = l; i < n; i++) { 091 b[i][l] = b[i][l] / s; 092 r = r + b[i][l] * b[i][l]; 093 } // L25 094 095 r = Math.copySign(Math.sqrt(r), b[l][l]); 096 b[l][l] = b[l][l] + r; 097 rho = r * b[l][l]; 098 099 for (j = l1; j < n; j++) { 100 t = 0.0; // dot product 101 for (i = l; i < n; i++) { 102 t = t + b[i][l] * b[i][j]; 103 } //L30 104 t = -t / rho; 105 106 for (i = l; i < n; i++) { 107 b[i][j] = b[i][j] + t * b[i][l]; 108 } // L40 109 110 } // L50 111 112 for (j = 0; j < n; j++) { 113 t = 0.0; // dot product 114 for (i = l; i < n; i++) { 115 t = t + b[i][l] * a[i][j]; 116 } // L60 117 118 t = -t / rho; 119 120 for (i = l; i < n; i++) { 121 a[i][j] = a[i][j] + t * b[i][l]; 122 } // L70 123 124 } // L80 125 126 b[l][l] = -s * r; 127 for (i = l1; i < n; i++) { 128 b[i][l] = 0.0; 129 } // L90 130 131 } // end if 132 } // L100 133 134 // reduce a to upper Hessenberg form, while keeping b triangular 135 if (n == 2) { 136 return; 137 } 138 nm2 = n - 2; 139 140 for (k = 0; k < nm2; k++) { 141 nk1 = nm1 - k - 1; // = n - 2 - k; 142 143 // for l=n-1 step -1 until k+1 do 144 for (lb = 0; lb < nk1; lb++) { 145 l = n - lb - 2; // TODO: check!! 146 l1 = l + 1; 147 148 // zero a(l+1,k) 149 s = Math.abs(a[l][k]) + Math.abs(a[l1][k]); 150 if (s == 0.0) { 151 continue; // goto 150 152 } 153 u1 = a[l][k] / s; 154 u2 = a[l1][k] / s; 155 r = Math.copySign(Math.hypot(u1, u2), u1); 156 v1 = -(u1 + r) / r; 157 v2 = -u2 / r; 158 u2 = v2 / v1; 159 160 for (j = k; j < n; j++) { 161 t = a[l][j] + u2 * a[l1][j]; 162 a[l][j] = a[l][j] + t * v1; 163 a[l1][j] = a[l1][j] + t * v2; 164 } // L110 165 166 a[l1][k] = 0.0; 167 168 for (j = l; j < n; j++) { 169 t = b[l][j] + u2 * b[l1][j]; 170 b[l][j] = b[l][j] + t * v1; 171 b[l1][j] = b[l1][j] + t * v2; 172 } // L120 173 174 // zero b(l+1,l) 175 s = Math.abs(b[l1][l1]) + Math.abs(b[l1][l]); 176 if (s == 0.0) { 177 continue; // go to 150 178 } 179 u1 = b[l1][l1] / s; 180 u2 = b[l1][l] / s; 181 r = Math.copySign(Math.hypot(u1, u2), u1); 182 v1 = -(u1 + r) / r; 183 v2 = -u2 / r; 184 u2 = v2 / v1; 185 186 for (i = 0; i <= l1; i++) { 187 t = b[i][l1] + u2 * b[i][l]; 188 b[i][l1] = b[i][l1] + t * v1; 189 b[i][l] = b[i][l] + t * v2; 190 } // L130 191 192 b[l1][l] = 0.0; 193 194 for (i = 0; i < n; i++) { 195 t = a[i][l1] + u2 * a[i][l]; 196 a[i][l1] = a[i][l1] + t * v1; 197 a[i][l] = a[i][l] + t * v2; 198 } // L140 199 200 if (matz) { 201 for (i = 0; i < n; i++) { 202 t = z[i][l1] + u2 * z[i][l]; 203 z[i][l1] = z[i][l1] + t * v1; 204 z[i][l] = z[i][l] + t * v2; 205 } // L145 206 } 207 } // L150 208 } // L160 209 210 } // end of qzhes 211 212}