MODULE skew_symmetric IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) CONTAINS SUBROUTINE trizd(n, a, e) ! **** ! FUNCTION - REDUCES A REAL SKEW-SYMMETRIC MATRIX TO A SKEW-SYMMETRIC ! TRIDIAGONAL MATRIX USING ORTHOGONAL SIMILARITY ! TRANSFORMATIONS ! PARAMETERS ! N - INPUT INTEGER SPECIFYING THE ORDER OF A ! A(NA,N) - ON INPUT, A CONTAINS THE REAL SKEW-SYMMETRIC MATRIX. ! ONLY THE STRICT LOWER TRIANGLE OF THE MATRIX NEED ! BE SUPPLIED. ! ON OUTPUT, A CONTAINS INFORMATION ABOUT THE ORTHOGONAL ! TRANSFORMATIONS USED IN THE REDUCTION IN ITS FULL LOWER ! TRIANGLE. THE STRICT UPPER TRIANGLE OF A IS UNALTERED. ! E(N) - OUTPUT ARRAY CONTAINING THE LOWER SUBDIAGONAL ELEMENTS ! OF THE TRIDIAGONAL MATRIX IN ITS LAST N-1 POSITIONS. ! E(1) IS SET TO ZERO. ! REQUIRED FUNCTIONS - ABS,SIGN,SQRT ! **** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: a(:,:) REAL (dp), INTENT(OUT) :: e(:) REAL (dp) :: f, g, h, scale INTEGER :: i, j, k, l, ii IF (n /= 1) THEN ! *** MAIN DO LOOP I=N STEP -1 UNTIL 2 DO ii = 2, n i = n + 2 - ii l = i - 1 h = 0. scale = SUM( ABS(a(i,1:l)) ) ! *** NORMALIZE ROW IF (scale == 0.) THEN e(i) = 0. ELSE ! *** COMPUTE ELEMENTS OF U VECTOR DO k = 1, l a(i,k) = a(i,k) / scale h = h + a(i,k) * a(i,k) END DO f = a(i,l) g = -SIGN(SQRT(h),f) e(i) = scale * g h = h - f * g a(i,l) = f - g IF (l /= 1) THEN ! *** COMPUTE ELEMENTS OF A*U/H DO j = 1, l g = DOT_PRODUCT( a(j,1:j-1), a(i,1:j-1) ) - & DOT_PRODUCT( a(j+1:l,j), a(i,j+1:l) ) e(j) = g / h END DO ! *** COMPUTE REDUCED A DO j = 2, l f = a(i,j) g = e(j) DO k = 1, j-1 a(j,k) = a(j,k) + f * e(k) - g * a(i,k) END DO END DO END IF a(i,1:l) = scale * a(i,1:l) END IF a(i,i) = scale * SQRT(h) END DO END IF e(1) = 0. RETURN END SUBROUTINE trizd SUBROUTINE imzd(n, e, matz, skew, z, ierr) ! **** ! FUNCTION - COMPUTE THE EIGENVALUES AND OPTIONALLY THE EIGENVECTORS OF A ! SYMMETRIC TRIDIAGONAL MATRIX WITH ZERO DIAGONALS OR A SKEW- ! SYMMETRIC TRIDIAGONAL MATRIX USING AN IMPLICIT QR-TYPE ! ITERATION ! PARAMETERS ! N - INPUT INTEGER SPECIFYING THE ORDER OF THE TRIDIAGONAL MATRIX ! E(N) - ON INPUT, ARRAY CONTAINING THE LOWER SUBDIAGONAL ELEMENTS OF ! THE TRIDIAGONAL MATRIX IN ITS LAST N-1 POSITIONS. ! E(1) IS ARBITRARY. ! ON OUTPUT, ARRAY CONTAINS THE EIGENVALUES. THE NON-ZERO ! EIGENVALUES OCCUR IN PAIRS WITH OPPOSITE SIGNS AND ARE FOUND ! IN ADJACENT LOCATIONS IN E. THE EIGENVALUES OF SYMMETRIC ! MATRICES ARE REAL AND THE EIGENVALUES OF SKEW-SYMMETRIC ! MATRICES ARE PURELY IMAGINARY COMPLEX NUMBERS. IF AN ERROR ! EXIT IS MADE, THE EIGENVALUES ARE CORRECT FOR INDICES ! IERR+1, IERR+2...N ! MATZ - INPUT LOGICAL VARIABLE SPECIFYING THE EIGENVECTOR ! OPTION ! = .TRUE. EIGENVECTORS ARE TO BE COMPUTED ! = .FALSE. EIGENVECTORS ARE NOT TO BE COMPUTED ! SKEW - INPUT LOGICAL VARIABLE SPECIFYING TYPE OF INPUT MATRIX ! = .TRUE. INPUT TRIDIAGONAL MATRIX IS SKEW-SYMMETRIC ! = .FALSE. INPUT TRIDIAGONAL MATRIX IS SYMMETRIC WITH ZERO ! DIAGONALS ! SKEW IS NOT REFERENCED IF MATZ = .FALSE. ! Z(NZ,N) - OUTPUT ARRAY CONTAINING THE ORTHOGONAL EIGENVECTORS OF THE ! INPUT TRIDIAGONAL MATRIX. EIGENVECTORS CORRESPONDING TO ZERO ! EIGENVALUES ARE NORMALIZE TO UNIT 2-NORM (LENGTH) AND THOSE ! CORRESPONDING TO NON-ZERO EIGENVALUES HAVE 2-NORM OF SQUARE ! ROOT 2. IF THE J-TH EIGENVALUE IS ZERO OR REAL (I.E. E(J)), ! ITS EIGENVECTOR IS FOUND IN THE J-TH COLUMN OF Z. ! IF THE J-TH EIGENVALUE IS IMAGINARY (I.E. E(J)*I) WITH ! E(J+1) = -E(J), THE REAL PART OF ITS EIGENVECTOR IS FOUND IN ! THE J-TH COLUMN OF Z AND ITS IMAGINARY PART FOUND IN THE ! (J+1)-TH COLUMN. IF AN ERROR EXIT IS MADE, Z CONTAINS THE ! EIGENVECTORS ASSOCIATED WITH THE STORED EIGENVALUES. ! Z IS NOT REFERENCED IF MATZ = .FALSE. ! IERR - OUTPUT ERROR CODE ! = 0 NORMAL RETURN (ALL EIGENVALUES/VECTORS FOUND) ! = J IF THE J-TH EIGENVALUE HAS NOT BEEN DETERMINED AFTER ! 30 ITERATIONS ! REQUIRED FUNCTIONS - ABS,SIGN,SQRT,MOD ! **** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: e(:) LOGICAL, INTENT(IN) :: matz LOGICAL, INTENT(IN) :: skew REAL (dp), INTENT(OUT) :: z(:,:) INTEGER, INTENT(OUT) :: ierr REAL (dp) :: f, g, q, c, s, r, p, test, tmag INTEGER :: i, j, k, l, m, l0, l0m1, m0, ls, im1, jp1, km1, kp1, lm1, & mm1, ip3, ieo, its, ip1 REAL (dp), PARAMETER :: eps = EPSILON( 0.0_dp ) IF (matz) THEN ! *** PLACE IDENTITY MATRIX IN Z DO i = 1, n z(i,1:n) = 0. z(i,i) = 1. END DO END IF ierr = 0 m = n mm1 = m - 1 e(1) = 0. its = 0 30 IF (m >= 2) THEN m0 = m ! *** SEARCH FOR NEXT SUBMATRIX TO SOLVE (MATRIX SPLITTING) f = 0. DO i = 1, mm1 j = m - i jp1 = j + 1 g = ABS(e(jp1)) tmag = ABS(e(j)) + f test = tmag + g IF (ABS(test - tmag) < eps*tmag) GO TO 50 f = g END DO jp1 = 1 50 l0 = jp1 + 1 l0m1 = jp1 IF (l0m1 == m) GO TO 160 IF (matz) THEN IF (skew) THEN ! *** PLACE CORRECT SIGN ON IDENTITY DIAGONALS DO i = l0m1, m, 4 z(i,i) = -z(i,i) ip3 = i + 3 IF (ip3 > m) GO TO 70 z(ip3,ip3) = -z(ip3,ip3) END DO END IF END IF 70 IF (l0 == m) GO TO 170 ieo = m - l0 ieo = MOD(ieo,2) l = l0 IF (ieo == 0) GO TO 130 ! *** FIND ZERO EIGENVALUE OF ODD ORDERED SUBMATRICES c = 0. s = -1. DO i = l0, mm1, 2 k = mm1 + l0 - i kp1 = k + 1 q = -s * e(kp1) e(kp1) = c * e(kp1) IF (ABS(e(k)) <= ABS(q)) THEN c = e(k) / q r = SQRT(c*c+1.) e(k) = q * r s = 1. / r c = c * s ELSE s = q / e(k) r = SQRT(1.+s*s) e(k) = e(k) * r c = 1. / r s = s * c END IF IF (matz) THEN ! *** ACCUMULATE TRANSFORMATIONS FOR EIGENVECTORS km1 = k - 1 z(km1,m) = -s * z(km1,km1) z(km1,km1) = c * z(km1,km1) DO j = kp1, m, 2 z(j,km1) = s * z(j,m) z(j,m) = c * z(j,m) END DO END IF END DO m = mm1 mm1 = m - 1 IF (l0 == m) GO TO 170 ! *** CHECK FOR CONVERGENCE OR SMALL SUBDIAGONAL ELEMENT 100 DO i = l0, mm1, 2 k = mm1 + l0 - i l = k + 1 tmag = ABS(e(l)) + ABS(e(k-1)) test = tmag + e(k) IF (ABS(test - tmag) < eps*tmag) GO TO 120 END DO l = l0 120 IF (l == m) GO TO 170 ! *** FORM SHIFT 130 its = its + 1 IF (its > 30) GO TO 220 f = e(m-3) g = e(m-2) c = e(mm1) s = e(m) p = ((c-f)*(c+f)+(s-g)*(s+g)) / (2.*g*c) r = SQRT(p*p+1.) q = (g/(p+SIGN(r,p))) - c f = e(l) lm1 = l - 1 e(lm1) = ((f-s)*(f+s)+c*q) / f ! *** PERFORM ONE IMPLICIT QR ITERATION ON CHOLESKY FACTOR ls = l0m1 c = 1. s = 1. DO i = l, mm1 ip1 = i + 1 im1 = i - 1 q = s * e(ip1) e(ip1) = c * e(ip1) IF (ABS(e(im1)) <= ABS(q)) THEN c = e(im1) / q r = SQRT(c*c+1.) e(im1) = q * r s = 1. / r c = c * s ELSE s = q / e(im1) r = SQRT(1.+s*s) e(im1) = e(im1) * r c = 1. / r s = s * c END IF f = e(ip1) e(ip1) = -s * e(i) + c * f e(i) = c * e(i) + s * f IF (matz) THEN ! *** ACCUMULATE TRANSFORMATIONS FOR EIGENVECTORS DO j = ls, m0, 2 f = z(j,ip1) z(j,ip1) = -s * z(j,im1) + c * f z(j,im1) = c * z(j,im1) + s * f END DO IF (ls /= l0m1) THEN ls = l0m1 ELSE ls = l0 END IF END IF END DO e(lm1) = 0. GO TO 100 ! *** ITERATION CONVERGED TO ONE ZERO EIGENVALUE 160 e(m) = 0. m = mm1 GO TO 180 ! *** ITERATION CONVERGED TO EIGENVALUE PAIR 170 e(mm1) = e(m) e(m) = -e(m) m = m - 2 180 its = 0 mm1 = m - 1 IF (m > l0) GO TO 100 IF (m == l0) GO TO 170 IF (.NOT.matz) GO TO 30 IF (skew) GO TO 30 ! *** COMPUTE EIGENVECTORS FROM ORTHONORMAL COLUMNS OF Z IF NOT SKEW 190 k = m0 200 IF (e(k) /= 0.) THEN km1 = k - 1 DO j = l0m1, m0, 2 z(j,k) = z(j,km1) f = z(j+1,k) z(j+1,km1) = f z(j+1,k) = -f END DO k = km1 END IF k = k - 1 IF (k > l0m1) GO TO 200 IF (ierr /= 0) GO TO 230 GO TO 30 ! *** ERROR EXIT 220 ierr = m IF (matz) THEN IF (.NOT.skew) GO TO 190 END IF END IF 230 RETURN END SUBROUTINE imzd SUBROUTINE tbakzd(n, a, m, z) ! **** ! FUNCTION - FORMS THE EIGENVECTORS OF A REAL SKEW-SYMMETRIC MATRIX ! BY BACK TRANSFORMING THOSE OF THE CORRESPONDING SKEW- ! SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRIZD ! PARAMETERS ! N - INPUT INTEGER SPECIFYING THE ORDER OF A ! A(NA,N) - INPUT ARRAY CONTAINING INFORMATION ABOUT THE ORTHOGONAL ! TRANSFORMATIONS USED IN THE REDUCTION BY TRIZD IN ! ITS FULL LOWER TRIANGLE ! M - INPUT INTEGER SPECIFYING THE NUMBER OF EIGENVECTORS TO ! BE BACK TRANSFORMED ! Z(NZ,M) - ON INPUT, Z CONTAINS THE REAL AND IMAGINARY (IF COMPLEX) ! PARTS OF THE EIGENVECTORS TO BE BACK TRANSFORMED IN ITS ! FIRST M COLUMNS ! ON OUTPUT, Z CONTAINS THE REAL AND IMAGINARY (IF COMPLEX) ! PARTS OF THE TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS ! **** INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: a(:,:) INTEGER, INTENT(IN) :: m REAL (dp), INTENT(IN OUT) :: z(:,:) REAL (dp) :: h, s INTEGER :: i, j, l IF (m /= 0) THEN IF (n /= 1) THEN DO i = 2, n l = i - 1 h = a(i,i) IF (h /= 0.0_dp) THEN DO j = 1, m s = DOT_PRODUCT( a(i,1:l), z(1:l,j) ) s = (s/h) / h z(1:l,j) = z(1:l,j) - s * a(i,1:l) END DO END IF END DO END IF END IF RETURN END SUBROUTINE tbakzd END MODULE skew_symmetric PROGRAM driver ! Code converted using TO_F90 by Alan Miller ! Date: 1999-10-29 Time: 18:49:08 ! **** ! FUNCTION - THIS IS THE MAIN PROGRAM (DRIVER) FOR ILLUSTRATING THE ! USE OF SUBROUTINES TRIZD, IMZD, AND TBAKZD. ! REFERENCES - EIGENSYSTEM COMPUTATION FOR SKEW-SYMMETRIC MATRICES AND ! A CLASS OF SYMMETRIC MATRICES, WARD,R C AND GRAY,L J. ! TO APPEAR IN MANUSCRIPT SECTION OF ACM-TOMS. ! AN ALGORITHM FOR COMPUTING THE EIGENSYSTEM OF SKEW- ! SYMMETRIC MATRICES AND A CLASS OF SYMMETRIC MATRICES, ! TO APPEAR IN ALGORITHM SECTION OF ACM-TOMS. ! REQUIRED FUNCTIONS FOR DRIVER AND SUBROUTINES - ABS,SIGN,SQRT,MOD ! **** USE skew_symmetric IMPLICIT NONE REAL (dp) :: a(6,6), z(6,6), e(6) REAL (dp) :: con INTEGER :: i, j, n, im1, jp1, ierr LOGICAL :: matz, skew matz = .true. ! *** SET UP SKEW-SYMMETRIC TEST CASE WRITE (6,5000) n = 5 skew = .true. ! *** READ AND PRINT TEST MATRIX WRITE (6,5100) DO i = 1, n READ (5,5200) a(i,1:n) WRITE (6,5200) a(i,1:n) END DO ! *** COMPUTE EIGENVALUES AND EIGENVECTORS CALL trizd(n, a, e) CALL imzd(n, e, matz, skew, z, ierr) IF (ierr /= 0) WRITE (6,5300) ierr CALL tbakzd(n, a, n, z) ! *** PRINT EIGENVALUES AND EIGENVECTORS WRITE (6,5400) j = 0 20 j = j + 1 WRITE (6,5500) IF (e(j) /= 0.d0) THEN jp1 = j + 1 WRITE (6,5600) e(j), z(1,j), z(1,jp1) WRITE (6,5700) (z(i,j), z(i,jp1), i = 2,n) WRITE (6,5500) con = -z(1,jp1) WRITE (6,5600) e(jp1), z(1,j), con DO i = 2, n con = -z(i,jp1) WRITE (6,5700) z(i,j), con END DO j = j + 1 ELSE WRITE (6,5800) e(j), z(1,j) WRITE (6,5900) z(2:n,j) END IF IF (j < n) GO TO 20 ! *** SET UP TRIDIAGONAL, SYMMETRIC, ZERO DIAGONAL TEST CASE WRITE (6,6000) n = 6 skew = .false. DO i = 1, n e(i) = 1. DO j = 1, n a(i,j) = 0. END DO END DO DO i = 2, n im1 = i - 1 a(i,im1) = e(i) a(im1,i) = e(i) END DO ! *** PRINT TEST MATRIX WRITE (6,6100) DO i = 1, n WRITE (6,5200) a(i,1:n) END DO ! *** COMPUTE EIGENVALUES AND EIGENVECTORS CALL imzd(n, e, matz, skew, z, ierr) IF (ierr /= 0) WRITE (6,5300) ierr ! *** PRINT EIGENVALUES AND EIGENVECTORS WRITE (6,6200) DO j = 1, n WRITE (6,5500) WRITE (6,5800) e(j), z(1,j) WRITE (6,5900) z(2:n,j) END DO STOP 5000 FORMAT ('1', 'EIGENSYSTEM COMPUTATION OF SKEW-SYMMETRIC TEST CASE'//) 5100 FORMAT ('0', ' TEST MATRIX'/) 5200 FORMAT (6F6.0) 5300 FORMAT ('0'/'0', 'IMZD IERR =', i5) 5400 FORMAT ('0'/'0', ' EIGENVALUES', t40, 'EIGENVECTORS') 5500 FORMAT (/) 5600 FORMAT (' ', e15.8, ' * I ', e15.8, ' + ', e15.8, ' * I') 5700 FORMAT (t26, e15.8, ' + ', e15.8, ' * I') 5800 FORMAT (' ', e15.8, t26, e15.8) 5900 FORMAT (t26, e15.8) 6000 FORMAT ('1', 'EIGENSYSTEM COMPUTATION OF TRIDIAGONAL, SYMMETRIC,', & ' ZERO DIAGONAL TEST CASE'//) 6100 FORMAT ('0', t14, 'TEST MATRIX'/) 6200 FORMAT ('0'/'0', ' EIGENVALUES EIGENVECTORS') END PROGRAM driver