subroutine bakvec ( n, t, e, m, z, ierr ) !*****************************************************************************80 ! !! BAKVEC determines eigenvectors by reversing the FIGI transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a nonsymmetric tridiagonal ! matrix by back transforming those of the corresponding symmetric ! matrix determined by FIGI. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, real ( kind = 8 ) T(N,3), contains the nonsymmetric matrix. Its ! subdiagonal is stored in the positions 2:N of the first column, ! its diagonal in positions 1:N of the second column, ! and its superdiagonal in positions 1:N-1 of the third column. ! T(1,1) and T(N,3) are arbitrary. ! ! Input/output, real ( kind = 8 ) E(N). On input, E(2:N) contains the ! subdiagonal elements of the symmetric matrix. E(1) is arbitrary. ! On output, the contents of E have been destroyed. ! ! Input, integer ( kind = 4 ) M, the number of eigenvectors to be back ! transformed. ! ! Input/output, real ( kind = 8 ) Z(N,M), contains the eigenvectors. ! On output, they have been transformed as requested. ! ! Output, integer ( kind = 4 ) IERR, an error flag. ! 0, for normal return, ! 2*N+I, if E(I) is zero with T(I,1) or T(I-1,3) non-zero. ! In this case, the symmetric matrix is not similar ! to the original matrix, and the eigenvectors ! cannot be found by this program. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) e(n) integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) j real ( kind = 8 ) t(n,3) real ( kind = 8 ) z(n,m) ierr = 0 if ( m == 0 ) then return end if e(1) = 1.0D+00 if ( n == 1 ) then return end if do i = 2, n if ( e(i) == 0.0D+00 ) then if ( t(i,1) /= 0.0D+00 .or. t(i-1,3) /= 0.0D+00 ) then ierr = 2 * n + i return end if e(i) = 1.0D+00 else e(i) = e(i-1) * e(i) / t(i-1,3) end if end do do j = 1, m z(2:n,j) = z(2:n,j) * e(2:n) end do return end subroutine balanc ( n, a, low, igh, scale ) !*****************************************************************************80 ! !! BALANC balances a real matrix before eigenvalue calculations. ! ! Discussion: ! ! This subroutine balances a real matrix and isolates eigenvalues ! whenever possible. ! ! Suppose that the principal submatrix in rows LOW through IGH ! has been balanced, that P(J) denotes the index interchanged ! with J during the permutation step, and that the elements ! of the diagonal matrix used are denoted by D(I,J). Then ! ! SCALE(J) = P(J), J = 1,...,LOW-1, ! = D(J,J), J = LOW,...,IGH, ! = P(J) J = IGH+1,...,N. ! ! The order in which the interchanges are made is N to IGH+1, ! then 1 to LOW-1. ! ! Note that 1 is returned for LOW if IGH is zero formally. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input/output, real ( kind = 8 ) A(N,N), the N by N matrix. On output, ! the matrix has been balanced. ! ! Output, integer ( kind = 4 ) LOW, IGH, indicate that A(I,J) is equal to zero if ! (1) I is greater than J and ! (2) J=1,...,LOW-1 or I=IGH+1,...,N. ! ! Output, real ( kind = 8 ) SCALE(N), contains information determining the ! permutations and scaling factors used. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n,n) real ( kind = 8 ) b2 real ( kind = 8 ) c real ( kind = 8 ) f real ( kind = 8 ) g integer ( kind = 4 ) i integer ( kind = 4 ) iexc integer ( kind = 4 ) igh integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) low integer ( kind = 4 ) m logical noconv real ( kind = 8 ) r real ( kind = 8 ) radix real ( kind = 8 ) s real ( kind = 8 ) scale(n) radix = 16.0D+00 iexc = 0 j = 0 m = 0 b2 = radix**2 k = 1 l = n go to 100 20 continue scale(m) = j if ( j /= m ) then do i = 1, l call r8_swap ( a(i,j), a(i,m) ) end do do i = k, n call r8_swap ( a(j,i), a(m,i) ) end do end if 50 continue if ( iexc == 2 ) then go to 130 end if ! ! Search for rows isolating an eigenvalue and push them down. ! 80 continue if ( l == 1 ) then low = k igh = l return end if l = l - 1 100 continue do j = l, 1, -1 do i = 1, l if ( i /= j ) then if ( a(j,i) /= 0.0D+00 ) then go to 120 end if end if end do m = l iexc = 1 go to 20 120 continue end do go to 140 ! ! Search for columns isolating an eigenvalue and push them left. ! 130 continue k = k + 1 140 continue do j = k, l do i = k, l if ( i /= j ) then if ( a(i,j) /= 0.0D+00 ) then go to 170 end if end if end do m = k iexc = 2 go to 20 170 continue end do ! ! Balance the submatrix in rows K to L. ! scale(k:l) = 1.0D+00 ! ! Iterative loop for norm reduction. ! noconv = .true. do while ( noconv ) noconv = .false. do i = k, l c = 0.0D+00 r = 0.0D+00 do j = k, l if ( j /= i ) then c = c + abs ( a(j,i) ) r = r + abs ( a(i,j) ) end if end do ! ! Guard against zero C or R due to underflow. ! if ( c /= 0.0D+00 .and. r /= 0.0D+00 ) then g = r / radix f = 1.0D+00 s = c + r do while ( c < g ) f = f * radix c = c * b2 end do g = r * radix do while ( g <= c ) f = f / radix c = c / b2 end do ! ! Balance. ! if ( ( c + r ) / f < 0.95D+00 * s ) then g = 1.0D+00 / f scale(i) = scale(i) * f noconv = .true. a(i,k:n) = a(i,k:n) * g a(1:l,i) = a(1:l,i) * f end if end if end do end do low = k igh = l return end subroutine balbak ( n, low, igh, scale, m, z ) !*****************************************************************************80 ! !! BALBAK determines eigenvectors by undoing the BALANC transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a real general matrix by ! back transforming those of the corresponding balanced matrix ! determined by BALANC. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! Parlett and Reinsch, ! Numerische Mathematik, ! Volume 13, pages 293-304, 1969. ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, column indices determined by BALANC. ! ! Input, real ( kind = 8 ) SCALE(N), contains information determining ! the permutations and scaling factors used by BALANC. ! ! Input, integer ( kind = 4 ) M, the number of columns of Z to be back-transformed. ! ! Input/output, real ( kind = 8 ) Z(N,M), contains the real and imaginary parts ! of the eigenvectors, which, on return, have been back-transformed. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) igh integer ( kind = 4 ) ii integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) low real ( kind = 8 ) s real ( kind = 8 ) scale(n) real ( kind = 8 ) z(n,m) if ( m <= 0 ) then return end if if ( igh /= low ) then do i = low, igh z(i,1:m) = z(i,1:m) * scale(i) end do end if do ii = 1, n i = ii if ( i < low .or. igh < i ) then if ( i < low ) then i = low - ii end if k = int ( scale(i) ) if ( k /= i ) then do j = 1, m call r8_swap ( z(i,j), z(k,j) ) end do end if end if end do return end subroutine bandr ( n, mb, a, d, e, e2, matz, z ) !*****************************************************************************80 ! !! BANDR reduces a symmetric band matrix to symmetric tridiagonal form. ! ! Discussion: ! ! This subroutine reduces a real symmetric band matrix ! to a symmetric tridiagonal matrix using and optionally ! accumulating orthogonal similarity transformations. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) MB, is the (half) band width of the matrix, defined as the ! number of adjacent diagonals, including the principal ! diagonal, required to specify the non-zero portion of the ! lower triangle of the matrix. ! ! Input/output, real ( kind = 8 ) A(N,MB). On input, contains the lower triangle of ! the symmetric band input matrix stored as an N by MB array. Its lowest ! subdiagonal is stored in the last N+1-MB positions of the first column, ! its next subdiagonal in the last N+2-MB positions of the second column, ! further subdiagonals similarly, and finally its principal diagonal in ! the N positions of the last column. Contents of storages not part of ! the matrix are arbitrary. On output, A has been destroyed, except for ! its last two columns which contain a copy of the tridiagonal matrix. ! ! Output, real ( kind = 8 ) D(N), the diagonal elements of the tridiagonal matrix. ! ! Output, real ( kind = 8 ) E(N), the subdiagonal elements of the tridiagonal ! matrix in E(2:N). E(1) is set to zero. ! ! Output, real ( kind = 8 ) E2(N), contains the squares of the corresponding elements ! of E. E2 may coincide with E if the squares are not needed. ! ! Input, logical MATZ, should be set to TRUE if the transformation matrix is ! to be accumulated, and to FALSE otherwise. ! ! Output, real ( kind = 8 ) Z(N,N), the orthogonal transformation matrix produced in ! the reduction if MATZ has been set to TRUE. Otherwise, Z is not ! referenced. ! implicit none integer ( kind = 4 ) mb integer ( kind = 4 ) n real ( kind = 8 ) a(n,mb) real ( kind = 8 ) b1 real ( kind = 8 ) b2 real ( kind = 8 ) c2 real ( kind = 8 ) d(n) real ( kind = 8 ) dmin real ( kind = 8 ) dminrt real ( kind = 8 ) e(n) real ( kind = 8 ) e2(n) real ( kind = 8 ) f1 real ( kind = 8 ) f2 real ( kind = 8 ) g integer ( kind = 4 ) i1 integer ( kind = 4 ) i2 integer ( kind = 4 ) j integer ( kind = 4 ) j1 integer ( kind = 4 ) j2 integer ( kind = 4 ) k integer ( kind = 4 ) kr integer ( kind = 4 ) l integer ( kind = 4 ) m1 logical matz integer ( kind = 4 ) maxl integer ( kind = 4 ) maxr integer ( kind = 4 ) mr integer ( kind = 4 ) r integer ( kind = 4 ) r1 real ( kind = 8 ) s2 real ( kind = 8 ) u integer ( kind = 4 ) ugl real ( kind = 8 ) z(n,n) dmin = epsilon ( dmin ) dminrt = sqrt ( dmin ) ! ! Initialize the diagonal scaling matrix. ! d(1:n) = 1.0D+00 if ( matz ) then call r8mat_identity ( n, z ) end if m1 = mb - 1 if ( m1 < 1 ) then d(1:n) = a(1:n,mb) e(1:n) = 0.0D+00 e2(1:n) = 0.0D+00 return end if if ( m1 == 1 ) then go to 800 end if do k = 1, n - 2 maxr = min ( m1, n-k ) do r1 = 2, maxr r = maxr + 2 - r1 kr = k + r mr = mb - r g = a(kr,mr) a(kr-1,1) = a(kr-1,mr+1) ugl = k do j = kr, n, m1 j1 = j - 1 j2 = j1 - 1 if ( g == 0.0D+00 ) then go to 600 end if b1 = a(j1,1) / g b2 = b1 * d(j1) / d(j) s2 = 1.0D+00 / ( 1.0D+00 + b1 * b2 ) if ( s2 < 0.5D+00 ) then b1 = g / a(j1,1) b2 = b1 * d(j) / d(j1) c2 = 1.0D+00 - s2 d(j1) = c2 * d(j1) d(j) = c2 * d(j) f1 = 2.0D+00 * a(j,m1) f2 = b1 * a(j1,mb) a(j,m1) = -b2 * ( b1 * a(j,m1) - a(j,mb) ) - f2 + a(j,m1) a(j1,mb) = b2 * ( b2 * a(j,mb) + f1 ) + a(j1,mb) a(j,mb) = b1 * ( f2 - f1 ) + a(j,mb) do l = ugl, j2 i2 = mb - j + l u = a(j1,i2+1) + b2 * a(j,i2) a(j,i2) = -b1 * a(j1,i2+1) + a(j,i2) a(j1,i2+1) = u end do ugl = j a(j1,1) = a(j1,1) + b2 * g if ( j /= n ) then maxl = min ( m1, n-j1 ) do l = 2, maxl i1 = j1 + l i2 = mb - l u = a(i1,i2) + b2 * a(i1,i2+1) a(i1,i2+1) = -b1 * a(i1,i2) + a(i1,i2+1) a(i1,i2) = u end do i1 = j + m1 if ( i1 <= n ) then g = b2 * a(i1,1) end if end if if ( matz ) then do l = 1, n u = z(l,j1) + b2 * z(l,j) z(l,j) = -b1 * z(l,j1) + z(l,j) z(l,j1) = u end do end if else u = d(j1) d(j1) = s2 * d(j) d(j) = s2 * u f1 = 2.0D+00 * a(j,m1) f2 = b1 * a(j,mb) u = b1 * ( f2 - f1 ) + a(j1,mb) a(j,m1) = b2 * ( b1 * a(j,m1) - a(j1,mb) ) + f2 - a(j,m1) a(j1,mb) = b2 * ( b2 * a(j1,mb) + f1 ) + a(j,mb) a(j,mb) = u do l = ugl, j2 i2 = mb - j + l u = b2 * a(j1,i2+1) + a(j,i2) a(j,i2) = -a(j1,i2+1) + b1 * a(j,i2) a(j1,i2+1) = u end do ugl = j a(j1,1) = b2 * a(j1,1) + g if ( j /= n ) then maxl = min ( m1, n-j1 ) do l = 2, maxl i1 = j1 + l i2 = mb - l u = b2 * a(i1,i2) + a(i1,i2+1) a(i1,i2+1) = -a(i1,i2) + b1 * a(i1,i2+1) a(i1,i2) = u end do i1 = j + m1 if ( i1 <= n ) then g = a(i1,1) a(i1,1) = b1 * a(i1,1) end if end if if ( matz ) then do l = 1, n u = b2 * z(l,j1) + z(l,j) z(l,j) = -z(l,j1) + b1 * z(l,j) z(l,j1) = u end do end if end if end do 600 continue end do ! ! Rescale to avoid underflow or overflow. ! if ( mod ( k, 64 ) == 0 ) then do j = k, n if ( d(j) < dmin ) then maxl = max ( 1, mb+1-j ) a(j,maxl:m1) = dminrt * a(j,maxl:m1) if ( j /= n ) then maxl = min ( m1, n-j ) do l = 1, maxl i1 = j + l i2 = mb - l a(i1,i2) = dminrt * a(i1,i2) end do end if if ( matz ) then z(1:n,j) = dminrt * z(1:n,j) end if a(j,mb) = dmin * a(j,mb) d(j) = d(j) / dmin end if end do end if end do ! ! Form square root of scaling matrix. ! 800 continue e(2:n) = sqrt ( d(2:n) ) if ( matz ) then do k = 2, n z(1:n,k) = z(1:n,k) * e(k) end do end if u = 1.0D+00 do j = 2, n a(j,m1) = u * e(j) * a(j,m1) u = e(j) e2(j) = a(j,m1)**2 a(j,mb) = d(j) * a(j,mb) d(j) = a(j,mb) e(j) = a(j,m1) end do d(1) = a(1,mb) e(1) = 0.0D+00 e2(1) = 0.0D+00 return end subroutine bandv ( n, mbw, a, e21, m, w, z, ierr ) !*****************************************************************************80 ! !! BANDV finds eigenvectors from eigenvalues, for a real symmetric band matrix. ! ! Discussion: ! ! This subroutine finds those eigenvectors of a real symmetric ! band matrix corresponding to specified eigenvalues, using inverse ! iteration. The subroutine may also be used to solve systems ! of linear equations with a symmetric or non-symmetric band ! coefficient matrix. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) MBW, the number of columns of the array A used to store the ! band matrix. If the matrix is symmetric, MBW is its (half) ! band width, denoted MB and defined as the number of adjacent ! diagonals, including the principal diagonal, required to ! specify the non-zero portion of the lower triangle of the ! matrix. If the subroutine is being used to solve systems ! of linear equations and the coefficient matrix is not ! symmetric, it must however have the same number of adjacent ! diagonals above the main diagonal as below, and in this ! case, MBW=2*MB-1. ! ! Input, real ( kind = 8 ) A(N,MBW), the lower triangle of the symmetric band input ! atrix stored as an N by MB array. Its lowest subdiagonal is stored ! in the last N+1-MB positions of the first column, its next subdiagonal ! in the last N+2-MB positions of the second column, further subdiagonals ! similarly, and finally its principal diagonal in the N positions of ! column MB. If the subroutine is being used to solve systems of linear ! equations, and the coefficient matrix is not symmetric, A is ! N by 2*MB-1 instead, with lower triangle as above and with its first ! superdiagonal stored in the first N-1 positions of column MB+1, its ! second superdiagonal in the first N-2 positions of column MB+2, further ! superdiagonals similarly, and finally its highest superdiagonal in ! the first N+1-MB positions of the last column. Contents of storages ! not part of the matrix are arbitrary. ! ! Input, real ( kind = 8 ) E21, specifies the ordering of the eigenvalues and contains ! 0.0D+00 if the eigenvalues are in ascending order, or 2.0D+00 if the ! eigenvalues are in descending order. If the subroutine is being used ! to solve systems of linear equations, E21 should be set to 1.0D+00 ! if the coefficient matrix is symmetric and to -1.0D+00 if not. ! ! Input, integer ( kind = 4 ) M, the number of specified eigenvalues or the number of ! systems of linear equations. ! ! Input, real ( kind = 8 ) W(M), contains the M eigenvalues in ascending or descending ! order. If the subroutine is being used to solve systems of linear ! equations (A-W(1:M)*I) * X(1:M) = B(1:M), where I is the identity matrix, ! W should be set accordingly. ! ! Input/output, real ( kind = 8 ) Z(N,M). On input, the constant matrix ! columns B(1:M), if the subroutine is used to solve systems of linear ! equations. On output, the associated set of orthogonal eigenvectors. ! Any vector which fails to converge is set to zero. If the ! routine is used to solve systems of linear equations, ! Z contains the solution matrix columns X(1:M). ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! -R, if the eigenvector corresponding to the R-th eigenvalue fails to ! converge, or if the R-th system of linear equations is nearly singular. ! implicit none integer ( kind = 4 ) mbw integer ( kind = 4 ) n real ( kind = 8 ) a(n,mbw) real ( kind = 8 ) e21 real ( kind = 8 ) eps2 real ( kind = 8 ) eps3 real ( kind = 8 ) eps4 integer ( kind = 4 ) group integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) ii integer ( kind = 4 ) ij integer ( kind = 4 ) ij1 integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) k integer ( kind = 4 ) kj integer ( kind = 4 ) kj1 integer ( kind = 4 ) m integer ( kind = 4 ) m1 integer ( kind = 4 ) m21 integer ( kind = 4 ) maxj integer ( kind = 4 ) maxk integer ( kind = 4 ) mb real ( kind = 8 ) norm real ( kind = 8 ) order real ( kind = 8 ) pythag integer ( kind = 4 ) r real ( kind = 8 ) rv(n*(2*mbw-1)) real ( kind = 8 ) rv6(n) real ( kind = 8 ) u real ( kind = 8 ) uk real ( kind = 8 ) v real ( kind = 8 ) w(m) real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) xu real ( kind = 8 ) z(n,m) ierr = 0 if ( m == 0 ) then return end if x0 = 0.0D+00 if ( e21 < 0.0D+00 ) then mb = ( mbw + 1 ) / 2 else mb = mbw end if m1 = mb - 1 m21 = m1 + mb order = 1.0D+00 - abs ( e21 ) ! ! Find vectors by inverse iteration. ! do r = 1, m its = 1 x1 = w(r) if ( r /= 1 ) go to 100 ! ! Compute norm of matrix. ! norm = 0.0D+00 do j = 1, mb jj = mb + 1 - j kj = jj + m1 ij = 1 v = 0.0D+00 do i = jj, n v = v + abs ( a(i,j) ) if ( e21 < 0.0D+00 ) then v = v + abs ( a(ij,kj) ) ij = ij + 1 end if end do norm = max ( norm, v ) end do if ( e21 < 0.0D+00 ) then norm = 0.5D+00 * norm end if ! ! EPS2 is the criterion for grouping, ! EPS3 replaces zero pivots and equal roots are modified by eps3, ! EPS4 is taken very small to avoid overflow. ! if ( norm == 0.0D+00 ) then norm = 1.0D+00 end if eps2 = 0.001D+00 * norm * abs ( order) eps3 = abs ( norm ) * epsilon ( norm ) uk = n uk = sqrt ( uk ) eps4 = uk * eps3 80 continue group = 0 go to 120 ! ! Look for close or coincident roots. ! 100 continue if ( eps2 <= abs ( x1 - x0 ) ) then go to 80 end if group = group + 1 if ( order * ( x1 - x0 ) <= 0.0D+00 ) then x1 = x0 + order * eps3 end if ! ! Expand matrix, subtract eigenvalue, and initialize vector. ! 120 continue do i = 1, n ij = i + min ( 0, i-m1 ) * n kj = ij + mb * n ij1 = kj + m1 * n if ( m1 == 0 ) go to 180 do j = 1, m1 if ( ij <= m1 ) then if ( ij <= 0 ) then rv(ij1) = 0.0D+00 ij1 = ij1 + n end if else rv(ij) = a(i,j) end if ij = ij + n ii = i + j if ( ii <= n ) then jj = mb - j if ( e21 < 0.0D+00 ) then ii = i jj = mb + j end if rv(kj) = a(ii,jj) kj = kj + n end if end do 180 continue rv(ij) = a(i,mb) - x1 rv6(i) = eps4 if ( order == 0.0D+00 ) then rv6(i) = z(i,r) end if end do if ( m1 /= 0 ) then ! ! Elimination with interchanges. ! do i = 1, n ii = i + 1 maxk = min ( i+m1-1, n ) maxj = min ( n-i, m21-2 ) * n do k = i, maxk kj1 = k j = kj1 + n jj = j + maxj do kj = j, jj, n rv(kj1) = rv(kj) kj1 = kj end do rv(kj1) = 0.0D+00 end do if ( i /= n ) then u = 0.0D+00 maxk = min ( i+m1, n ) maxj = min ( n-ii, m21-2 ) * n do j = i, maxk if ( abs ( u ) <= abs ( rv(j) ) ) then u = rv(j) k = j end if end do j = i + n jj = j + maxj if ( k /= i ) then kj = k do ij = i, jj, n call r8_swap ( rv(ij), rv(kj) ) kj = kj + n end do if ( order == 0.0D+00 ) then call r8_swap ( rv6(i), rv6(k) ) end if end if if ( u /= 0.0D+00 ) then do k = ii, maxk v = rv(k) / u kj = k do ij = j, jj, n kj = kj + n rv(kj) = rv(kj) - v * rv(ij) end do if ( order == 0.0D+00 ) then rv6(k) = rv6(k) - v * rv6(i) end if end do end if end if end do end if ! ! Back substitution. ! 600 continue do ii = 1, n i = n + 1 - ii maxj = min ( ii, m21 ) if ( maxj /= 1 ) then ij1 = i j = ij1 + n jj = j + (maxj - 2) * n do ij = j, jj, n ij1 = ij1 + 1 rv6(i) = rv6(i) - rv(ij) * rv6(ij1) end do end if v = rv(i) ! ! Error: nearly singular linear system. ! if ( abs ( v ) < eps3) then if ( order == 0.0D+00 ) then ierr = -r end if v = sign ( eps3, v ) end if rv6(i) = rv6(i) / v end do xu = 1.0D+00 if ( order == 0.0D+00 ) go to 870 ! ! Orthogonalize with respect to previous members of group. ! do jj = 1, group j = r - group - 1 + jj xu = dot_product ( rv6(1:n), z(1:n,j) ) rv6(1:n) = rv6(1:n) - xu * z(1:n,j) end do norm = sum ( abs ( rv6(1:n) ) ) ! ! Choose a new starting vector. ! if ( norm < 0.1D+00 ) then if ( its < n ) then its = its + 1 xu = eps4 / ( uk + 1.0D+00 ) rv6(1) = eps4 rv6(2:n) = xu rv6(its) = rv6(its) - eps4 * uk go to 600 else ierr = -r xu = 0.0D+00 go to 870 end if end if ! ! Normalize so that sum of squares is 1 and expand to full order. ! u = 0.0D+00 do i = 1, n u = pythag ( u, rv6(i) ) end do xu = 1.0D+00 / u 870 continue z(1:n,r) = rv6(1:n) * xu x0 = x1 end do return end subroutine bisect ( n, eps1, d, e, e2, lb, ub, mm, m, w, ind, ierr ) !*****************************************************************************80 ! !! BISECT computes some eigenvalues of a real symmetric tridiagonal matrix. ! ! Discussion: ! ! This subroutine finds those eigenvalues of a real symmetric ! tridiagonal matrix which lie in a specified interval, using bisection. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input/output, real ( kind = 8 ) EPS1, is an absolute error tolerance for the computed ! eigenvalues. If the input EPS1 is non-positive, it is reset for each ! submatrix to a default value, namely, minus the product of the relative ! machine precision and the 1-norm of the submatrix. ! ! Input, real ( kind = 8 ) D(N), the diagonal elements of the input matrix. ! ! Input, real ( kind = 8 ) E(N), contains in E(2:N) the subdiagonal elements of the ! matrix. E(1) is arbitrary. ! ! Input/output, real ( kind = 8 ) E2(N). On input, the squares of the corresponding ! elements of E. E2(1) is arbitrary. On output, elements of E2, ! corresponding to elements of E regarded as negligible, have been ! replaced by zero, causing the matrix to split into a direct sum of ! submatrices. E2(1) is also set to zero. ! ! Input, real ( kind = 8 ) LB, UB, define the interval to be searched for eigenvalues. ! If LB is not less than UB, no eigenvalues will be found. ! ! Input, integer ( kind = 4 ) MM, an upper bound for the number of eigenvalues in the ! interval. Warning: if more than MM eigenvalues are determined to lie ! in the interval, an error return is made with no eigenvalues found. ! ! Output, integer ( kind = 4 ) M, the number of eigenvalues determined to lie ! in (LB,UB). ! ! Output, real ( kind = 8 ) W(M), the eigenvalues in ascending order. ! ! Output, integer ( kind = 4 ) IND(MM), contains in its first M positions the submatrix ! indices associated with the corresponding eigenvalues in W: ! 1 for eigenvalues belonging to the first submatrix from the top, 2 for ! those belonging to the second submatrix, and so on. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! 3*N+1, if M exceeds MM. ! implicit none integer ( kind = 4 ) mm integer ( kind = 4 ) n real ( kind = 8 ) d(n) real ( kind = 8 ) e(n) real ( kind = 8 ) e2(n) real ( kind = 8 ) eps1 integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) ii integer ( kind = 4 ) ind(mm) integer ( kind = 4 ) isturm integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l real ( kind = 8 ) lb integer ( kind = 4 ) m integer ( kind = 4 ) m1 integer ( kind = 4 ) m2 integer ( kind = 4 ) p integer ( kind = 4 ) q integer ( kind = 4 ) r real ( kind = 8 ) rv4(n) real ( kind = 8 ) rv5(n) integer ( kind = 4 ) s real ( kind = 8 ) t1 real ( kind = 8 ) t2 integer ( kind = 4 ) tag real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) u real ( kind = 8 ) ub real ( kind = 8 ) v real ( kind = 8 ) w(mm) real ( kind = 8 ) x0 real ( kind = 8 ) x1 real ( kind = 8 ) xu ierr = 0 s = 0 tag = 0 t1 = lb t2 = ub ! ! Look for small sub-diagonal entries. ! e2(1) = 0.0D+00 do i = 2, n tst1 = abs ( d(i) ) + abs ( d(i-1) ) tst2 = tst1 + abs ( e(i) ) if ( tst2 <= tst1 ) then e2(i) = 0.0D+00 end if end do ! ! Determine the number of eigenvalues in the interval. ! p = 1 q = n x1 = ub isturm = 1 go to 320 60 continue m = s x1 = lb isturm = 2 go to 320 80 continue m = m - s if ( mm < m ) then go to 980 end if q = 0 r = 0 ! ! Establish and process next submatrix, refining ! interval by the Gerschgorin bounds. ! 100 continue if ( r == m ) go to 1001 tag = tag + 1 p = q + 1 xu = d(p) x0 = d(p) u = 0.0D+00 do q = p, n x1 = u u = 0.0D+00 v = 0.0D+00 if ( q /= n ) then u = abs ( e(q+1) ) v = e2(q+1) end if xu = min ( d(q) - ( x1 + u ), xu ) x0 = max ( d(q) + ( x1 + u ), x0 ) if ( v == 0.0D+00 ) then exit end if end do x1 = max ( abs ( xu ), abs ( x0 ) ) * epsilon ( x1 ) if ( eps1 <= 0.0D+00 ) then eps1 = -x1 end if if ( p /= q ) go to 180 ! ! Check for an isolated root within interval. ! if ( d(p) < t1 .or. t2 <= d(p) ) then go to 940 end if m1 = p m2 = p rv5(p) = d(p) go to 900 180 continue x1 = x1 * ( q - p + 1 ) lb = max ( t1, xu - x1 ) ub = min ( t2, x0 + x1 ) x1 = lb isturm = 3 go to 320 200 continue m1 = s + 1 x1 = ub isturm = 4 go to 320 220 continue m2 = s if ( m2 < m1 ) then go to 940 end if ! ! Find roots by bisection. ! x0 = ub isturm = 5 rv5(m1:m2) = ub rv4(m1:m2) = lb ! ! Loop for the K-th eigenvalue. ! k = m2 250 continue xu = lb do ii = m1, k i = m1 + k - ii if ( xu < rv4(i) ) then xu = rv4(i) go to 280 end if end do 280 continue x0 = min ( x0, rv5(k) ) ! ! Next bisection step. ! 300 continue x1 = ( xu + x0 ) * 0.5D+00 if ( (x0 - xu) <= abs ( eps1 ) ) go to 420 tst1 = 2.0D+00 * ( abs ( xu ) + abs ( x0 ) ) tst2 = tst1 + ( x0 - xu ) if ( tst2 == tst1 ) go to 420 ! ! Sturm sequence. ! 320 continue s = p - 1 u = 1.0D+00 do i = p, q if ( u == 0.0D+00 ) then v = abs ( e(i) ) / epsilon ( v ) if ( e2(i) == 0.0D+00 ) v = 0.0D+00 else v = e2(i) / u end if u = d(i) - x1 - v if ( u < 0.0D+00 ) then s = s + 1 end if end do go to (60,80,200,220,360), isturm ! ! Refine intervals. ! 360 continue if ( k <= s ) then go to 400 end if xu = x1 if ( s < m1 ) then rv4(m1) = x1 go to 300 end if 380 continue rv4(s+1) = x1 if ( x1 < rv5(s) ) then rv5(s) = x1 end if go to 300 400 continue x0 = x1 go to 300 ! ! K-th eigenvalue found. ! 420 continue rv5(k) = x1 k = k - 1 if ( k >= m1 ) go to 250 ! ! Order eigenvalues tagged with their submatrix associations. ! 900 continue s = r r = r + m2 - m1 + 1 j = 1 k = m1 do l = 1, r if ( j <= s ) then if ( k > m2 ) then exit end if if ( rv5(k) >= w(l) ) then j = j + 1 cycle end if do ii = j, s i = l + s - ii w(i+1) = w(i) ind(i+1) = ind(i) end do end if w(l) = rv5(k) ind(l) = tag k = k + 1 end do 940 continue if ( q < n ) then go to 100 end if go to 1001 ! ! Set error: underestimate of number of eigenvalues in interval. ! 980 continue ierr = 3 * n + 1 1001 continue lb = t1 ub = t2 return end subroutine bqr ( n, mb, a, t, r, ierr ) !*****************************************************************************80 ! !! BQR finds the smallest eigenvalue of a real symmetric band matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalue of smallest magnitude of a real ! symmetric band matrix using the QR algorithm with shifts of origin. ! Consecutive calls can be made to find further eigenvalues. ! ! Note that for a subsequent call, N should be replaced by N-1, but ! MB should not be altered even when it exceeds the current N. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) MB, the (half) band width of the matrix, defined as the ! number of adjacent diagonals, including the principal ! diagonal, required to specify the non-zero portion of the ! lower triangle of the matrix. ! ! Input/output, real ( kind = 8 ) A(N,MB). On input, A contains the lower triangle ! of the symmetric band input matrix stored as an N by MB array. Its ! lowest subdiagonal is stored in the last N+1-MB positions of the first ! column, its next subdiagonal in the last N+2-MB positions of the ! second column, further subdiagonals similarly, and finally its principal ! diagonal in the N positions of the last column. Contents of storages ! not part of the matrix are arbitrary. On a subsequent call, its output ! contents from the previous call should be passed. On output, A contains ! the transformed band matrix. The matrix A+T*I derived from the output ! parameters is similar to the input A+T*I to within rounding errors. ! Its last row and column are null as long as IERR is zero. ! ! Input/output, real ( kind = 8 ) T. On input, T specifies the shift (of eigenvalues) ! applied to the diagonal of A in forming the input matrix. What is ! actually determined is the eigenvalue nearest to T of A+T*I, where I ! is the identity matrix. On a subsequent call, the output value of T ! from the previous call should be passed if the next nearest eigenvalue ! is sought. On output, T contains the computed eigenvalue of A+T*I, ! as long as IERR is zero. ! ! Input/output, real ( kind = 8 ) R. On input for the first call, R should be ! specified as zero, and as its output value from the previous call ! on a subsequent call. It is used to determine when the last row and ! column of the transformed band matrix can be regarded as negligible. ! On output, R contains the maximum of its input value and the norm of the ! last column of the input matrix A. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, normal return. ! N, if the eigenvalue has not been determined after 30 iterations. ! implicit none integer ( kind = 4 ) mb integer ( kind = 4 ) n real ( kind = 8 ) a(n,mb) real ( kind = 8 ) f real ( kind = 8 ) g integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) ii integer ( kind = 4 ) ik integer ( kind = 4 ) imult integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) jk integer ( kind = 4 ) jm integer ( kind = 4 ) k integer ( kind = 4 ) kj integer ( kind = 4 ) kj1 integer ( kind = 4 ) kk integer ( kind = 4 ) km integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) m integer ( kind = 4 ) m1 integer ( kind = 4 ) m2 integer ( kind = 4 ) m21 integer ( kind = 4 ) m3 integer ( kind = 4 ) m31 integer ( kind = 4 ) m4 integer ( kind = 4 ) mk integer ( kind = 4 ) mn integer ( kind = 4 ) mz integer ( kind = 4 ) ni real ( kind = 8 ) pythag real ( kind = 8 ) q real ( kind = 8 ) r real ( kind = 8 ) rv(2*mb*mb+4*mb-3) real ( kind = 8 ) s real ( kind = 8 ) scale real ( kind = 8 ) t real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 ierr = 0 m1 = min ( mb, n ) m = m1 - 1 m2 = m + m m21 = m2 + 1 m3 = m21 + m m31 = m3 + 1 m4 = m31 + m2 mn = m + n mz = mb - m1 its = 0 ! ! Test for convergence. ! 40 continue g = a(n,mb) if ( m == 0 ) go to 360 f = 0.0D+00 do k = 1, m mk = k + mz f = f + abs ( a(n,mk) ) end do if ( its == 0 .and. r < f ) then r = f end if tst1 = r tst2 = tst1 + f if ( tst2 <= tst1 ) go to 360 if ( 30 <= its ) then ierr = n return end if its = its + 1 ! ! Form shift from bottom 2 by 2 minor. ! if ( f <= 0.25D+00 * r .or. its >= 5 ) then f = a(n,mb-1) if ( f /= 0.0D+00 ) then q = ( a(n-1,mb) - g ) / ( 2.0D+00 * f ) s = pythag ( q, 1.0D+00 ) g = g - f / ( q + sign ( s, q ) ) end if t = t + g a(1:n,mb) = a(1:n,mb) - g end if rv(m31:m4) = 0.0D+00 do ii = 1, mn i = ii - m ni = n - ii if ( ni < 0 ) go to 230 ! ! Form column of shifted matrix A-G*I. ! l = max ( 1, 2-i ) rv(1:m3) = 0.0D+00 do k = l, m1 km = k + m mk = k + mz rv(km) = a(ii,mk) end do ll = min ( m, ni ) do k = 1, ll km = k + m21 ik = ii + k mk = mb - k rv(km) = a(ik,mk) end do ! ! Pre-multiply with Householder reflections. ! ll = m2 imult = 0 ! ! Multiplication procedure. ! 140 continue kj = m4 - m1 do j = 1, ll kj = kj + m1 jm = j + m3 if ( rv(jm) /= 0.0D+00 ) then f = 0.0D+00 do k = 1, m1 kj = kj + 1 jk = j + k - 1 f = f + rv(kj) * rv(jk) end do f = f / rv(jm) kj = kj - m1 do k = 1, m1 kj = kj + 1 jk = j + k - 1 rv(jk) = rv(jk) - rv(kj) * f end do kj = kj - m1 end if end do if ( imult /= 0 ) go to 280 ! ! Householder reflection. ! f = rv(m21) s = 0.0D+00 rv(m4) = 0.0D+00 scale = sum ( abs ( rv(m21:m3) ) ) if ( scale == 0.0D+00 ) then go to 210 end if do k = m21, m3 s = s + ( rv(k) / scale )**2 end do s = scale * scale * s g = - sign ( sqrt ( s ), f ) rv(m21) = g rv(m4) = s - f * g kj = m4 + m2 * m1 + 1 rv(kj) = f - g do k = 2, m1 kj = kj + 1 km = k + m2 rv(kj) = rv(km) end do ! ! Save column of triangular factor R. ! 210 continue do k = l, m1 km = k + m mk = k + mz a(ii,mk) = rv(km) end do 230 continue l = max ( 1, m1+1-i ) if ( i <= 0 ) go to 300 ! ! Perform additional steps. ! rv(1:m21) = 0.0D+00 ll = min ( m1, ni+m1 ) ! ! Get row of triangular factor R. ! do kk = 1, ll k = kk - 1 km = k + m1 ik = i + k mk = mb - k rv(km) = a(ik,mk) end do ! ! Post-multiply with Householder reflections. ! ll = m1 imult = 1 go to 140 ! ! Store column of new a matrix. ! 280 continue do k = l, m1 mk = k + mz a(i,mk) = rv(k) end do ! ! Update Householder reflections. ! 300 continue if ( 1 < l ) then l = l - 1 end if kj1 = m4 + l * m1 do j = l, m2 jm = j + m3 rv(jm) = rv(jm+1) do k = 1, m1 kj1 = kj1 + 1 kj = kj1 - m1 rv(kj) = rv(kj1) end do end do end do go to 40 ! ! Convergence. ! 360 continue t = t + g a(1:n,mb) = a(1:n,mb) - g do k = 1, m1 mk = k + mz a(n,mk) = 0.0D+00 end do return end subroutine cbabk2 ( n, low, igh, scale, m, zr, zi ) !*****************************************************************************80 ! !! CBABK2 finds eigenvectors by undoing the CBAL transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a complex general ! matrix by back transforming those of the corresponding ! balanced matrix determined by CBAL. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, values determined by CBAL. ! ! Input, real ( kind = 8 ) SCALE(N), information determining the permutations ! and scaling factors used by CBAL. ! ! Input, integer ( kind = 4 ) M, the number of eigenvectors to be back transformed. ! ! Input/output, real ( kind = 8 ) ZR(N,M), ZI(N,M). On input, the real and imaginary ! parts, respectively, of the eigenvectors to be back transformed in ! their first M columns. On output, the transformed eigenvectors. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) i integer ( kind = 4 ) igh integer ( kind = 4 ) ii integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) low real ( kind = 8 ) s real ( kind = 8 ) scale(n) real ( kind = 8 ) zi(n,m) real ( kind = 8 ) zr(n,m) if ( m == 0 ) then return end if if ( igh /= low ) then do i = low, igh s = scale(i) zr(i,1:m) = zr(i,1:m) * s zi(i,1:m) = zi(i,1:m) * s end do end if do ii = 1, n i = ii if ( i < low .or. i > igh ) then if ( i < low ) then i = low - ii end if k = scale(i) if ( k /= i ) then do j = 1, m call r8_swap ( zr(i,j), zr(k,j) ) call r8_swap ( zi(i,j), zi(k,j) ) end do end if end if end do return end subroutine cbal ( n, ar, ai, low, igh, scale ) !*****************************************************************************80 ! !! CBAL balances a complex matrix before eigenvalue calculations. ! ! Discussion: ! ! This subroutine balances a complex matrix and isolates ! eigenvalues whenever possible. ! ! Suppose that the principal submatrix in rows low through igh ! has been balanced, that P(J) denotes the index interchanged ! with J during the permutation step, and that the elements ! of the diagonal matrix used are denoted by D(I,J). Then ! SCALE(J) = P(J), for J = 1,...,LOW-1 ! = D(J,J) J = LOW,...,IGH ! = P(J) J = IGH+1,...,N. ! The order in which the interchanges are made is N to IGH+1, ! then 1 to LOW-1. ! ! Note that 1 is returned for IGH if IGH is zero formally. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input/output, real ( kind = 8 ) AR(N,N), AI(N,N). On input, the real and ! imaginary parts of the complex matrix to be balanced. On output, ! the real and imaginary parts of the balanced matrix. ! ! Output, integer ( kind = 4 ) LOW, IGH, are values such that AR(I,J) and AI(I,J) ! are zero if I is greater than J and either J=1,...,LOW-1 or ! I=IGH+1,...,N. ! ! Output, real ( kind = 8 ) SCALE(N), information determining the ! permutations and scaling factors used. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) real ( kind = 8 ) b2 real ( kind = 8 ) c real ( kind = 8 ) f real ( kind = 8 ) g integer ( kind = 4 ) i integer ( kind = 4 ) iexc integer ( kind = 4 ) igh integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) low integer ( kind = 4 ) m logical noconv real ( kind = 8 ) r real ( kind = 8 ) radix real ( kind = 8 ) s real ( kind = 8 ) scale(n) radix = 16.0D+00 iexc = 0 j = 0 m = 0 b2 = radix * radix k = 1 l = n go to 100 20 continue scale(m) = j if ( j /= m ) then do i = 1, l call r8_swap ( ar(i,j), ar(i,m) ) call r8_swap ( ai(i,j), ai(i,m) ) end do do i = k, n call r8_swap ( ar(j,i), ar(m,i) ) call r8_swap ( ai(j,i), ai(m,i) ) end do end if if ( iexc == 2 ) then go to 130 end if ! ! Search for rows isolating an eigenvalue and push them down. ! 80 continue if ( l == 1 ) then go to 280 end if l = l - 1 100 continue do jj = 1, l j = l + 1 - jj do i = 1, l if ( i /= j ) then if ( ar(j,i) /= 0.0D+00 .or. ai(j,i) /= 0.0D+00 ) go to 120 end if end do m = l iexc = 1 go to 20 120 continue end do go to 140 ! ! Search for columns isolating an eigenvalue and push them left. ! 130 continue k = k + 1 140 continue do j = k, l do i = k, l if ( i /= j ) then if ( ar(i,j) /= 0.0D+00 .or. ai(i,j) /= 0.0D+00 ) go to 170 end if end do m = k iexc = 2 go to 20 170 continue end do ! ! Now balance the submatrix in rows k to l. ! scale(k:l) = 1.0D+00 ! ! Iterative loop for norm reduction. ! 190 continue noconv = .false. do i = k, l c = 0.0D+00 r = 0.0D+00 do j = k, l if ( j /= i ) then c = c + abs ( ar(j,i) ) + abs ( ai(j,i) ) r = r + abs ( ar(i,j) ) + abs ( ai(i,j) ) end if end do ! ! Guard against zero C or R due to underflow. ! if ( c == 0.0D+00 .or. r == 0.0D+00 ) go to 270 g = r / radix f = 1.0D+00 s = c + r do while ( c < g ) f = f * radix c = c * b2 end do g = r * radix do while ( c >= g ) f = f / radix c = c / b2 end do ! ! Now balance. ! if ( ( c + r ) / f < 0.95D+00 * s ) then g = 1.0D+00 / f scale(i) = scale(i) * f noconv = .true. ar(i,k:n) = ar(i,k:n) * g ai(i,k:n) = ai(i,k:n) * g ar(1:l,i) = ar(1:l,i) * f ai(1:l,i) = ai(1:l,i) * f end if 270 continue end do if ( noconv ) go to 190 280 continue low = k igh = l return end subroutine cdiv ( ar, ai, br, bi, cr, ci ) !*****************************************************************************80 ! !! CDIV emulates complex division, using real arithmetic. ! ! Discussion: ! ! This routine performs complex division: ! ! (CR,CI) = (AR,AI) / (BR,BI) ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, real ( kind = 8 ) AR, AI, the real and imaginary parts of the numerator. ! ! Input, real ( kind = 8 ) BR, BI, the real and imaginary parts of the denominator. ! ! Output, real ( kind = 8 ) CR, CI, the real and imaginary parts of the result. ! implicit none real ( kind = 8 ) ai real ( kind = 8 ) ais real ( kind = 8 ) ar real ( kind = 8 ) ars real ( kind = 8 ) bi real ( kind = 8 ) bis real ( kind = 8 ) br real ( kind = 8 ) brs real ( kind = 8 ) ci real ( kind = 8 ) cr real ( kind = 8 ) s s = abs ( br ) + abs ( bi ) ars = ar / s ais = ai / s brs = br / s bis = bi / s s = brs**2 + bis**2 cr = ( ars * brs + ais * bis ) / s ci = ( ais * brs - ars * bis ) / s return end subroutine cg ( n, ar, ai, wr, wi, matz, zr, zi, ierr ) !*****************************************************************************80 ! !! CG gets eigenvalues and eigenvectors of a complex general matrix. ! ! Discussion: ! ! This subroutine calls the recommended sequence of EISPACK subroutines ! to find the eigenvalues and eigenvectors (if desired) ! of a complex general matrix. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input/output, real ( kind = 8 ) AR(N,N), AI(N,N). On input, the real and ! imaginary parts of the complex matrix. On output, AR and AI ! have been overwritten by other information. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts ! of the eigenvalues. ! ! Input, integer ( kind = 4 ) MATZ, is 0 if only eigenvalues are desired, and ! nonzero if both eigenvalues and eigenvectors are to be computed. ! ! Output, real ( kind = 8 ) ZR(N,N), ZI(N,N), the real and imaginary parts, ! respectively, of the eigenvectors, if MATZ is not zero. ! ! Output, integer ( kind = 4 ) IERR, an error completion code described in the ! documentation for COMQR and COMQR2. The normal completion code is zero. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) real ( kind = 8 ) fv1(n) real ( kind = 8 ) fv2(n) real ( kind = 8 ) fv3(n) integer ( kind = 4 ) ierr integer ( kind = 4 ) is1 integer ( kind = 4 ) is2 integer ( kind = 4 ) matz real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) zi(n,n) real ( kind = 8 ) zr(n,n) call cbal ( n, ar, ai, is1, is2, fv1 ) call corth ( n, is1, is2, ar, ai, fv2, fv3 ) if ( matz == 0 ) then call comqr ( n, is1, is2, ar, ai, wr, wi, ierr ) if ( ierr /= 0 ) then return end if else call comqr2 ( n, is1, is2, fv2, fv3, ar, ai, wr, wi, zr, zi, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CG - Fatal error!' write ( *, '(a)' ) ' Nonzero error return from COMQR2.' return end if call cbabk2 ( n, is1, is2, fv1, n, zr, zi ) end if return end subroutine ch ( n, ar, ai, w, matz, zr, zi, ierr ) !*****************************************************************************80 ! !! CH gets eigenvalues and eigenvectors of a complex Hermitian matrix. ! ! Discussion: ! ! This subroutine calls the recommended sequence of subroutines from the ! EISPACK eigensystem package to find the eigenvalues and eigenvectors ! of a complex hermitian matrix. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input/output, real ( kind = 8 ) AR(N,N), AI(N,N). On input, the real and ! imaginary parts of the complex matrix. On output, AR and AI ! have been overwritten by other information. ! ! Output, real ( kind = 8 ) W(N), the eigenvalues in ascending order. ! ! Input, integer ( kind = 4 ) MATZ, is 0 if only eigenvalues are desired, and ! nonzero if both eigenvalues and eigenvectors are to be computed. ! ! Output, real ( kind = 8 ) ZR(N,N), ZI(N,N), the real and imaginary parts, ! respectively, of the eigenvectors, if MATZ is not zero. ! ! Output, integer ( kind = 4 ) IERR, an error completion code described in the ! documentation for TQLRAT and TQL2. The normal completion code is zero. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) real ( kind = 8 ) fm1(2,n) real ( kind = 8 ) fv1(n) real ( kind = 8 ) fv2(n) integer ( kind = 4 ) ierr integer ( kind = 4 ) matz real ( kind = 8 ) w(n) real ( kind = 8 ) zi(n,n) real ( kind = 8 ) zr(n,n) call htridi ( n, ar, ai, w, fv1, fv2, fm1 ) if ( matz == 0 ) then call tqlrat ( n, w, fv2, ierr ) else call r8mat_identity ( n, zr ) call tql2 ( n, w, fv1, zr, ierr ) if ( ierr /= 0 ) then return end if call htribk ( n, ar, ai, fm1, n, zr, zi ) end if return end subroutine cinvit ( n, ar, ai, wr, wi, select, mm, m, zr, zi, ierr ) !*****************************************************************************80 ! !! CINVIT gets eigenvectors from eigenvalues, for a complex Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds those eigenvectors of a complex upper ! Hessenberg matrix corresponding to specified eigenvalues, ! using inverse iteration. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, real ( kind = 8 ) AR(N,N), AI(N,N), the real and imaginary parts of ! the complex Hessenberg matrix. ! ! Input/output, real ( kind = 8 ) WR(N), WI(N). On input, the real and imaginary parts ! of the eigenvalues of the matrix. The eigenvalues must be stored in a ! manner identical to that of subroutine COMLR, which recognizes possible ! splitting of the matrix. On output, WR may have been altered since ! close eigenvalues are perturbed slightly in searching for independent ! eigenvectors. ! ! Input, logical SELECT(N), specifies the eigenvectors to be found. The ! eigenvector corresponding to the J-th eigenvalue is specified by ! setting SELECT(J) to TRUE. ! ! Input, integer ( kind = 4 ) MM, an upper bound for the number of eigenvectors ! to be found. ! ! Output, integer ( kind = 4 ) M, the number of eigenvectors actually found. ! ! Output, real ( kind = 8 ) ZR(N,MM), ZI(N,MM), the real and imaginary parts ! of the eigenvectors. The eigenvectors are normalized so that the ! component of largest magnitude is 1. ! Any vector which fails the acceptance test is set to zero. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! -(2*N+1), if more than MM eigenvectors have been specified, ! -K, if the iteration corresponding to the K-th value fails, ! -(N+K), if both error situations occur. ! implicit none integer ( kind = 4 ) mm integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) real ( kind = 8 ) eps3 real ( kind = 8 ) growto integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) ii real ( kind = 8 ) ilambd integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) km1 integer ( kind = 4 ) m integer ( kind = 4 ) mp real ( kind = 8 ) norm real ( kind = 8 ) normv real ( kind = 8 ) pythag real ( kind = 8 ) rlambd real ( kind = 8 ) rm1(n,n) real ( kind = 8 ) rm2(n,n) real ( kind = 8 ) rv1(n) real ( kind = 8 ) rv2(n) integer ( kind = 4 ) s logical select(n) integer ( kind = 4 ) uk real ( kind = 8 ) ukroot real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) zi(n,mm) real ( kind = 8 ) zr(n,mm) ierr = 0 uk = 0 s = 1 do k = 1, n if ( .not. select(k) ) then cycle end if if ( s > mm ) go to 1000 if ( uk >= k ) go to 200 ! ! Check for possible splitting. ! do uk = k, n - 1 if ( ar(uk+1,uk) == 0.0D+00 .and. ai(uk+1,uk) == 0.0D+00 ) then exit end if end do ! ! Compute infinity norm of leading UK by UK (Hessenberg) matrix. ! norm = 0.0D+00 mp = 1 do i = 1, uk x = 0.0D+00 do j = mp, uk x = x + pythag ( ar(i,j), ai(i,j) ) end do norm = max ( norm, x ) mp = i end do ! ! EPS3 replaces zero pivot in decomposition ! and close roots are modified by EPS3. ! if ( norm == 0.0D+00 ) norm = 1.0D+00 eps3 = abs ( norm ) * epsilon ( eps3 ) ! ! GROWTO is the criterion for growth. ! ukroot = uk ukroot = sqrt ( ukroot ) growto = 0.1D+00 / ukroot 200 continue rlambd = wr(k) ilambd = wi(k) if ( k == 1 ) go to 280 km1 = k - 1 go to 240 ! ! Perturb eigenvalue if it is close to any previous eigenvalue. ! 220 continue rlambd = rlambd + eps3 240 continue do ii = 1, km1 i = k - ii if ( select(i) .and. abs ( wr(i)-rlambd) < eps3 .and. & abs ( wi(i)-ilambd) < eps3 ) then go to 220 end if end do wr(k) = rlambd ! ! Form upper Hessenberg (ar,ai)-(rlambd,ilambd) * I ! and initial complex vector. ! 280 continue mp = 1 do i = 1, uk do j = mp, uk rm1(i,j) = ar(i,j) rm2(i,j) = ai(i,j) end do rm1(i,i) = rm1(i,i) - rlambd rm2(i,i) = rm2(i,i) - ilambd mp = i rv1(i) = eps3 end do ! ! Triangular decomposition with interchanges, replacing zero pivots by eps3. ! do i = 2, uk mp = i - 1 if ( pythag ( rm1(i,mp), rm2(i,mp) ) > & pythag ( rm1(mp,mp),rm2(mp,mp) ) ) then do j = mp, uk call r8_swap ( rm1(i,j), rm1(mp,j) ) call r8_swap ( rm2(i,j), rm2(mp,j) ) end do end if if ( rm1(mp,mp) == 0.0D+00 .and. rm2(mp,mp) == 0.0D+00 ) then rm1(mp,mp) = eps3 end if call cdiv ( rm1(i,mp), rm2(i,mp), rm1(mp,mp), rm2(mp,mp), x, y ) if ( x /= 0.0D+00 .or. y /= 0.0D+00 ) then do j = i, uk rm1(i,j) = rm1(i,j) - x * rm1(mp,j) + y * rm2(mp,j) rm2(i,j) = rm2(i,j) - x * rm2(mp,j) - y * rm1(mp,j) end do end if end do if ( rm1(uk,uk) == 0.0D+00 .and. rm2(uk,uk) == 0.0D+00 ) then rm1(uk,uk) = eps3 end if its = 0 ! ! Back substitution. ! 660 continue do ii = 1, uk i = uk + 1 - ii x = rv1(i) y = 0.0D+00 do j = i+1, uk x = x - rm1(i,j) * rv1(j) + rm2(i,j) * rv2(j) y = y - rm1(i,j) * rv2(j) - rm2(i,j) * rv1(j) end do call cdiv ( x, y, rm1(i,i), rm2(i,i), rv1(i), rv2(i) ) end do ! ! Acceptance test for eigenvector and normalization. ! its = its + 1 norm = 0.0D+00 normv = 0.0D+00 do i = 1, uk x = pythag ( rv1(i), rv2(i) ) if ( normv < x ) then normv = x j = i end if norm = norm + x end do if ( norm < growto ) go to 840 ! ! Accept vector. ! x = rv1(j) y = rv2(j) do i = 1, uk call cdiv ( rv1(i), rv2(i), x, y, zr(i,s), zi(i,s) ) end do if ( uk == n ) then go to 940 end if j = uk + 1 go to 900 ! ! Choose a new starting vector. ! 840 continue if ( its < uk ) then x = ukroot y = eps3 / ( x + 1.0D+00 ) rv1(1) = eps3 rv1(2:uk) = y j = uk - its + 1 rv1(j) = rv1(j) - eps3 * x go to 660 end if ! ! Error: unaccepted eigenvector. ! 880 continue j = 1 ierr = -k ! ! Set remaining vector components to zero. ! 900 continue zr(j:n,s) = 0.0D+00 zi(j:n,s) = 0.0D+00 940 continue s = s + 1 end do go to 1001 ! ! Set error: underestimate of eigenvector space required. ! 1000 continue if ( ierr /= 0 ) ierr = ierr - n if ( ierr == 0 ) ierr = -(2 * n + 1) 1001 continue m = s - 1 return end subroutine combak ( n, low, igh, ar, ai, int, m, zr, zi ) !*****************************************************************************80 ! !! COMBAK determines eigenvectors by undoing the COMHES transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a complex general ! matrix by back transforming those of the corresponding ! upper Hessenberg matrix determined by COMHES. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = to the order of the matrix. ! ! Input, real ( kind = 8 ) AR(N,IGH), AI(N,IGH), the multipliers which were used in the ! reduction by COMHES in their lower triangles below the subdiagonal. ! ! Input, integer ( kind = 4 ) INT(IGH), information on the rows and columns interchanged ! in the reduction by COMHES. ! ! Input, integer ( kind = 4 ) M, the number of eigenvectors to be back transformed. ! ! Input/output, real ( kind = 8 ) ZR(N,M), ZI(N,M). On input, the real and imaginary ! parts of the eigenvectors to be back transformed. On output, the real ! and imaginary parts of the transformed eigenvectors. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) ai(n,igh) real ( kind = 8 ) ar(n,igh) integer ( kind = 4 ) i integer ( kind = 4 ) int(igh) integer ( kind = 4 ) j integer ( kind = 4 ) la integer ( kind = 4 ) low integer ( kind = 4 ) mm integer ( kind = 4 ) mp real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) zi(n,m) real ( kind = 8 ) zr(n,m) if ( m == 0 ) then return end if la = igh - 1 if ( igh - 1 < low + 1 ) then return end if do mm = low + 1, la mp = low + igh - mm do i = mp+1, igh xr = ar(i,mp-1) xi = ai(i,mp-1) if ( xr /= 0.0D+00 .or. xi /= 0.0D+00 ) then zr(i,1:m) = zr(i,1:m) + xr * zr(mp,1:m) - xi * zi(mp,1:m) zi(i,1:m) = zi(i,1:m) + xr * zi(mp,1:m) + xi * zr(mp,1:m) end if end do i = int(mp) if ( i /= mp ) then do j = 1, m call r8_swap ( zr(i,j), zr(mp,j) ) call r8_swap ( zi(i,j), zi(mp,j) ) end do end if end do return end subroutine comhes ( n, low, igh, ar, ai, int ) !*****************************************************************************80 ! !! COMHES transforms a complex general matrix to upper Hessenberg form. ! ! Discussion: ! ! Given a complex general matrix, this subroutine ! reduces a submatrix situated in rows and columns ! LOW through IGH to upper Hessenberg form by ! stabilized elementary similarity transformations. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input/output, real ( kind = 8 ) AR(N,N), AI(N,N). On input, the real and imaginary ! parts of the complex input matrix. On output, the real and imaginary ! parts of the Hessenberg matrix. The multipliers which were used in the ! reduction are stored in the remaining triangles under the ! Hessenberg matrix. ! ! Output, integer ( kind = 4 ) INT(IGH), information on the rows and columns ! interchanged in the reduction. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) int(igh) integer ( kind = 4 ) j integer ( kind = 4 ) la integer ( kind = 4 ) low integer ( kind = 4 ) m real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr la = igh - 1 do m = low + 1, la xr = 0.0D+00 xi = 0.0D+00 i = m do j = m, igh if ( abs ( ar(j,m-1) ) + abs ( ai(j,m-1) ) > & abs ( xr ) + abs ( xi ) ) then xr = ar(j,m-1) xi = ai(j,m-1) i = j end if end do int(m) = i ! ! Interchange rows and columns of AR and AI. ! if ( i /= m ) then do j = m-1, n call r8_swap ( ar(i,j), ar(m,j) ) call r8_swap ( ai(i,j), ai(m,j) ) end do do j = 1, igh call r8_swap ( ar(j,i), ar(j,m) ) call r8_swap ( ai(j,i), ai(j,m) ) end do end if if ( xr /= 0.0D+00 .or. xi /= 0.0D+00 ) then do i = m+1, igh yr = ar(i,m-1) yi = ai(i,m-1) if ( yr /= 0.0D+00 .or. yi /= 0.0D+00 ) then call cdiv ( yr, yi, xr, xi, yr, yi ) ar(i,m-1) = yr ai(i,m-1) = yi do j = m, n ar(i,j) = ar(i,j) - yr * ar(m,j) + yi * ai(m,j) ai(i,j) = ai(i,j) - yr * ai(m,j) - yi * ar(m,j) end do ar(1:igh,m) = ar(1:igh,m) + yr * ar(1:igh,i) - yi * ai(1:igh,i) ai(1:igh,m) = ai(1:igh,m) + yr * ai(1:igh,i) + yi * ar(1:igh,i) end if end do end if end do return end subroutine comlr ( n, low, igh, hr, hi, wr, wi, ierr ) !*****************************************************************************80 ! !! COMLR gets all eigenvalues of a complex upper Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalues of a complex upper Hessenberg ! matrix by the modified LR method. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input/output, real ( kind = 8 ) HR(N,N), HI(N,N). On input, the real and imaginary ! parts of the complex upper Hessenberg matrix. Their lower triangles ! below the subdiagonal contain the multipliers which were used in the ! reduction by COMHES if performed. On output, the upper Hessenberg ! portions of HR and HI have been destroyed. Therefore, they must be ! saved before calling COMLR if subsequent calculation of eigenvectors ! is to be performed. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts of the ! eigenvalues. If an error exit is made, the eigenvalues should be correct ! for indices IERR+1,...,N. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! J, if the limit of 30*N iterations is exhausted while the J-th ! eigenvalue is being sought. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) en integer ( kind = 4 ) enm1 real ( kind = 8 ) hi(n,n) real ( kind = 8 ) hr(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) igh integer ( kind = 4 ) itn integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) low integer ( kind = 4 ) m integer ( kind = 4 ) mm real ( kind = 8 ) si real ( kind = 8 ) sr real ( kind = 8 ) ti real ( kind = 8 ) tr real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr real ( kind = 8 ) zzi real ( kind = 8 ) zzr ierr = 0 ! ! Store roots isolated by CBAL. ! do i = 1, n if ( i < low .or. i > igh ) then wr(i) = hr(i,i) wi(i) = hi(i,i) end if end do en = igh tr = 0.0D+00 ti = 0.0D+00 itn = 30 * n ! ! Search for next eigenvalue. ! 220 continue if ( en < low ) then return end if its = 0 enm1 = en - 1 ! ! Look for single small sub-diagonal element. ! 240 continue do ll = low, en l = en + low - ll if ( l == low ) go to 300 tst1 = abs ( hr(l-1,l-1) ) + abs ( hi(l-1,l-1) ) + abs ( hr(l,l) ) & + abs ( hi(l,l) ) tst2 = tst1 + abs ( hr(l,l-1) ) + abs ( hi(l,l-1) ) if ( tst2 == tst1) go to 300 end do ! ! Form shift. ! 300 continue if ( l == en ) then go to 660 end if if ( itn == 0 ) then ierr = en return end if if ( its == 10 .or. its == 20 ) go to 320 sr = hr(en,en) si = hi(en,en) xr = hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1) xi = hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1) if ( xr == 0.0D+00 .and. xi == 0.0D+00 ) go to 340 yr = (hr(enm1,enm1) - sr) / 2.0D+00 yi = (hi(enm1,enm1) - si) / 2.0D+00 call csroot ( yr**2-yi**2+xr, 2.0D+00*yr*yi+xi, zzr, zzi ) if ( yr * zzr + yi * zzi < 0.0D+00 ) then zzr = -zzr zzi = -zzi end if call cdiv ( xr, xi, yr+zzr, yi+zzi, xr, xi ) sr = sr - xr si = si - xi go to 340 ! ! Form exceptional shift. ! 320 continue sr = abs ( hr(en,enm1) ) + abs ( hr(enm1,en-2) ) si = abs ( hi(en,enm1) ) + abs ( hi(enm1,en-2) ) 340 continue do i = low, en hr(i,i) = hr(i,i) - sr hi(i,i) = hi(i,i) - si end do tr = tr + sr ti = ti + si its = its + 1 itn = itn - 1 ! ! Look for two consecutive small sub-diagonal elements. ! xr = abs ( hr(enm1,enm1) ) + abs ( hi(enm1,enm1) ) yr = abs ( hr(en,enm1) ) + abs ( hi(en,enm1) ) zzr = abs ( hr(en,en) ) + abs ( hi(en,en) ) do mm = l, enm1 m = enm1 + l - mm if ( m == l ) then exit end if yi = yr yr = abs ( hr(m,m-1) ) + abs ( hi(m,m-1) ) xi = zzr zzr = xr xr = abs ( hr(m-1,m-1) ) + abs ( hi(m-1,m-1) ) tst1 = zzr / yi * (zzr + xr + xi) tst2 = tst1 + yr if ( tst2 == tst1 ) then exit end if end do ! ! Triangular decomposition H=L*R. ! do i = m+1, en xr = hr(i-1,i-1) xi = hi(i-1,i-1) yr = hr(i,i-1) yi = hi(i,i-1) if ( abs ( xr ) + abs ( xi ) >= abs ( yr ) + abs ( yi ) ) go to 460 ! ! Interchange rows of HR and HI. ! do j = i-1, en call r8_swap ( hr(i-1,j), hr(i,j) ) call r8_swap ( hi(i-1,j), hi(i,j) ) end do call cdiv ( xr, xi, yr, yi, zzr, zzi ) wr(i) = 1.0D+00 go to 480 460 continue call cdiv ( yr, yi, xr, xi, zzr, zzi ) wr(i) = -1.0D+00 480 continue hr(i,i-1) = zzr hi(i,i-1) = zzi do j = i, en hr(i,j) = hr(i,j) - zzr * hr(i-1,j) + zzi * hi(i-1,j) hi(i,j) = hi(i,j) - zzr * hi(i-1,j) - zzi * hr(i-1,j) end do end do ! ! Composition R*L=H. ! do j = m+1, en xr = hr(j,j-1) xi = hi(j,j-1) hr(j,j-1) = 0.0D+00 hi(j,j-1) = 0.0D+00 ! ! Interchange columns of HR and HI, if necessary. ! if ( wr(j) > 0.0D+00 ) then do i = l, j call r8_swap ( hr(i,j-1), hr(i,j) ) call r8_swap ( hi(i,j-1), hi(i,j) ) end do end if do i = l, j hr(i,j-1) = hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j) hi(i,j-1) = hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j) end do end do go to 240 ! ! A root found. ! 660 continue wr(en) = hr(en,en) + tr wi(en) = hi(en,en) + ti en = enm1 go to 220 end subroutine comlr2 ( n, low, igh, int, hr, hi, wr, wi, zr, zi, ierr ) !*****************************************************************************80 ! !! COMLR2 gets eigenvalues/vectors of a complex upper Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalues and eigenvectors of a complex ! upper Hessenberg matrix by the modified LR method. The eigenvectors ! of a complex general matrix can also be found if COMHES has been used ! to reduce this general matrix to Hessenberg form. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input, integer ( kind = 4 ) INT(IGH), information on the rows and columns interchanged ! in the reduction by COMHES, if performed. If the eigenvectors of the ! Hessenberg matrix are desired, set INT(J)=J for these elements. ! ! Input/output, real ( kind = 8 ) HR(N,N), HI(N,N). On input, the real and imaginary ! parts of the complex upper Hessenberg matrix. Their lower triangles ! below the subdiagonal contain the multipliers which were used in the ! reduction by COMHES, if performed. If the eigenvectors of the Hessenberg ! matrix are desired, these elements must be set to zero. On output, ! the upper Hessenberg portions of HR and HI have been destroyed, but the ! location HR(1,1) contains the norm of the triangularized matrix. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts of the ! eigenvalues. If an error exit is made, the eigenvalues should be ! correct for indices IERR+1,...,N. ! ! Output, real ( kind = 8 ) ZR(N,N), ZI(N,N), the real and imaginary parts of the ! eigenvectors. The eigenvectors are unnormalized. If an error exit ! is made, none of the eigenvectors has been found. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! J, if the limit of 30*N iterations is exhausted while the J-th ! eigenvalue is being sought. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) en integer ( kind = 4 ) enm1 real ( kind = 8 ) hi(n,n) real ( kind = 8 ) hr(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) iend integer ( kind = 4 ) ierr integer ( kind = 4 ) igh integer ( kind = 4 ) ii integer ( kind = 4 ) int(igh) integer ( kind = 4 ) itn integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) low integer ( kind = 4 ) m integer ( kind = 4 ) mm integer ( kind = 4 ) nn real ( kind = 8 ) norm real ( kind = 8 ) si real ( kind = 8 ) sr real ( kind = 8 ) ti real ( kind = 8 ) tr real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr real ( kind = 8 ) zi(n,n) real ( kind = 8 ) zr(n,n) real ( kind = 8 ) zzi real ( kind = 8 ) zzr ierr = 0 ! ! Initialize the eigenvector matrix. ! call r8mat_identity ( n, zr ) zi(1:n,1:n) = 0.0D+00 ! ! Form the matrix of accumulated transformations from the information left ! by COMHES. ! iend = igh - low - 1 do ii = 1, iend i = igh - ii do k = i+1, igh zr(k,i) = hr(k,i-1) zi(k,i) = hi(k,i-1) end do j = int(i) if ( i /= j ) then do k = i, igh zr(i,k) = zr(j,k) zi(i,k) = zi(j,k) zr(j,k) = 0.0D+00 zi(j,k) = 0.0D+00 end do zr(j,i) = 1.0D+00 end if end do ! ! Store roots isolated by CBAL. ! do i = 1, n if ( i < low .or. i > igh ) then wr(i) = hr(i,i) wi(i) = hi(i,i) end if end do en = igh tr = 0.0D+00 ti = 0.0D+00 itn = 30 * n ! ! Search for next eigenvalue. ! 220 continue if ( en < low ) then go to 680 end if its = 0 enm1 = en - 1 ! ! Look for single small sub-diagonal element. ! 240 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if tst1 = abs ( hr(l-1,l-1) ) + abs ( hi(l-1,l-1) ) + abs ( hr(l,l) ) & + abs ( hi(l,l) ) tst2 = tst1 + abs ( hr(l,l-1) ) + abs ( hi(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! if ( l == en ) go to 660 if ( itn == 0 ) go to 1000 if ( its == 10 .or. its == 20 ) go to 320 sr = hr(en,en) si = hi(en,en) xr = hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1) xi = hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1) if ( xr == 0.0D+00 .and. xi == 0.0D+00 ) go to 340 yr = (hr(enm1,enm1) - sr) / 2.0D+00 yi = (hi(enm1,enm1) - si) / 2.0D+00 call csroot ( yr**2-yi**2+xr, 2.0D+00*yr*yi+xi, zzr, zzi ) if ( yr * zzr + yi * zzi < 0.0D+00 ) then zzr = -zzr zzi = -zzi end if call cdiv ( xr, xi, yr+zzr, yi+zzi, xr, xi ) sr = sr - xr si = si - xi go to 340 ! ! Form exceptional shift. ! 320 continue sr = abs ( hr(en,enm1) ) + abs ( hr(enm1,en-2) ) si = abs ( hi(en,enm1) ) + abs ( hi(enm1,en-2) ) 340 continue do i = low, en hr(i,i) = hr(i,i) - sr hi(i,i) = hi(i,i) - si end do tr = tr + sr ti = ti + si its = its + 1 itn = itn - 1 ! ! Look for two consecutive small sub-diagonal elements. ! xr = abs ( hr(enm1,enm1) ) + abs ( hi(enm1,enm1) ) yr = abs ( hr(en,enm1) ) + abs ( hi(en,enm1) ) zzr = abs ( hr(en,en) ) + abs ( hi(en,en) ) do mm = l, enm1 m = enm1 + l - mm if ( m == l ) then exit end if yi = yr yr = abs ( hr(m,m-1) ) + abs ( hi(m,m-1) ) xi = zzr zzr = xr xr = abs ( hr(m-1,m-1) ) + abs ( hi(m-1,m-1) ) tst1 = zzr / yi * (zzr + xr + xi) tst2 = tst1 + yr if ( tst2 == tst1 ) then exit end if end do ! ! Triangular decomposition H=L*R. ! do i = m+1, en xr = hr(i-1,i-1) xi = hi(i-1,i-1) yr = hr(i,i-1) yi = hi(i,i-1) if ( abs ( xr ) + abs ( xi) >= abs ( yr ) + abs ( yi ) ) go to 460 ! ! Interchange rows of HR and HI. ! do j = i-1, n call r8_swap ( hr(i-1,j), hr(i,j) ) call r8_swap ( hi(i-1,j), hi(i,j) ) end do call cdiv ( xr, xi, yr, yi, zzr, zzi ) wr(i) = 1.0D+00 go to 480 460 continue call cdiv ( yr, yi, xr, xi, zzr, zzi ) wr(i) = -1.0D+00 480 continue hr(i,i-1) = zzr hi(i,i-1) = zzi do j = i, n hr(i,j) = hr(i,j) - zzr * hr(i-1,j) + zzi * hi(i-1,j) hi(i,j) = hi(i,j) - zzr * hi(i-1,j) - zzi * hr(i-1,j) end do end do ! ! Composition R*L=H. ! do j = m+1, en xr = hr(j,j-1) xi = hi(j,j-1) hr(j,j-1) = 0.0D+00 hi(j,j-1) = 0.0D+00 ! ! Interchange columns of HR, HI, ZR, and ZI. ! if ( wr(j) > 0.0D+00 ) then do i = 1, j call r8_swap ( hr(i,j-1), hr(i,j) ) call r8_swap ( hi(i,j-1), hi(i,j) ) end do do i = low, igh call r8_swap ( zr(i,j-1), zr(i,j) ) call r8_swap ( zi(i,j-1), zi(i,j) ) end do end if do i = 1, j hr(i,j-1) = hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j) hi(i,j-1) = hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j) end do ! ! Accumulate transformations. ! do i = low, igh zr(i,j-1) = zr(i,j-1) + xr * zr(i,j) - xi * zi(i,j) zi(i,j-1) = zi(i,j-1) + xr * zi(i,j) + xi * zr(i,j) end do end do go to 240 ! ! A root found. ! 660 continue hr(en,en) = hr(en,en) + tr wr(en) = hr(en,en) hi(en,en) = hi(en,en) + ti wi(en) = hi(en,en) en = enm1 go to 220 ! ! All roots found. ! Backsubstitute to find vectors of upper triangular form. ! 680 continue norm = 0.0D+00 do i = 1, n do j = i, n tr = abs ( hr(i,j) ) + abs ( hi(i,j) ) if ( tr > norm ) norm = tr end do end do hr(1,1) = norm if ( n == 1 ) then return end if if ( norm == 0.0D+00 ) then return end if do nn = 2, n en = n + 2 - nn xr = wr(en) xi = wi(en) hr(en,en) = 1.0D+00 hi(en,en) = 0.0D+00 enm1 = en - 1 do ii = 1, enm1 i = en - ii zzr = 0.0D+00 zzi = 0.0D+00 do j = i+1, en zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) end do yr = xr - wr(i) yi = xi - wi(i) if ( yr == 0.0D+00 .and. yi == 0.0D+00 ) then tst1 = norm yr = tst1 do yr = 0.01D+00 * yr tst2 = norm + yr if ( tst2 <= tst1 ) then exit end if end do end if call cdiv ( zzr, zzi, yr, yi, hr(i,en), hi(i,en) ) ! ! Overflow control. ! tr = abs ( hr(i,en) ) + abs ( hi(i,en) ) if ( tr /= 0.0D+00 ) then tst1 = tr tst2 = tst1 + 1.0D+00 / tst1 if ( tst2 <= tst1 ) then hr(i:en,en) = hr(i:en,en) / tr hi(i:en,en) = hi(i:en,en) / tr end if end if end do end do ! ! End backsubstitution. ! enm1 = n - 1 ! ! Vectors of isolated roots. ! do i = 1, n - 1 if ( i < low .or. i > igh ) then zr(i,i+1:n) = hr(i,i+1:n) zi(i,i+1:n) = hi(i,i+1:n) end if end do ! ! Multiply by transformation matrix to give vectors of original full matrix. ! do jj = low, n - 1 j = n + low - jj m = min ( j, igh ) do i = low, igh zzr = 0.0D+00 zzi = 0.0D+00 do k = low, m zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) end do zr(i,j) = zzr zi(i,j) = zzi end do end do return ! ! Set error: all eigenvalues have not converged after 30*N iterations. ! 1000 continue ierr = en return end subroutine comqr ( n, low, igh, hr, hi, wr, wi, ierr ) !*****************************************************************************80 ! !! COMQR gets eigenvalues of a complex upper Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalues of a complex ! upper Hessenberg matrix by the QR method. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input/output, real ( kind = 8 ) HR(N,N), HI(N,N). On input, the real and imaginary ! parts of the complex upper Hessenberg matrix. Their lower triangles ! below the subdiagonal contain information about the unitary ! transformations used in the reduction by CORTH, if performed. On output, ! the upper Hessenberg portions of HR and HI have been destroyed. ! Therefore, they must be saved before calling COMQR if subsequent ! calculation of eigenvectors is to be performed. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts of the ! eigenvalues. If an error exit is made, the eigenvalues should be ! correct for indices IERR+1,...,N. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! J, if the limit of 30*N iterations is exhausted while the J-th ! eigenvalue is being sought. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) en integer ( kind = 4 ) enm1 real ( kind = 8 ) hi(n,n) real ( kind = 8 ) hr(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) igh integer ( kind = 4 ) itn integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) low real ( kind = 8 ) norm real ( kind = 8 ) pythag real ( kind = 8 ) si real ( kind = 8 ) sr real ( kind = 8 ) ti real ( kind = 8 ) tr real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr real ( kind = 8 ) zzi real ( kind = 8 ) zzr ierr = 0 ! ! Create real subdiagonal elements. ! l = low + 1 do i = l, igh ll = min ( i+1, igh ) if ( hi(i,i-1) /= 0.0D+00 ) then norm = pythag ( hr(i,i-1), hi(i,i-1) ) yr = hr(i,i-1) / norm yi = hi(i,i-1) / norm hr(i,i-1) = norm hi(i,i-1) = 0.0D+00 do j = i, igh si = yr * hi(i,j) - yi * hr(i,j) hr(i,j) = yr * hr(i,j) + yi * hi(i,j) hi(i,j) = si end do do j = low, ll si = yr * hi(j,i) + yi * hr(j,i) hr(j,i) = yr * hr(j,i) - yi * hi(j,i) hi(j,i) = si end do end if end do ! ! Store roots isolated by CBAL. ! do i = 1, n if ( i < low .or. i > igh ) then wr(i) = hr(i,i) wi(i) = hi(i,i) end if end do en = igh tr = 0.0D+00 ti = 0.0D+00 itn = 30 * n ! ! Search for next eigenvalue. ! 220 continue if ( en < low ) then return end if its = 0 enm1 = en - 1 ! ! Look for single small sub-diagonal element. ! 240 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if tst1 = abs ( hr(l-1,l-1) ) + abs ( hi(l-1,l-1) ) + abs ( hr(l,l) ) & + abs ( hi(l,l) ) tst2 = tst1 + abs ( hr(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! if ( l == en ) then go to 660 end if if ( itn == 0 ) go to 1000 if ( its == 10 .or. its == 20 ) go to 320 sr = hr(en,en) si = hi(en,en) xr = hr(enm1,en) * hr(en,enm1) xi = hi(enm1,en) * hr(en,enm1) if ( xr == 0.0D+00 .and. xi == 0.0D+00 ) go to 340 yr = (hr(enm1,enm1) - sr) / 2.0D+00 yi = (hi(enm1,enm1) - si) / 2.0D+00 call csroot ( yr**2-yi**2+xr, 2.0D+00*yr*yi+xi, zzr, zzi ) if ( yr * zzr + yi * zzi < 0.0D+00 ) then zzr = -zzr zzi = -zzi end if call cdiv ( xr, xi, yr+zzr, yi+zzi, xr, xi ) sr = sr - xr si = si - xi go to 340 ! ! Form exceptional shift. ! 320 continue sr = abs ( hr(en,enm1) ) + abs ( hr(enm1,en-2) ) si = 0.0D+00 340 continue do i = low, en hr(i,i) = hr(i,i) - sr hi(i,i) = hi(i,i) - si end do tr = tr + sr ti = ti + si its = its + 1 itn = itn - 1 ! ! Reduce to triangle (rows). ! do i = l+1, en sr = hr(i,i-1) hr(i,i-1) = 0.0D+00 norm = pythag ( pythag ( hr(i-1,i-1), hi(i-1,i-1) ), sr ) xr = hr(i-1,i-1) / norm wr(i-1) = xr xi = hi(i-1,i-1) / norm wi(i-1) = xi hr(i-1,i-1) = norm hi(i-1,i-1) = 0.0D+00 hi(i,i-1) = sr / norm do j = i, en yr = hr(i-1,j) yi = hi(i-1,j) zzr = hr(i,j) zzi = hi(i,j) hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi end do end do si = hi(en,en) if ( si /= 0.0D+00 ) then norm = pythag ( hr(en,en), si ) sr = hr(en,en) / norm si = si / norm hr(en,en) = norm hi(en,en) = 0.0D+00 end if ! ! Inverse operation (columns). ! do j = l+1, en xr = wr(j-1) xi = wi(j-1) do i = l, j yr = hr(i,j-1) yi = 0.0D+00 zzr = hr(i,j) zzi = hi(i,j) if ( i /= j ) then yi = hi(i,j-1) hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi end if hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi end do end do if ( si /= 0.0D+00 ) then do i = l, en yr = hr(i,en) yi = hi(i,en) hr(i,en) = sr * yr - si * yi hi(i,en) = sr * yi + si * yr end do end if go to 240 ! ! A root found. ! 660 continue wr(en) = hr(en,en) + tr wi(en) = hi(en,en) + ti en = enm1 go to 220 ! ! Set error: all eigenvalues have not converged after 30*n iterations. ! 1000 continue ierr = en return end subroutine comqr2 ( n, low, igh, ortr, orti, hr, hi, wr, wi, zr, zi, ierr ) !*****************************************************************************80 ! !! COMQR2 gets eigenvalues/vectors of a complex upper Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalues and eigenvectors ! of a complex upper Hessenberg matrix by the QR ! method. The eigenvectors of a complex general matrix ! can also be found if CORTH has been used to reduce ! this general matrix to Hessenberg form. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input/output, real ( kind = 8 ) ORTR(N), ORTI(N). On input, information about the ! unitary transformations used in the reduction by CORTH, if performed. ! If the eigenvectors of the Hessenberg matrix are desired, set ORTR(J) and ! ORTI(J) to 0.0D+00 for these elements. On output, these arrays ! have been overwritten. ! ! Input/output, real ( kind = 8 ) HR(N,N), HI(N,N). On input, the real and imaginary ! parts of the complex upper Hessenberg matrix. Their lower triangles ! below the subdiagonal contain further information about the ! transformations which were used in the reduction by CORTH, if performed. ! If the eigenvectors of the Hessenberg matrix are desired, these elements ! may be arbitrary. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts of the ! eigenvalues. If an error exit is made, the eigenvalues should be ! correct for indices IERR+1,...,N. ! ! Output, real ( kind = 8 ) ZR(N,N), ZI(N,N), the real and imaginary parts of the ! eigenvectors. The eigenvectors are unnormalized. If an error exit ! is made, none of the eigenvectors has been found. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! J, if the limit of 30*N iterations is exhausted while the J-th ! eigenvalue is being sought. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) n integer ( kind = 4 ) en integer ( kind = 4 ) enm1 real ( kind = 8 ) hi(n,n) real ( kind = 8 ) hr(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) iend integer ( kind = 4 ) ierr integer ( kind = 4 ) ii integer ( kind = 4 ) itn integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) low integer ( kind = 4 ) m integer ( kind = 4 ) nn real ( kind = 8 ) norm real ( kind = 8 ) orti(igh) real ( kind = 8 ) ortr(igh) real ( kind = 8 ) pythag real ( kind = 8 ) si real ( kind = 8 ) sr real ( kind = 8 ) ti real ( kind = 8 ) tr real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr real ( kind = 8 ) zi(n,n) real ( kind = 8 ) zr(n,n) real ( kind = 8 ) zzi real ( kind = 8 ) zzr ierr = 0 ! ! Initialize eigenvector matrix. ! call r8mat_identity ( n, zr ) zi(1:n,1:n) = 0.0D+00 ! ! Form the matrix of accumulated transformations from the information ! left by CORTH. ! iend = igh - low - 1 if ( iend ) 180, 150, 105 105 continue do ii = 1, iend i = igh - ii if ( ortr(i) == 0.0D+00 .and. orti(i) == 0.0D+00 ) go to 140 if ( hr(i,i-1) == 0.0D+00 .and. hi(i,i-1) == 0.0D+00 ) go to 140 ! ! Norm below is negative of H formed in CORTH. ! norm = hr(i,i-1) * ortr(i) + hi(i,i-1) * orti(i) do k = i+1, igh ortr(k) = hr(k,i-1) orti(k) = hi(k,i-1) end do do j = i, igh sr = 0.0D+00 si = 0.0D+00 do k = i, igh sr = sr + ortr(k) * zr(k,j) + orti(k) * zi(k,j) si = si + ortr(k) * zi(k,j) - orti(k) * zr(k,j) end do sr = sr / norm si = si / norm do k = i, igh zr(k,j) = zr(k,j) + sr * ortr(k) - si * orti(k) zi(k,j) = zi(k,j) + sr * orti(k) + si * ortr(k) end do end do 140 continue end do ! ! Create real subdiagonal elements. ! 150 continue l = low + 1 do i = l, igh ll = min ( i+1, igh ) if ( hi(i,i-1) == 0.0D+00 ) then go to 170 end if norm = pythag ( hr(i,i-1), hi(i,i-1) ) yr = hr(i,i-1) / norm yi = hi(i,i-1) / norm hr(i,i-1) = norm hi(i,i-1) = 0.0D+00 do j = i, n si = yr * hi(i,j) - yi * hr(i,j) hr(i,j) = yr * hr(i,j) + yi * hi(i,j) hi(i,j) = si end do do j = 1, ll si = yr * hi(j,i) + yi * hr(j,i) hr(j,i) = yr * hr(j,i) - yi * hi(j,i) hi(j,i) = si end do do j = low, igh si = yr * zi(j,i) + yi * zr(j,i) zr(j,i) = yr * zr(j,i) - yi * zi(j,i) zi(j,i) = si end do 170 continue end do ! ! Store roots isolated by CBAL. ! 180 continue do i = 1, n if ( i < low .or. i > igh) then wr(i) = hr(i,i) wi(i) = hi(i,i) end if end do en = igh tr = 0.0D+00 ti = 0.0D+00 itn = 30 * n ! ! Search for next eigenvalue. ! 220 continue if ( en < low ) go to 680 its = 0 enm1 = en - 1 ! ! Look for single small sub-diagonal element. ! 240 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if tst1 = abs ( hr(l-1,l-1) ) + abs ( hi(l-1,l-1) ) + abs ( hr(l,l) ) & + abs ( hi(l,l) ) tst2 = tst1 + abs ( hr(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! if ( l == en ) go to 660 if ( itn == 0 ) go to 1000 if ( its == 10 .or. its == 20 ) go to 320 sr = hr(en,en) si = hi(en,en) xr = hr(enm1,en) * hr(en,enm1) xi = hi(enm1,en) * hr(en,enm1) if ( xr == 0.0D+00 .and. xi == 0.0D+00 ) go to 340 yr = ( hr(enm1,enm1) - sr ) / 2.0D+00 yi = ( hi(enm1,enm1) - si ) / 2.0D+00 call csroot ( yr**2-yi**2+xr, 2.0D+00*yr*yi+xi, zzr, zzi ) if ( yr * zzr + yi * zzi < 0.0D+00 ) then zzr = -zzr zzi = -zzi end if call cdiv ( xr, xi, yr+zzr, yi+zzi, xr, xi ) sr = sr - xr si = si - xi go to 340 ! ! Form exceptional shift. ! 320 continue sr = abs ( hr(en,enm1) ) + abs ( hr(enm1,en-2) ) si = 0.0D+00 340 continue do i = low, en hr(i,i) = hr(i,i) - sr hi(i,i) = hi(i,i) - si end do tr = tr + sr ti = ti + si its = its + 1 itn = itn - 1 ! ! Reduce to triangle (rows). ! do i = l+1, en sr = hr(i,i-1) hr(i,i-1) = 0.0D+00 norm = pythag ( pythag ( hr(i-1,i-1), hi(i-1,i-1) ), sr ) xr = hr(i-1,i-1) / norm wr(i-1) = xr xi = hi(i-1,i-1) / norm wi(i-1) = xi hr(i-1,i-1) = norm hi(i-1,i-1) = 0.0D+00 hi(i,i-1) = sr / norm do j = i, n yr = hr(i-1,j) yi = hi(i-1,j) zzr = hr(i,j) zzi = hi(i,j) hr(i-1,j) = xr * yr + xi * yi + hi(i,i-1) * zzr hi(i-1,j) = xr * yi - xi * yr + hi(i,i-1) * zzi hr(i,j) = xr * zzr - xi * zzi - hi(i,i-1) * yr hi(i,j) = xr * zzi + xi * zzr - hi(i,i-1) * yi end do end do si = hi(en,en) if ( si /= 0.0D+00 ) then norm = pythag ( hr(en,en), si ) sr = hr(en,en) / norm si = si / norm hr(en,en) = norm hi(en,en) = 0.0D+00 do j = en+1, n yr = hr(en,j) yi = hi(en,j) hr(en,j) = sr * yr + si * yi hi(en,j) = sr * yi - si * yr end do end if ! ! Inverse operation (columns). ! do j = l+1, en xr = wr(j-1) xi = wi(j-1) do i = 1, j yr = hr(i,j-1) yi = 0.0D+00 zzr = hr(i,j) zzi = hi(i,j) if ( i /= j ) then yi = hi(i,j-1) hi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi end if hr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr hr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr hi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi end do do i = low, igh yr = zr(i,j-1) yi = zi(i,j-1) zzr = zr(i,j) zzi = zi(i,j) zr(i,j-1) = xr * yr - xi * yi + hi(j,j-1) * zzr zi(i,j-1) = xr * yi + xi * yr + hi(j,j-1) * zzi zr(i,j) = xr * zzr + xi * zzi - hi(j,j-1) * yr zi(i,j) = xr * zzi - xi * zzr - hi(j,j-1) * yi end do end do if ( si /= 0.0D+00 ) then do i = 1, en yr = hr(i,en) yi = hi(i,en) hr(i,en) = sr * yr - si * yi hi(i,en) = sr * yi + si * yr end do do i = low, igh yr = zr(i,en) yi = zi(i,en) zr(i,en) = sr * yr - si * yi zi(i,en) = sr * yi + si * yr end do end if go to 240 ! ! A root found. ! 660 continue hr(en,en) = hr(en,en) + tr wr(en) = hr(en,en) hi(en,en) = hi(en,en) + ti wi(en) = hi(en,en) en = enm1 go to 220 ! ! All roots found. ! Backsubstitute to find vectors of upper triangular form. ! 680 continue norm = 0.0D+00 do i = 1, n do j = i, n tr = abs ( hr(i,j) ) + abs ( hi(i,j) ) norm = max ( norm, tr ) end do end do if ( n == 1 ) then return end if if ( norm == 0.0D+00 ) then return end if do nn = 2, n en = n + 2 - nn xr = wr(en) xi = wi(en) hr(en,en) = 1.0D+00 hi(en,en) = 0.0D+00 enm1 = en - 1 do ii = 1, enm1 i = en - ii zzr = 0.0D+00 zzi = 0.0D+00 do j = i+1, en zzr = zzr + hr(i,j) * hr(j,en) - hi(i,j) * hi(j,en) zzi = zzi + hr(i,j) * hi(j,en) + hi(i,j) * hr(j,en) end do yr = xr - wr(i) yi = xi - wi(i) if ( yr == 0.0D+00 .and. yi == 0.0D+00 ) then tst1 = norm yr = tst1 do yr = 0.01D+00 * yr tst2 = norm + yr if ( tst2 <= tst1 ) then exit end if end do end if call cdiv ( zzr, zzi, yr, yi, hr(i,en), hi(i,en) ) ! ! Overflow control. ! tr = abs ( hr(i,en) ) + abs ( hi(i,en) ) if ( tr /= 0.0D+00 ) then tst1 = tr tst2 = tst1 + 1.0D+00 / tst1 if ( tst2 <= tst1 ) then do j = i, en hr(j,en) = hr(j,en)/tr hi(j,en) = hi(j,en)/tr end do end if end if end do end do ! ! End backsubstitution. ! enm1 = n - 1 ! ! Vectors of isolated roots. ! do i = 1, n - 1 if ( i < low .or. i > igh ) then do j = i+1, n zr(i,j) = hr(i,j) zi(i,j) = hi(i,j) end do end if end do ! ! Multiply by transformation matrix to give vectors of original full matrix. ! do jj = low, n - 1 j = n + low - jj m = min ( j, igh ) do i = low, igh zzr = 0.0D+00 zzi = 0.0D+00 do k = low, m zzr = zzr + zr(i,k) * hr(k,j) - zi(i,k) * hi(k,j) zzi = zzi + zr(i,k) * hi(k,j) + zi(i,k) * hr(k,j) end do zr(i,j) = zzr zi(i,j) = zzi end do end do return ! ! Set error: all eigenvalues have not converged after 30*n iterations. ! 1000 continue ierr = en return end subroutine cortb ( n, low, igh, ar, ai, ortr, orti, m, zr, zi ) !*****************************************************************************80 ! !! CORTB determines eigenvectors by undoing the CORTH transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a complex general ! matrix by back transforming those of the corresponding ! upper Hessenberg matrix determined by CORTH. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH to the order of the matrix. ! ! Input, real ( kind = 8 ) AR(N,IGH), AI(N,IGH), information about the unitary ! transformations used in the reduction by CORTH in their strict lower ! triangles. ! ! Input/output, real ( kind = 8 ) ORTR(IGH), ORTI(IGH). On input, further information ! about the transformations used in the reduction by CORTH. On output, ! ORTR and ORTI have been further altered. ! ! Input, integer ( kind = 4 ) M, the number of columns of ZR and ZI to be back ! transformed. ! ! Input/output, real ( kind = 8 ) ZR(N,M), ZI(N,M). On input, the real and imaginary ! parts of the eigenvectors to be back transformed. On output, the real ! and imaginary parts of the transformed eigenvectors. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) ai(n,igh) real ( kind = 8 ) ar(n,igh) real ( kind = 8 ) gi real ( kind = 8 ) gr real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) la integer ( kind = 4 ) low integer ( kind = 4 ) mm integer ( kind = 4 ) mp real ( kind = 8 ) orti(igh) real ( kind = 8 ) ortr(igh) real ( kind = 8 ) zi(n,m) real ( kind = 8 ) zr(n,m) if ( m == 0 ) then return end if la = igh - 1 if ( igh - 1 < low + 1 ) then return end if do mm = low + 1, la mp = low + igh - mm if ( ar(mp,mp-1) /= 0.0D+00 .or. ai(mp,mp-1) /= 0.0D+00 ) then h = ar(mp,mp-1) * ortr(mp) + ai(mp,mp-1) * orti(mp) ortr(mp+1:igh) = ar(mp+1:igh,mp-1) orti(mp+1:igh) = ai(mp+1:igh,mp-1) do j = 1, m gr = ( dot_product ( ortr(mp:igh), zr(mp:igh,j) ) & + dot_product ( orti(mp:igh), zi(mp:igh,j) ) ) / h gi = ( dot_product ( ortr(mp:igh), zi(mp:igh,j) ) & - dot_product ( orti(mp:igh), zr(mp:igh,j) ) ) / h do i = mp, igh zr(i,j) = zr(i,j) + gr * ortr(i) - gi * orti(i) zi(i,j) = zi(i,j) + gr * orti(i) + gi * ortr(i) end do end do end if end do return end subroutine corth ( n, low, igh, ar, ai, ortr, orti ) !*****************************************************************************80 ! !! CORTH transforms a complex general matrix to upper Hessenberg form. ! ! Discussion: ! ! Given a complex general matrix, this subroutine ! reduces a submatrix situated in rows and columns ! LOW through IGH to upper Hessenberg form by ! unitary similarity transformations. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine CBAL. ! If CBAL is not used, set LOW = 1 and IGH = N. ! ! Input/output, real ( kind = 8 ) AR(N,N), AI(N,N). On input, the real and imaginary ! parts of the complex input matrix. On output, the real and imaginary ! parts of the Hessenberg matrix. Information about the unitary ! transformations used in the reduction is stored in the remaining ! triangles under the Hessenberg matrix. ! ! Output, real ( kind = 8 ) ORTR(IGH), ORTI(IGH), further information about the ! transformations. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) n real ( kind = 8 ) ai(n,n) real ( kind = 8 ) ar(n,n) real ( kind = 8 ) f real ( kind = 8 ) fi real ( kind = 8 ) fr real ( kind = 8 ) g real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) ii integer ( kind = 4 ) j integer ( kind = 4 ) jj integer ( kind = 4 ) la integer ( kind = 4 ) m,mp,low real ( kind = 8 ) orti(igh) real ( kind = 8 ) ortr(igh) real ( kind = 8 ) pythag real ( kind = 8 ) scale la = igh - 1 if ( igh - 1 < low + 1 ) then return end if do m = low + 1, la h = 0.0D+00 ortr(m) = 0.0D+00 orti(m) = 0.0D+00 scale = 0.0D+00 ! ! Scale column. ! do i = m, igh scale = scale + abs ( ar(i,m-1) ) + abs ( ai(i,m-1) ) end do if ( scale == 0.0D+00 ) then cycle end if mp = m + igh do ii = m, igh i = mp - ii ortr(i) = ar(i,m-1) / scale orti(i) = ai(i,m-1) / scale h = h + ortr(i) * ortr(i) + orti(i) * orti(i) end do g = sqrt ( h ) f = pythag ( ortr(m), orti(m) ) if ( f /= 0.0D+00 ) then h = h + f * g g = g / f ortr(m) = ( 1.0D+00 + g ) * ortr(m) orti(m) = ( 1.0D+00 + g ) * orti(m) else ortr(m) = g ar(m,m-1) = scale end if ! ! Form (I-(U*Ut)/h) * A. ! do j = m, n fr = 0.0D+00 fi = 0.0D+00 do ii = m, igh i = mp - ii fr = fr + ortr(i) * ar(i,j) + orti(i) * ai(i,j) fi = fi + ortr(i) * ai(i,j) - orti(i) * ar(i,j) end do fr = fr / h fi = fi / h ar(m:igh,j) = ar(m:igh,j) - fr * ortr(m:igh) + fi * orti(m:igh) ai(m:igh,j) = ai(m:igh,j) - fr * orti(m:igh) - fi * ortr(m:igh) end do ! ! Form (I-(U*Ut)/h) * A * (I-(U*Ut)/h) ! do i = 1, igh fr = 0.0D+00 fi = 0.0D+00 do jj = m, igh j = mp - jj fr = fr + ortr(j) * ar(i,j) - orti(j) * ai(i,j) fi = fi + ortr(j) * ai(i,j) + orti(j) * ar(i,j) end do fr = fr / h fi = fi / h ar(i,m:igh) = ar(i,m:igh) - fr * ortr(m:igh) - fi * orti(m:igh) ai(i,m:igh) = ai(i,m:igh) + fr * orti(m:igh) - fi * ortr(m:igh) end do ortr(m) = scale * ortr(m) orti(m) = scale * orti(m) ar(m,m-1) = - g * ar(m,m-1) ai(m,m-1) = - g * ai(m,m-1) end do return end subroutine csroot ( xr, xi, yr, yi ) !*****************************************************************************80 ! !! CSROOT computes the complex square root of a complex quantity. ! ! Discussion: ! ! The branch of the square function is chosen so that ! YR >= 0.0D+00 ! and ! sign ( YI ) == sign ( XI ) ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, real ( kind = 8 ) XR, XI, the real and imaginary parts of the quantity ! whose square root is desired. ! ! Output, real ( kind = 8 ) YR, YI, the real and imaginary parts of the square root. ! implicit none real ( kind = 8 ) pythag real ( kind = 8 ) s real ( kind = 8 ) ti real ( kind = 8 ) tr real ( kind = 8 ) xi real ( kind = 8 ) xr real ( kind = 8 ) yi real ( kind = 8 ) yr tr = xr ti = xi s = sqrt ( 0.5D+00 * ( pythag ( tr, ti ) + abs ( tr ) ) ) if ( tr >= 0.0D+00 ) yr = s if ( ti < 0.0D+00 ) s = -s if ( tr <= 0.0D+00 ) yi = s if ( tr < 0.0D+00 ) then yr = 0.5D+00 * ( ti / yi ) else if ( tr > 0.0D+00 ) then yi = 0.5D+00 * ( ti / yr ) end if return end subroutine elmbak ( n, low, igh, a, ind, m, z ) !*****************************************************************************80 ! !! ELMBAK determines eigenvectors by undoing the ELMHES transformation. ! ! Discussion: ! ! This subroutine forms the eigenvectors of a real general ! matrix by back transforming those of the corresponding ! upper Hessenberg matrix determined by ELMHES. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, integers determined by the balancing ! routine BALANC. If BALANC has not been used, set LOW = 1 and ! IGH equal to the order of the matrix. ! ! Input, real ( kind = 8 ) A(N,IGH), the multipliers which were used in the ! reduction by ELMHES in its lower triangle below the subdiagonal. ! ! Input, integer ( kind = 4 ) IND(IGH), information on the rows and columns ! interchanged in the reduction by ELMHES. ! ! Input, integer ( kind = 4 ) M, the number of columns of Z to be back transformed. ! ! Input/output, real ( kind = 8 ) Z(N,M). On input, the real and imaginary parts ! of the eigenvectors to be back transformed. On output, the real and ! imaginary parts of the transformed eigenvectors. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) a(n,igh) integer ( kind = 4 ) i integer ( kind = 4 ) ind(igh) integer ( kind = 4 ) j integer ( kind = 4 ) la integer ( kind = 4 ) low integer ( kind = 4 ) mm integer ( kind = 4 ) mp real ( kind = 8 ) x real ( kind = 8 ) z(n,m) if ( m == 0 ) then return end if la = igh - 1 if ( la < low + 1 ) then return end if do mm = low + 1, la mp = low + igh - mm do i = mp+1, igh x = a(i,mp-1) if ( x /= 0.0D+00 ) then do j = 1, m z(i,j) = z(i,j) + x * z(mp,j) end do end if end do i = ind(mp) if ( i /= mp ) then do j = 1, m call r8_swap ( z(i,j), z(mp,j) ) end do end if end do return end subroutine elmhes ( n, low, igh, a, ind ) !*****************************************************************************80 ! !! ELMHES transforms a real general matrix to upper Hessenberg form. ! ! Discussion: ! ! Given a real general matrix, this subroutine reduces a submatrix ! situated in rows and columns LOW through IGH to upper Hessenberg ! form by stabilized elementary similarity transformations. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! Martin and Wilkinson, ! ELMHES, ! Numerische Mathematik, ! Volume 12, pages 349-368, 1968. ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine ! BALANC. If BALANC has not been used, set LOW = 1, IGH = N. ! ! Input/output, real ( kind = 8 ) A(N,N). On input, the matrix to be reduced. ! On output, the Hessenberg matrix. The multipliers ! which were used in the reduction are stored in the ! remaining triangle under the Hessenberg matrix. ! ! Output, integer ( kind = 4 ) IND(N), contains information on the rows and columns ! interchanged in the reduction. Only elements LOW through IGH are used. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) n real ( kind = 8 ) a(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) ind(igh) integer ( kind = 4 ) j integer ( kind = 4 ) la integer ( kind = 4 ) low integer ( kind = 4 ) m real ( kind = 8 ) x real ( kind = 8 ) y la = igh - 1 do m = low + 1, la x = 0.0D+00 i = m do j = m, igh if ( abs ( a(j,m-1) ) > abs ( x ) ) then x = a(j,m-1) i = j end if end do ind(m) = i ! ! Interchange rows and columns of the matrix. ! if ( i /= m ) then do j = m-1, n call r8_swap ( a(i,j), a(m,j) ) end do do j = 1, igh call r8_swap ( a(j,i), a(j,m) ) end do end if if ( x /= 0.0D+00 ) then do i = m+1, igh y = a(i,m-1) if ( y /= 0.0D+00 ) then y = y / x a(i,m-1) = y do j = m, n a(i,j) = a(i,j) - y * a(m,j) end do a(1:igh,m) = a(1:igh,m) + y * a(1:igh,i) end if end do end if end do return end subroutine eltran ( n, low, igh, a, ind, z ) !*****************************************************************************80 ! !! ELTRAN accumulates similarity transformations used by ELMHES. ! ! Discussion: ! ! This subroutine accumulates the stabilized elementary ! similarity transformations used in the reduction of a ! real general matrix to upper Hessenberg form by ELMHES. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! Peters and WIlkinson, ! ELMTRANS, ! Numerische Mathematik, ! Volume 16, pages 181-204, 1970. ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, are determined by the balancing routine ! BALANC. If BALANC has not been used, set LOW = 1, IGH = N. ! ! Input, real ( kind = 8 ) A(N,IGH), the multipliers which were used in the ! reduction by ELMHES in its lower triangle below the subdiagonal. ! ! Input, integer ( kind = 4 ) IND(IGH), information on the rows and columns ! interchanged in the reduction by ELMHES. ! ! Output, real ( kind = 8 ) Z(N,N), the transformation matrix produced in the ! reduction by ELMHES. ! implicit none integer ( kind = 4 ) igh integer ( kind = 4 ) n real ( kind = 8 ) a(n,igh) integer ( kind = 4 ) i integer ( kind = 4 ) ind(igh) integer ( kind = 4 ) j integer ( kind = 4 ) kl integer ( kind = 4 ) low integer ( kind = 4 ) mm integer ( kind = 4 ) mp real ( kind = 8 ) z(n,n) ! ! Initialize Z to the identity matrix. ! call r8mat_identity ( n, z ) kl = igh - low - 1 if ( kl < 1 ) then return end if do mm = 1, kl mp = igh - mm do i = mp+1, igh z(i,mp) = a(i,mp-1) end do i = ind(mp) if ( i /= mp ) then z(mp,mp:igh) = z(i,mp:igh) z(i,mp) = 1.0D+00 z(i,mp+1:igh) = 0.0D+00 end if end do return end subroutine figi ( n, t, d, e, e2, ierr ) !*****************************************************************************80 ! !! FIGI transforms a real nonsymmetric tridiagonal matrix to symmetric form. ! ! Discussion: ! ! Given a nonsymmetric tridiagonal matrix such that the products ! of corresponding pairs of off-diagonal elements are all ! non-negative, this subroutine reduces it to a symmetric ! tridiagonal matrix with the same eigenvalues. If, further, ! a zero product only occurs when both factors are zero, ! the reduced matrix is similar to the original matrix. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, real ( kind = 8 ) T(N,3) contains the input matrix. Its subdiagonal is ! stored in the last N-1 positions of the first column, its diagonal in ! the N positions of the second column, and its superdiagonal in the ! first N-1 positions of the third column. T(1,1) and T(N,3) are arbitrary. ! ! Output, real ( kind = 8 ) D(N), the diagonal elements of the symmetric matrix. ! ! Output, real ( kind = 8 ) E(N), contains the subdiagonal elements of the symmetric ! matrix in E(2:N). E(1) is not set. ! ! Output, real ( kind = 8 ) E2(N), the squares of the corresponding elements of E. ! E2 may coincide with E if the squares are not needed. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! N+I, if T(I,1) * T(I-1,3) is negative, ! -(3*N+I), if T(I,1) * T(I-1,3) is zero with one factor non-zero. In ! this case, the eigenvectors of the symmetric matrix are not simply ! related to those of T and should not be sought. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) d(n) real ( kind = 8 ) e(n) real ( kind = 8 ) e2(n) integer ( kind = 4 ) i integer ( kind = 4 ) ierr real ( kind = 8 ) t(n,3) ierr = 0 do i = 1, n if ( i >= 1 ) then e2(i) = t(i,1) * t(i-1,3) if ( e2(i) < 0.0D+00 ) then ierr = n + i return else if ( e2(i) == 0.0D+00 ) then if ( t(i,1) /= 0.0D+00 .or. t(i-1,3) /= 0.0D+00 ) then ierr = - 3 * n - i return end if e(i) = 0.0D+00 else e(i) = sqrt ( e2(i) ) end if end if d(i) = t(i,2) end do return end subroutine figi2 ( n, t, d, e, z, ierr ) !*****************************************************************************80 ! !! FIGI2 transforms a real nonsymmetric tridiagonal matrix to symmetric form. ! ! Discussion: ! ! Given a nonsymmetric tridiagonal matrix such that the products ! of corresponding pairs of off-diagonal elements are all ! non-negative, and zero only when both factors are zero, this ! subroutine reduces it to a symmetric tridiagonal matrix ! using and accumulating diagonal similarity transformations. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, real ( kind = 8 ) T(N,3) contains the input matrix. Its subdiagonal is ! stored in the last N-1 positions of the first column, its diagonal in ! the N positions of the second column, and its superdiagonal in the ! first N-1 positions of the third column. T(1,1) and T(N,3) are arbitrary. ! ! Output, real ( kind = 8 ) D(N), the diagonal elements of the symmetric matrix. ! ! Output, real ( kind = 8 ) E(N), contains the subdiagonal elements of the symmetric ! matrix in E(2:N). E(1) is not set. ! ! Output, real ( kind = 8 ) Z(N,N), contains the transformation matrix produced in ! the reduction. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, for normal return, ! N+I, if T(I,1) * T(I-1,3) is negative, ! 2*N+I, if T(I,1) * T(I-1,3) is zero with one factor non-zero. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) d(n) real ( kind = 8 ) e(n) real ( kind = 8 ) h integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) j real ( kind = 8 ) t(n,3) real ( kind = 8 ) z(n,n) ierr = 0 do i = 1, n z(i,1:n) = 0.0D+00 if ( i == 1 ) then z(i,i) = 1.0D+00 else h = t(i,1) * t(i-1,3) if ( h < 0.0D+00 ) then ierr = n + i return else if ( h == 0 ) then if ( t(i,1) /= 0.0D+00 .or. t(i-1,3) /= 0.0D+00 ) then ierr = 2 * n + i return end if e(i) = 0.0D+00 z(i,i) = 1.0D+00 else e(i) = sqrt ( h ) z(i,i) = z(i-1,i-1) * e(i) / t(i-1,3) end if end if d(i) = t(i,2) end do return end subroutine hqr ( n, low, igh, h, wr, wi, ierr ) !*****************************************************************************80 ! !! HQR computes all eigenvalues of a real upper Hessenberg matrix. ! ! Discussion: ! ! This subroutine finds the eigenvalues of a real ! upper Hessenberg matrix by the QR method. ! ! Modified: ! ! 04 February 2003 ! ! Reference: ! ! Martin, Peters, and Wilkinson, ! HQR, ! Numerische Mathematik, ! Volume 14, pages 219-231, 1970. ! ! James Wilkinson, Christian Reinsch, ! Handbook for Automatic Computation, ! Volume II, Linear Algebra, Part 2, ! Springer, 1971, ! ISBN: 0387054146, ! LC: QA251.W67. ! ! B Smith, J Boyle, J Dongarra, B Garbow, Y Ikebe, V Klema, C Moler, ! Matrix Eigensystem Routines, EISPACK Guide, ! Lecture Notes in Computer Science, Volume 6, ! Springer Verlag, 1976. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order of the matrix. ! ! Input, integer ( kind = 4 ) LOW, IGH, two integers determined by the routine ! BALANC. If BALANC is not used, set LOW=1, IGH=N. ! ! Input/output, real ( kind = 8 ) H(N,N), the N by N upper Hessenberg matrix. ! Information about the transformations used in the reduction to ! Hessenberg form by ELMHES or ORTHES, if performed, is stored ! in the remaining triangle under the Hessenberg matrix. ! On output, the information in H has been destroyed. ! ! Output, real ( kind = 8 ) WR(N), WI(N), the real and imaginary parts of the ! eigenvalues. The eigenvalues are unordered, except that complex ! conjugate pairs of values appear consecutively, with the eigenvalue ! having positive imaginary part listed first. If an error exit ! occurred, then the eigenvalues should be correct for indices ! IERR+1 through N. ! ! Output, integer ( kind = 4 ) IERR, error flag. ! 0, no error. ! J, the limit of 30*N iterations was reached while searching for ! the J-th eigenvalue. ! implicit none integer ( kind = 4 ) n integer ( kind = 4 ) en integer ( kind = 4 ) enm2 real ( kind = 8 ) h(n,n) integer ( kind = 4 ) i integer ( kind = 4 ) ierr integer ( kind = 4 ) igh integer ( kind = 4 ) itn integer ( kind = 4 ) its integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) l integer ( kind = 4 ) ll integer ( kind = 4 ) low integer ( kind = 4 ) m integer ( kind = 4 ) mm integer ( kind = 4 ) na real ( kind = 8 ) norm logical notlas real ( kind = 8 ) p real ( kind = 8 ) q real ( kind = 8 ) r real ( kind = 8 ) s real ( kind = 8 ) t real ( kind = 8 ) tst1 real ( kind = 8 ) tst2 real ( kind = 8 ) w real ( kind = 8 ) wi(n) real ( kind = 8 ) wr(n) real ( kind = 8 ) x real ( kind = 8 ) y real ( kind = 8 ) zz ierr = 0 norm = 0.0D+00 k = 1 ! ! Store roots isolated by BALANC and compute matrix norm. ! do i = 1, n do j = k, n norm = norm + abs ( h(i,j) ) end do k = i if ( i < low .or. i > igh ) then wr(i) = h(i,i) wi(i) = 0.0D+00 end if end do en = igh t = 0.0D+00 itn = 30 * n ! ! Search for next eigenvalues. ! 60 continue if ( en < low ) then return end if its = 0 na = en - 1 enm2 = na - 1 ! ! Look for a single small sub-diagonal element. ! 70 continue do ll = low, en l = en + low - ll if ( l == low ) then exit end if s = abs ( h(l-1,l-1) ) + abs ( h(l,l) ) if ( s == 0.0D+00 ) then s = norm end if tst1 = s tst2 = tst1 + abs ( h(l,l-1) ) if ( tst2 == tst1 ) then exit end if end do ! ! Form shift. ! x = h(en,en) if ( l == en ) then go to 270 end if y = h(na,na) w = h(en,na) * h(na,en) if ( l == na ) then go