C ALGORITHM 800, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 26, NO. 1, March, 2000, P. 49-- 77. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc/ # Doc/PackingList # Doc/README # Fortran77/ # Fortran77/Dp/ # Fortran77/Dp/Drivers/ # Fortran77/Dp/Drivers/data1 # Fortran77/Dp/Drivers/data2 # Fortran77/Dp/Drivers/data4 # Fortran77/Dp/Drivers/data6 # Fortran77/Dp/Drivers/data8 # Fortran77/Dp/Drivers/driver1.f # Fortran77/Dp/Drivers/driver2.f # Fortran77/Dp/Drivers/driver3.f # Fortran77/Dp/Drivers/driver4.f # Fortran77/Dp/Drivers/driver5.f # Fortran77/Dp/Drivers/driver6.f # Fortran77/Dp/Drivers/driver7.f # Fortran77/Dp/Drivers/driver8.f # Fortran77/Dp/Drivers/driver9.f # Fortran77/Dp/Drivers/res1 # Fortran77/Dp/Drivers/res2 # Fortran77/Dp/Drivers/res3 # Fortran77/Dp/Drivers/res4 # Fortran77/Dp/Drivers/res5 # Fortran77/Dp/Drivers/res6 # Fortran77/Dp/Drivers/res7 # Fortran77/Dp/Drivers/res8 # Fortran77/Dp/Drivers/res9 # Fortran77/Dp/Src/ # Fortran77/Dp/Src/src.f # This archive created: Mon Sep 4 12:17:42 2000 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'PackingList' then echo shar: will not over-write existing file "'PackingList'" else cat << SHAR_EOF > 'PackingList' Algorithm:999 Language:Fortran77 Precision:Dp Src:src.f SrcLibs:lapack blas Driver_0:driver1.f @Src Data_0:stdin=data1 Res_0:stdout=res1 Driver_1:driver2.f @Src Data_1:stdin=data2 Res_1:stdout=res2 Driver_2:driver3.f @Src Res_2:stdout=res3 Driver_3:driver4.f @Src Data_3:stdin=data4 Res_3:stdout=res4 Driver_4:driver5.f @Src Res_4:stdout=res5 Driver_5:driver6.f @Src Data_5:stdin=data6 Res_5:stdout=res6 Driver_6:driver7.f @Src Res_6:stdout=res7 Driver_7:driver8.f @Src Data_7:stdin=data8 Res_7:stdout=res8 Driver_8:driver9.f @Src Res_8:stdout=res9 SHAR_EOF fi # end of overwriting check if test -f 'README' then echo shar: will not over-write existing file "'README'" else cat << SHAR_EOF > 'README' DHABL/DHASRD/DHAEVS/DHAEVD ========================== The FORTRAN 77 subroutines DHABL, DHASRD, DHAEVS, and DHAEVD are designed to compute the eigenvalues of a Hamiltonian matrix ( A G ) (*) H = ( T ) ( Q -A ) by Van Loan's square reduced method as given in [1,2]. Here, A, G, and Q are real n-by-n matrices where G and Q are symmetric. DHABL applies one of two different kinds of balancing/scaling transformations. The first, symplectic scaling, is of the form -1 ( A' G' ) (**) H' = D H D = ( T ), ( Q' -A' ) -1 where D = diag(D0,D0 ) is diagonal (and hence symplectic). D0 is chosen to give the rows and columns of A' approximately equal 1-norms and so that G' and Q' have equal 1-norms. The latter condition implies that the bound on the right hand side of T || H' || <= max(||A'||,||A' ||) + max(||G'/r||,||r*Q'||) is minimized over all r > 0. Alternatively, DHABL may apply a norm-scaling transformation 1 ( A G/tau) (***) H' = --- ( T), tau (tau*Q -A ) where tau is a real scalar chosen as tau = max ( 1, ||A||, ||G||, ||Q|| ). Note that this is not a similarity transformation. The eigenvalues of H are tau times the eigenvalues of H'. DHASRD transforms a Hamiltonian matrix H' to square-reduced form by a symplectic orthogonal similarity transformation, T ( A'' G'' ) (#) H'' = U H' U = ( T ), ( Q'' -A'' ) where U is orthogonal and symplectic chosen so that 2 ( K1 K2 ) (##) H'' = K = ( T ), ( 0 K1 ) 2 where K1 = A'' + G'' Q'' is upper Hessenberg and K2 is skew-symmetric. The eigenvalues of H are the positive and negative square roots of the eigenvalues of K1. DHAEVS completes the eigenvalue calculation by accepting the output of DHASRD, forming the Hessenberg matrix K1 and calculating this eigenvalues by the Hessenberg QR algorithm. DHAEVD is an easy-to-use driver which calls DHASRD and DHAEVS. In order to make it as simple to use as possible, no scaling options are provided. --------------------------------------------------------------------------- IMPLEMENTATION: =============== DHABL, DHASRD, DHAEVS, and DHAEVD are implemented in the style of LAPACK [3], following the SLICOT (Subroutine Library in Control Theory) Implementation and Documentation Standard [4]. To date, only DOUBLE PRECISION versions of these subroutines are available. The subroutines require the DOUBLE PRECISION versions of the BLAS and LAPACK [3]. Thus, to run the codes, you have to link the BLAS and LAPACK libraries. If those are not available on your computer, they can be retrieved through netlib, http://www.netlib.org. These are the calling sequences for DHABL, DHASRD, DHAEVS, and DHAEVD. SUBROUTINE DHABL(JOBSCL, N, A, LDA, QG, LDQG, D, RWORK, IERR) INTEGER IERR, LDA, LDQG, N CHARACTER*1 JOBSCL DOUBLE PRECISION A(LDA,N), D(*), QG(LDQG,N+1), RWORK(N) SUBROUTINE DHASRD(COMPU, N, A, LDA, QG, LDQG, U, LDU, RWORK, IERR) INTEGER IERR, LDA, LDQG, LDU, N CHARACTER*1 COMPU DOUBLE PRECISION A(LDA,N), QG(LDQG,N+1), RWORK(2*N), U(LDU,*) SUBROUTINE DHAEVS(JOBSCL, N, A, LDA, QG, LDQG, WR, WI, RWORK, $ IERR) INTEGER IERR, LDA, LDQG, N CHARACTER*1 JOBSCL DOUBLE PRECISION A(LDA,N), QG(LDQG,N+1), RWORK(N*(N+1)), WI(N), $ WR(N) SUBROUTINE DHAEVD(N, A, LDA, QG, LDQG, WR, WI, RWORK, IERR) INTEGER IERR, LDA, LDQG, N DOUBLE PRECISION A(LDA,N), QG(LDQG,N+1), RWORK(N*(N+1)), WI(N), $ WR(N) Detailed descriptions of argument lists appear in prolog of the subroutines source and in [1]. The Hamiltonian matrix H is presented by its blocks A, G, and Q of (*). The symmetric matrices G and Q are packed together in the n-by-(n+1) array QG. The lower triangle of Q occupies the lower triangle of QG, and the upper triangle of G occupies the upper right triangle of QG. Thus, if i >= j, then Q(i,j) is stored in QG(i,j) and G(j,i) is stored in QG(j,i+1). DHABL returns the scaled Hamiltonian matrix (**) or (***) as directed by JOBSCL. DHASRD returns the transformed Hamiltonian matrix H' of (#) represented by its blocks A', G', and Q' which are returned in arrays A and QG. DHAEVS accepts the output from DHASRD and calculates the eigenvalues of (*) as the square roots of the eigenvalues of K1 in (##). The eigenvalues are returned in the arrays WR (real parts) and WI (imaginary parts) making use of the symmetry of the spectrum of Hamiltonian matrices, i.e., only N eigenvalues with nonnegative parts are returned. The other N eigenvalues are the negatives of the returned eigenvalues. DHAEVD is an easy-to-use driver. Sample and validation programs appear along with DHABL, DHASRD, DHAEVS, and DHAEVD in directory "test". --------------------------------------------------------------------------- CONTENTS: ========= README - This file. source - Directory of fortran source codes. dcroot.f - The subroutine DCROOT computes the square root of a complex number using real arithmetic. This was adapted from the EISPACK subroutine CSROOT. dhaevs.f - Fortran source of DHAEVS. dhaevd.f - Fortran source of DHAEVD. dhasrd.f - Fortran source of DHASRD. dhabl.f - Fortran source of DHABL. test - Directory of fortran source codes, sample programs, and validation programs. The source codes from directory "source" are duplicated. Makefile - A sample makefile for 1) compiling DHABL, DHASRD, DHAEVS, DHAEVD; 2) compiling the example programs and validation programs; 3) Running the example and validation programs. The output overwrites the ...output files described below. example.data - Data read by example.f. example.f - Fortran source of a simple example demonstrating the use of DHABL, DHASRD, and DHAEVS to calculate the eigenvalues of a Hamiltonian matrix. example.output - Sample output from example produced by example < example.data. dhaevs_example.data - Data read by dhaevs_example. dhaevs_example.f - Fortran source of a simple example demonstrating the use of DHAEVS. dhaevs_example.output - Sample output from dhaevs_example produced by dhaevs_example < dhaevs_example.data. dhaevs_test.f - Fortran source of a program to validate DHAEVS. The output is self explanatory. dhaevs_test.output - Sample output from a run of dhaevs_test. dhaevd_example.f - Fortran source of a simple example demonstrating the use of DHAEVD. dhaevd_example.output - Sample output from dhaevd_example produced by dhaevd_example < dhaevd_example.data dhaevd_test.f - Fortran source of a program to validate DHAEVD. The output is self explanatory. dhaevd_test.output - Sample output from a run of dhaevd_test. dhasrd_example.data - Data read by dhasrd_example. dhasrd_example.f - Fortran source of a simple example demonstrating the use of DHASRD. dhasrd_example.output - Sample output from a run of dhasrd_example produced by dhasrd_example < dhasrd_example.data. dhasrd_test.f - Fortran source of a program to validate DHASRD. The output is self explanatory. dhasrd_test.output - Sample output from a run of dhasrd_test. dhabl_example.data - Data read by dhabl_example. dhabl_example.f - Fortran source of a simple example demonstrating the use of DHABL. dhabl_example.output - Sample output from dhabl_example produced by dhabl_example < dhabl_example.data. dhabl_test.f - Fortran source of a program to validate DHABL. The output is self explanatory. dhabl_test.output - Sample output from dhabl_test. --------------------------------------------------------------------------- HELP AND BUGS: ============== If you have trouble either compiling or running the codes or if you find any bugs please send an e-mail message reporting the problem or bug to benner@math.uni-bremen.de We will get in touch with you as soon as possible. --------------------------------------------------------------------------- REFERENCES: =========== [1] P. Benner, R. Byers, and E. Barth: 'Fortran 77 Subroutines for Computing the Eigenvalues of Hamiltonian Matrices {I}: The Square Reduced Method', submitted to ACM Trans. Math. Software, 1998. [2] C. Van Loan: 'A Symplectic Method for Approximating all the Eigenvalues of a Hamiltonian Matrix', Linear Algebra and its Applications, Vol. 16, pp. 233-251 (1984). [3] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov and D. Sorensen: 'LAPACK Users' Guide', 2nd edition, SIAM, Philadelphia, PA, 1994. [4] The Working Group on Software: WGS 'SLICOT Implementation and Documentation Standards 2.1', WGS Report 96-1, available from http://www.win.tue.nl/wgs/reports.html. --------------------------------------------------------------------------- CONTRIBUTORS: ============= Peter Benner Zentrum f"ur Technomathematik Fachbereich 3 - Mathematik und Informatik Universit"at Bremen 28334 Bremen (FRG) e-mail: benner@math.uni-bremen.de Ralph Byers Department of Mathematics University of Kansas 405 Snow Hall Lawrence, KS 66045 (USA) e-mail: byers@math.ukans.edu Eric Barth Department of Mathematics Kalamazoo College 1200 Academy Street Kalamazoo, MI 49006 (USA) e-mail: barth@kzoo.edu --------------------------------------------------------------------------- SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << SHAR_EOF > 'data1' DHABL/DHASRD/DHAEVS EXAMPLE PROGRAM DATA 3 'A' 2.0 0.0 0.0 0.0 1.0 2.0 0.0 -1.0 3.0 1.0 0.0 0.0 2.0 3.0 4.0 -2.0 10.0 10.0 10.0 10.0 10.0 SHAR_EOF fi # end of overwriting check if test -f 'data2' then echo shar: will not over-write existing file "'data2'" else cat << SHAR_EOF > 'data2' DHABL EXAMPLE PROGRAM DATA 3 'A' -0.4 0.05 0.0007 -4.7 0.8 0.025 81.0 29.0 -0.9 0.0034 0.0014 0.00077 -0.005 0.0004 0.003 -18.0 -12.0 43.0 99.0 420.0 -200.0 SHAR_EOF fi # end of overwriting check if test -f 'data4' then echo shar: will not over-write existing file "'data4'" else cat << SHAR_EOF > 'data4' DHAEV EXAMPLE PROGRAM DATA 3 2.0 0.0 0.0 0.0 1.0 2.0 0.0 -1.0 3.0 1.0 0.0 0.0 2.0 3.0 4.0 -2.0 10.0 10.0 10.0 10.0 10.0 SHAR_EOF fi # end of overwriting check if test -f 'data6' then echo shar: will not over-write existing file "'data6'" else cat << SHAR_EOF > 'data6' DHAEVS EXAMPLE PROGRAM DATA 3 'S' 2.0 0.0 0.0 0.0 1.0 2.0 0.0 -1.0 3.0 1.0 0.0 0.0 2.0 3.0 4.0 -2.0 0.0 0.0 0.0 0.0 0.0 SHAR_EOF fi # end of overwriting check if test -f 'data8' then echo shar: will not over-write existing file "'data8'" else cat << SHAR_EOF > 'data8' DHASRD EXAMPLE PROGRAM DATA 3 'N' 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 1.0 1.0 1.0 2.0 2.0 3.0 7.0 6.0 5.0 8.0 4.0 9.0 SHAR_EOF fi # end of overwriting check if test -f 'driver1.f' then echo shar: will not over-write existing file "'driver1.f'" else cat << SHAR_EOF > 'driver1.f' * DHABL/DHASRD/DHAEVS EXAMPLE PROGRAM TEXT. * * Find the eigenvalues of a Hamiltonian matrix by the square reduced * method using symplectic balancing, and norm balancing. * C .. Parameters .. * NIN is the standard input. * NOUT is the standard output, e.g., a terminal INTEGER NIN,NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDA,LDQG PARAMETER (LDA=NMAX,LDQG=NMAX) C .. C .. Local Scalars .. INTEGER I,IERR,IGNORE,J,N CHARACTER JOBSCL C .. C .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX),D(NMAX),DUMMY(1),QG(LDQG,NMAX+1), + RWORK(NMAX* (NMAX+1)),WI(NMAX),WR(NMAX) C .. C .. External Subroutines .. * ************************************************************************ EXTERNAL DHABL,DHAEVS,DHASRD C .. * WRITE (NOUT,FMT=9000) * Skip the heading in the data file and read the data. READ (NIN,FMT='()') READ (NIN,FMT=*) N,JOBSCL IF (N.LE.0 .OR. N.GT.NMAX) THEN WRITE (NOUT,FMT=9010) N ELSE READ (NIN,FMT=*) ((A(I,J),J=1,N),I=1,N) READ (NIN,FMT=*) ((QG(J,I+1),I=J,N),J=1,N) READ (NIN,FMT=*) ((QG(I,J),I=J,N),J=1,N) * .. CALL DHABL(JOBSCL,N,A,LDA,QG,LDQG,D,RWORK,IGNORE) CALL DHASRD('N',N,A,LDA,QG,LDQG,DUMMY,1,RWORK,IGNORE) CALL DHAEVS('S',N,A,LDA,QG,LDQG,WR,WI,RWORK,IERR) IF (IERR.NE.0) THEN WRITE (NOUT,FMT=9020) IERR ELSE * .. * .. Show the computed eigenvalues. .. * Show the computed eigenvalues. WRITE (NOUT,FMT=9030) DO 10 I = 1,N WRITE (NOUT,FMT=9040) WR(I),' + (',WI(I),')i' 10 CONTINUE DO 20 I = N,1,-1 WRITE (NOUT,FMT=9040) - WR(I),' + (',-WI(I),')i' 20 CONTINUE * END IF END IF STOP * 9000 FORMAT (' DHABL/DHASRD/DHAEVS EXAMPLE PROGRAM RESULTS',/,1X) 9010 FORMAT (/,' N is out of range.',/,' N = ',I5) 9020 FORMAT (' IERR on exit from DHAEVS = ',I2) 9030 FORMAT (/,' The eigenvalues are ') 9040 FORMAT (1X,F8.4,A,F8.4,A) END SHAR_EOF fi # end of overwriting check if test -f 'driver2.f' then echo shar: will not over-write existing file "'driver2.f'" else cat << SHAR_EOF > 'driver2.f' * DHABL EXAMPLE PROGRAM TEXT. * C .. Parameters .. * NIN is the standard input. * NOUT is the standard output, e.g., a terminal INTEGER NIN,NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDA,LDQG PARAMETER (LDA=NMAX,LDQG=NMAX) INTEGER LRWORK PARAMETER (LRWORK=NMAX) C .. C .. Local Scalars .. INTEGER I,IERR,J,N CHARACTER JOBSCL C .. C .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX),D(NMAX),QG(LDQG,NMAX+1),RWORK(LRWORK) C .. C .. External Subroutines .. EXTERNAL DHABL C .. C .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * WRITE (NOUT,FMT=9000) * Skip the heading in the data file and read the data. READ (NIN,FMT='()') READ (NIN,FMT=*) N,JOBSCL IF (N.LE.0 .OR. N.GT.NMAX) THEN WRITE (NOUT,FMT=9010) N ELSE READ (NIN,FMT=*) ((A(I,J),J=1,N),I=1,N) READ (NIN,FMT=*) ((QG(J,I+1),I=J,N),J=1,N) READ (NIN,FMT=*) ((QG(I,J),I=J,N),J=1,N) END IF * .. scale the Hamiltonian matrix .. CALL DHABL(JOBSCL,N,A,LDA,QG,LDQG,D,RWORK,IERR) IF (IERR.NE.0) THEN WRITE (NOUT,FMT=9020) IERR ELSE * .. show the scaled Hamiltonian matrix.. WRITE (NOUT,FMT=9030) DO 10 I = 1,N WRITE (NOUT,FMT=9060) (A(I,J),J=1,N), + (QG(J,I+1),J=1,I-1), (QG(I,J+1),J=I,N) 10 CONTINUE DO 20 I = 1,N WRITE (NOUT,FMT=9060) (QG(I,J),J=1,I-1), (QG(J,I),J=I,N), + (-A(J,I),J=1,N) 20 CONTINUE * .. show the scaling factors.. IF (LSAME(JOBSCL,'A')) THEN WRITE (NOUT,FMT=9040) WRITE (NOUT,FMT=9060) (D(I),I=1,N) ELSE IF (LSAME(JOBSCL,'B')) THEN WRITE (NOUT,FMT=9050) WRITE (NOUT,FMT=9060) D(1) END IF END IF STOP * 9000 FORMAT (' DHABL EXAMPLE PROGRAM RESULTS',/,1X) 9010 FORMAT (/,' N is out of range.',/,' N = ',I5) 9020 FORMAT (' IERR on exit from DHABL = ',I2) 9030 FORMAT (/,' The scaled Hamiltonian is ') 9040 FORMAT (/,' The scaling factors are ') 9050 FORMAT (/,' The scaling factor tau is ') 9060 FORMAT (1X,8 (F10.4)) END SHAR_EOF fi # end of overwriting check if test -f 'driver3.f' then echo shar: will not over-write existing file "'driver3.f'" else cat << SHAR_EOF > 'driver3.f' * test driver for DHABL.f * ***************************************************************** * This program scales several ad-hoc Hamiltonian matrices, then * verifies that the scaled matrices are similar to the originals * as they should be. * ***************************************************************** C .. Parameters .. * NOUT is the standard output, e.g., a terminal INTEGER NOUT PARAMETER (NOUT=6) INTEGER NMAX PARAMETER (NMAX=15) INTEGER LDA,LDQG PARAMETER (LDA=NMAX,LDQG=NMAX) DOUBLE PRECISION ZERO,ONE,TWO PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0) CHARACTER*(*) EQULS PARAMETER (EQULS='=============================================') C .. C .. Local Scalars .. * DOUBLE PRECISION ANRM,BASE,EPS,ERRA,ERRG,ERRQ,GNRM,GPNRM,QNRM, + QPNRM,R,SIMERR,TAU,Y INTEGER I,IERR,IFUDGE,INFO,IREP,J,N,NFAILS LOGICAL BBAL,FAILED,NORM,SYMP CHARACTER JOBSCL C .. C .. Local Arrays .. * DOUBLE PRECISION A(LDA,LDA),AP(LDA,LDA),D(NMAX),QG(LDQG,LDQG+1), + QPGP(LDQG,LDQG+1),RWORK(NMAX) INTEGER ISEED(4) C .. C .. External Functions .. * DOUBLE PRECISION DLAMCH,DLANGE,DLANSY,DLAPY2 EXTERNAL DLAMCH,DLANGE,DLANSY,DLAPY2 C .. C .. External Subroutines .. * EXTERNAL DHABL,DLACPY,DLARNV,DLASCL,DRSCL,DSCAL C .. C .. Intrinsic Functions .. * INTRINSIC DBLE,MAX C .. C .. Data statements .. * . (seed for random number generation) . * * ****************************************************************** DATA ISEED/86,1967,2001,1995/ C .. EPS = DLAMCH('P') BASE = DLAMCH('B') * WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) 'TESTING SYMPLECTIC SCALING' NFAILS = 0 DO 100 IREP = 1,2 DO 90 N = -1,NMAX WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) 'DHAEV TESTS: ORDER N = ',N IF (N.LT.0) THEN WRITE (NOUT,FMT=*) + 'OF COURSE, THIS IS NOT A VALID TEST CASE.' WRITE (NOUT,FMT=*) + 'IT CHECKS THAT DHAEV RETURNS WITH AN ERROR.' END IF * * .. set up Hamiltonian .. DO 10 I = 1,N CALL DLARNV(3,ISEED,N,A(1,I)) CALL DLARNV(3,ISEED,N,QG(1,I)) 10 CONTINUE DO 20 I = 1,N CALL DRSCL(N,TWO** (-I),A(I,1),LDA) CALL DSCAL(N,TWO** (-I),A(1,I),1) 20 CONTINUE IF (N.GT.0) THEN CALL DLARNV(3,ISEED,N,QG(1,N+1)) CALL DLASCL('U',0,0,DBLE(10*N),ONE,N,N,QG(1,2),LDQG, + INFO) CALL DLASCL('L',0,0,ONE,DBLE(10*N),N,N,QG,LDQG,INFO) END IF CALL DLACPY('A',N,N,A,LDA,AP,LDA) CALL DLACPY('A',N,N+1,QG,LDQG,QPGP,LDQG) * * .. symplectic scaling .. IF (IREP.EQ.1) THEN WRITE (NOUT,FMT=*) 'SYMPLECTIC SCALING' JOBSCL = 'A' SYMP = .TRUE. NORM = .FALSE. ELSE WRITE (NOUT,FMT=*) 'NORM SCALING' JOBSCL = 'B' SYMP = .FALSE. NORM = .TRUE. END IF CALL DHABL(JOBSCL,N,AP,LDA,QPGP,LDQG,D,RWORK,IERR) * * ********************************************************** * .. Has DHABL something to tell? ===== IF (IERR.EQ.0) THEN FAILED = .FALSE. WRITE (NOUT,FMT=*) 'NO ERRORS DETECTED.' IF (N.GT.0) TAU = D(1) ELSE WRITE (NOUT,FMT=*) 'ERROR: DHABL RETURNS IERR = ',IERR FAILED = .TRUE. END IF IF (.NOT.FAILED) THEN * ***************************************************** * .. Nonsingularity check .. DO 30 I = 1,N IF (D(I).EQ.ZERO) THEN FAILED = .TRUE. WRITE (NOUT,FMT=*) 'ERROR IN DHABL: D(',I, + ') = 0' END IF IF (NORM) GO TO 40 30 CONTINUE 40 CONTINUE END IF IF (.NOT.FAILED .AND. N.GT.0) THEN * ***************************************************** * .. block scaling check .. * IF (NORM) THEN ANRM = DLANGE('1',N,N,A,LDA,RWORK) GNRM = DLANSY('1','U',N,QG(1,2),LDQG,RWORK) QNRM = DLANSY('1','L',N,QG,LDQG,RWORK) Y = MAX(ONE,ANRM,GNRM,QNRM) IF (Y.LE.TAU*BASE .AND. Y.GE.TAU/BASE) THEN WRITE (NOUT,FMT=*) 'SCALE OK:' ELSE WRITE (NOUT,FMT=*) + 'ERROR IN DHABL: OUT OF SCALE' FAILED = .TRUE. END IF WRITE (NOUT,FMT=*) ' TAU = ',TAU WRITE (NOUT,FMT=*) ' ANRM = ',ANRM WRITE (NOUT,FMT=*) ' GNRM = ',GNRM WRITE (NOUT,FMT=*) ' QNRM = ',QNRM ELSE GPNRM = DLANSY('1','U',N,QPGP(1,2),LDQG,RWORK) QPNRM = DLANSY('1','L',N,QPGP,LDQG,RWORK) BBAL = (GPNRM.LT.QPNRM*BASE) .AND. + (QPNRM/BASE.LT.GPNRM) BBAL = BBAL .OR. (QPNRM.EQ.ZERO) .OR. + (GPNRM.EQ.ZERO) IF (BBAL) THEN WRITE (NOUT,FMT=*) 'SCALE OK: GPNRM = ', + GPNRM WRITE (NOUT,FMT=*) ' QPNRM = ', + QPNRM ELSE WRITE (NOUT,FMT=*) + 'ERROR IN DHABL: OUT OF SCALE' WRITE (NOUT,FMT=*) 'GP NORM = ',GPNRM WRITE (NOUT,FMT=*) 'QP NORM = ',QPNRM FAILED = .TRUE. END IF END IF END IF IF (.NOT.FAILED) THEN * ***************************************************** * .. Similarity check (modulo rounding) .. * ERRA = ZERO DO 60 J = 1,N DO 50 I = 1,N IF (NORM) R = A(I,J) - AP(I,J)*TAU IF (SYMP) R = A(I,J) - AP(I,J)*D(I)/D(J) IF (A(I,J).NE.ZERO) THEN R = R/A(I,J) ELSE R = R/ANRM END IF ERRA = DLAPY2(ERRA,R) 50 CONTINUE 60 CONTINUE ERRQ = ZERO ERRG = ZERO DO 80 J = 1,N DO 70 I = J,N IF (NORM) R = QG(I,J) - QPGP(I,J) IF (SYMP) R = QG(I,J) - QPGP(I,J)/D(I)/D(J) IF (QG(I,J).NE.ZERO) THEN R = R/QG(I,J) ELSE R = R/QNRM END IF ERRQ = DLAPY2(ERRQ,R) IF (NORM) R = QG(J,I+1) - QPGP(J,I+1)*TAU*TAU IF (SYMP) R = QG(J,I+1) - + QPGP(J,I+1)*D(I)*D(J) IF (QG(J,I+1).NE.ZERO) THEN R = R/QG(J,I+1) ELSE R = R/GNRM END IF ERRG = DLAPY2(ERRG,R) 70 CONTINUE 80 CONTINUE SIMERR = DLAPY2(DLAPY2(ERRA,ERRG),DLAPY2(ERRQ,ERRA)) IFUDGE = 2*N * WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) 'SIMILARITY CHECK:' WRITE (NOUT,FMT=*) 'RELATIVE ERROR = ',SIMERR * IF (SIMERR.LE.DBLE(IFUDGE)*EPS) THEN WRITE (NOUT,FMT=*) + 'OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE' WRITE (NOUT,FMT=*) + ' INPUT ''MODULO ROUNDING ERRORS'' AS IT ' WRITE (NOUT,FMT=*) ' SHOULD BE.' ELSE FAILED = .TRUE. WRITE (NOUT,FMT=*) + 'SUSPICIOUS: THE RELATIVE ERROR OF SIMILARITY' WRITE (NOUT,FMT=*) + ' IS LARGER THAN EXPECTED. ' END IF END IF IF (FAILED .AND. N.GE.0) NFAILS = NFAILS + 1 WRITE (NOUT,FMT=*) 90 CONTINUE WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) 100 CONTINUE * ****************************************************************** WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS IF (NFAILS.EQ.0) THEN WRITE (NOUT,FMT=*) 'NO SUSPICIOUS TEST RESULTS.' ELSE WRITE (NOUT,FMT=*) 'SIGH..., ',NFAILS, + ' SUSPICIOUS TEST RESULTS.' END IF STOP END SHAR_EOF fi # end of overwriting check if test -f 'driver4.f' then echo shar: will not over-write existing file "'driver4.f'" else cat << SHAR_EOF > 'driver4.f' * DHAEVD EXAMPLE PROGRAM TEXT. * C .. Parameters .. * NIN is the standard input. * NOUT is the standard output, e.g., a terminal INTEGER NIN,NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDA,LDQG PARAMETER (LDA=NMAX,LDQG=NMAX) INTEGER LRWORK PARAMETER (LRWORK=NMAX* (NMAX+1)) C .. C .. Local Scalars .. INTEGER I,IERR,J,N C .. C .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX),QG(LDQG,NMAX+1),RWORK(LRWORK), + WI(NMAX),WR(NMAX) C .. C .. External Subroutines .. EXTERNAL DHAEVD C .. * WRITE (NOUT,FMT=9000) * Skip the heading in the data file and read the data. * NOTE: input must define a square-reduced Hamiltonian matrix. READ (NIN,FMT='()') READ (NIN,FMT=*) N IF (N.LE.0 .OR. N.GT.NMAX) THEN WRITE (NOUT,FMT=9010) N ELSE READ (NIN,FMT=*) ((A(I,J),J=1,N),I=1,N) READ (NIN,FMT=*) ((QG(J,I+1),I=J,N),J=1,N) READ (NIN,FMT=*) ((QG(I,J),I=J,N),J=1,N) * Compute the eigenvalues. CALL DHAEVD(N,A,LDA,QG,LDQG,WR,WI,RWORK,IERR) * IF (IERR.NE.0) THEN WRITE (NOUT,FMT=9020) IERR ELSE * Show the computed eigenvalues. WRITE (NOUT,FMT=9030) DO 10 I = 1,N WRITE (NOUT,FMT=9040) WR(I),' + (',WI(I),')i' 10 CONTINUE DO 20 I = N,1,-1 WRITE (NOUT,FMT=9040) - WR(I),' + (',-WI(I),')i' 20 CONTINUE END IF END IF STOP * 9000 FORMAT (' DHAEVD EXAMPLE PROGRAM RESULTS',/,1X) 9010 FORMAT (/,' N is out of range.',/,' N = ',I5) 9020 FORMAT (' IERR on exit from DHAEVD = ',I2) 9030 FORMAT (/,' The eigenvalues are ') 9040 FORMAT (1X,F8.4,A,F8.4,A) END SHAR_EOF fi # end of overwriting check if test -f 'driver5.f' then echo shar: will not over-write existing file "'driver5.f'" else cat << SHAR_EOF > 'driver5.f' * test driver for DHAEVD.f * ***************************************************************** * * This program sets up several ad-hoc Hamiltonian matrices and then * calls DHAEVD to find their eigenvalues. The program verifies * that the nonzero smallest singular values sigma_min(H - lambda I) * are can be accounted for by rounding errors. * * ****************************************************************** C .. Parameters .. * NOUT is the standard output, e.g., a terminal INTEGER NOUT PARAMETER (NOUT=6) INTEGER NMAX,NMIN PARAMETER (NMAX=15,NMIN=-1) INTEGER LD PARAMETER (LD=NMAX) INTEGER LDA,LDQG,LDH PARAMETER (LDA=LD,LDQG=LD,LDH=2*LD) INTEGER LWORK PARAMETER (LWORK=5*LDH) DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) CHARACTER*(*) EQULS PARAMETER (EQULS=' =============================================') C .. C .. Local Scalars .. DOUBLE PRECISION ANRM,GNRM,HNRM,QNRM,RATMAX,RTEPS INTEGER I,IERR,INFO,J,K,MULT,N,NFAILS,NUNSRT C .. C .. Local Arrays .. DOUBLE COMPLEX HP(LDH,LDH),U(1),VT(1),ZWORK(LWORK) DOUBLE PRECISION A(LDA,LDA),H(LDH,LDH),QG(LDQG,LDQG+1), + RWORK(10*LD),S(LDH),WI(2*LD),WORK(LD,LD),WR(2*LD) INTEGER ISEED(4) C .. C .. External Functions .. DOUBLE PRECISION DLAMCH,DLANGE,DLANSY,DLAPY2 EXTERNAL DLAMCH,DLANGE,DLANSY,DLAPY2 C .. C .. External Subroutines .. EXTERNAL DCOPY,DHAEVD,DLACPY,DLARNV,DLASCL,ZGESVD C .. C .. Intrinsic Functions .. INTRINSIC DCMPLX,MAX,SQRT C .. C .. Data statements .. * . (seed for random number generation) . * DATA ISEED/86,1967,2001,1995/ C .. * ****************************************************************** RTEPS = SQRT(DLAMCH('P')) * NFAILS = 0 NUNSRT = 0 MULT = 0 * ==== 2 test runs: with and without Hessenberg scaling ===== WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS DO 90 N = NMIN,NMAX WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) 'DHAEVD TEST: ORDER N = ',N IF (N.LT.0) THEN WRITE (NOUT,FMT=*) + 'OF COURSE, THIS IS NOT A VALID TEST CASE.' WRITE (NOUT,FMT=*) + 'IT CHECKS THAT DHAEVD RETURNS WITH ERROR IERR = -1.' END IF * * ==== set up Hamiltonian ==== DO 10 I = 1,N CALL DLARNV(3,ISEED,N,A(1,I)) CALL DLARNV(3,ISEED,N,QG(1,I)) 10 CONTINUE IF (N.GE.0) THEN CALL DLARNV(3,ISEED,N,QG(1,N+1)) CALL DLACPY('A',N,N,A,LDA,H,LDH) CALL DLACPY('L',N,N,QG,LDQG,H(N+1,1),LDH) CALL DLACPY('U',N,N,QG(1,2),LDQG,H(1,N+1),LDH) END IF DO 20 J = 1,N CALL DCOPY(N,A(1,J),1,H(N+J,N+1),LDH) CALL DCOPY(J,QG(J,1),LDQG,H(N+1,J),1) CALL DCOPY(J,QG(1,J+1),1,H(J,N+1),LDH) 20 CONTINUE IF (N.GT.0) CALL DLASCL('G',0,0,ONE,-ONE,N,N,H(N+1,N+1),LDH, + INFO) * * ==== Frobenius norm of Hamiltonian Matrix ==== ANRM = DLANGE('F',N,N,A,LD,RWORK) GNRM = DLANSY('F','L',N,QG,LDQG,RWORK) QNRM = DLANSY('F','U',N,QG(1,2),LDQG,RWORK) HNRM = DLAPY2(DLAPY2(ANRM,GNRM),DLAPY2(ANRM,QNRM)) * * ==== find eigenvalues ==== CALL DHAEVD(N,A,LDA,QG,LDQG,WR,WI,WORK,IERR) * * ********************************************************** * ==== Test 0: Error exit from DHAEVD? ===== IF (IERR.NE.0) WRITE (NOUT,FMT=*) 'DHAEVD RETURNS IERR = ', + IERR IF (IERR.GT.0) WRITE (NOUT,FMT=* + ) 'DHSEQR FAILED TO CONVERGE WHILE COMPUTING THE',J, + 'TH EIGENVALUE. (THIS IS VERY RARE.)' IF (IERR.EQ.-1) WRITE (NOUT,FMT=*) 'N = ',N,' IS NEGATIVE.' IF (IERR.EQ.-3) WRITE (NOUT,FMT=*) 'LDA = ',LDA, + ' IS LESS THAN N = ',N,'.' IF (IERR.EQ.-5) WRITE (NOUT,FMT=*) 'LDQG = ',LDQG, + ' IS LESS THAN N = ',N,'.' IF (IERR.EQ.0) THEN WRITE (NOUT,FMT=*) 'NO ERRORS WERE DETECTED.' * ****************************************************** * ==== Test 1: Sort check ===== DO 30 I = 1,N - 1 IF ((WR(I).LT.WR(I+1)) .OR. + ((WR(I).EQ.ZERO).AND. (WR(I+1).EQ.ZERO).AND. + (WI(I).LT.WI(I+1)))) THEN WRITE (NOUT,FMT=*) 'EIGENVALUE OUT OF ORDER!' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ',WR(I+1) NUNSRT = NUNSRT + 1 GO TO 50 END IF 30 CONTINUE DO 40 I = 1,N - 1 IF (WR(I).EQ.WR(I+1) .AND. WR(I).NE.ZERO) THEN IF (WI(I).GT.ZERO .AND. WI(I).NE.-WI(I+1)) THEN WRITE (NOUT,FMT=*) + 'COMPLEX CONJUGATES NOT ADJACENT OR ', + 'OUT OF ORDER.' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I), + ' WI(',I,') = ',WI(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ',WR(I+1), + ' WI(',I + 1,') = ',WI(I+1) NUNSRT = NUNSRT + 1 GO TO 50 END IF IF (WR(I).EQ.WR(I+1) .AND. WI(I).EQ.WI(I+1)) THEN WRITE (NOUT,FMT=*) + 'PECULIAR! THERE ARE MULTIPLE EIGENVALUES.' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I), + ' WI(',I,') = ',WI(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ',WR(I+1), + ' WI(',I + 1,') = ',WI(I+1) MULT = MULT + 1 END IF END IF 40 CONTINUE 50 CONTINUE * ==== Test 2: eigenvalue check ==== RATMAX = ZERO DO 80 K = 1,N DO 70 J = 1,2*N DO 60 I = 1,2*N HP(I,J) = H(I,J) 60 CONTINUE HP(J,J) = HP(J,J) - DCMPLX(WR(K),WI(K)) 70 CONTINUE CALL ZGESVD('N','N',2*N,2*N,HP,LDH,S,U,1,VT,1,ZWORK, + LWORK,RWORK,INFO) RATMAX = MAX(RATMAX,S(2*N)/HNRM) IF (S(2*N).GT.RTEPS*HNRM) THEN WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) 'SUSPICIOUS RATIO' WRITE (NOUT,FMT=*) 'COMPUTED EIGENVALUE = ',WR(K), + ' + ',WI(K),'I' WRITE (NOUT,FMT=*) 'SIGMA-MIN(H-LAMBDA*I) = ', + S(2*N) WRITE (NOUT,FMT=*) END IF 80 CONTINUE WRITE (NOUT,FMT=*) 'MAXIMUM RATIO = ',RATMAX IF (RATMAX.GT.RTEPS) NFAILS = NFAILS + 1 END IF 90 CONTINUE * ****************************************************************** WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) IF (NFAILS.EQ.0 .AND. NUNSRT.EQ.0) WRITE (NOUT, + FMT=*) 'NO SUSPICIOUS TEST RESULTS.' IF (MULT.NE.0) WRITE (NOUT,FMT=*) 'PECULIAR! THERE WERE ',MULT, + ' MULTIPLE EIGENVALUES.' IF (NUNSRT.NE.0) WRITE (NOUT,FMT=*) 'SIGH..., THERE WERE ',NUNSRT, + ' EIGENVALUE SORTING FAILURES.' IF (NFAILS.NE.0) WRITE (NOUT,FMT=*) 'SIGH..., THERE WERE ',NFAILS, + ' SUSPICIOUS TEST RESULTS.' STOP END SHAR_EOF fi # end of overwriting check if test -f 'driver6.f' then echo shar: will not over-write existing file "'driver6.f'" else cat << SHAR_EOF > 'driver6.f' * DHAEVS EXAMPLE PROGRAM TEXT. * C .. Parameters .. * NIN is the standard input. * NOUT is the standard output, e.g., a terminal INTEGER NIN,NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDA,LDQG PARAMETER (LDA=NMAX,LDQG=NMAX) INTEGER LRWORK PARAMETER (LRWORK=NMAX* (NMAX+1)) C .. C .. Local Scalars .. INTEGER I,IERR,J,N CHARACTER JOBSCL C .. C .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX),QG(LDQG,NMAX+1),RWORK(LRWORK), + WI(NMAX),WR(NMAX) C .. C .. External Subroutines .. EXTERNAL DHAEVS C .. * WRITE (NOUT,FMT=9000) * Skip the heading in the data file and read the data. * NOTE: input must define a square-reduced Hamiltonian matrix. READ (NIN,FMT='()') READ (NIN,FMT=*) N,JOBSCL IF (N.LE.0 .OR. N.GT.NMAX) THEN WRITE (NOUT,FMT=9010) N ELSE READ (NIN,FMT=*) ((A(I,J),J=1,N),I=1,N) READ (NIN,FMT=*) ((QG(J,I+1),I=J,N),J=1,N) READ (NIN,FMT=*) ((QG(I,J),I=J,N),J=1,N) * Compute the eigenvalues. CALL DHAEVS(JOBSCL,N,A,LDA,QG,LDQG,WR,WI,RWORK,IERR) * IF (IERR.NE.0) THEN WRITE (NOUT,FMT=9020) IERR ELSE * Show the computed eigenvalues. WRITE (NOUT,FMT=9030) DO 10 I = 1,N WRITE (NOUT,FMT=9040) WR(I),' + (',WI(I),')i' 10 CONTINUE DO 20 I = N,1,-1 WRITE (NOUT,FMT=9040) - WR(I),' + (',-WI(I),')i' 20 CONTINUE END IF END IF STOP * 9000 FORMAT (' DHAEVS EXAMPLE PROGRAM RESULTS',/,1X) 9010 FORMAT (/,' N is out of range.',/,' N = ',I5) 9020 FORMAT (' IERR on exit from DHAEVS = ',I2) 9030 FORMAT (/,' The eigenvalues are ') 9040 FORMAT (1X,F8.4,A,F8.4,A) END SHAR_EOF fi # end of overwriting check if test -f 'driver7.f' then echo shar: will not over-write existing file "'driver7.f'" else cat << SHAR_EOF > 'driver7.f' * test driver for DHAEVS.f * ***************************************************************** * * This program sets up several ad-hoc Hamiltonian matrices and * transforms them to square-reduced form by Van Loan's method, then * calls DHAEVS to find their eigenvalues. The program verifies * that the nonzero smallest singular values sigma_min(H - lambda I) * are can be accounted for by rounding errors. * * ****************************************************************** C .. Parameters .. * NOUT is the standard output, e.g., a terminal INTEGER NOUT PARAMETER (NOUT=6) INTEGER NMAX,NMIN PARAMETER (NMAX=15,NMIN=-1) INTEGER LD PARAMETER (LD=NMAX) INTEGER LDA,LDQG,LDH PARAMETER (LDA=LD,LDQG=LD,LDH=2*LD) INTEGER LWORK PARAMETER (LWORK=5*LDH) DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) CHARACTER*(*) EQULS PARAMETER (EQULS=' =============================================') C .. C .. Local Scalars .. DOUBLE PRECISION ANRM,GNRM,HNRM,QNRM,RATMAX,RTEPS INTEGER I,IERRH,IERRS,INFO,ITER,J,K,LDU,MULT,N,NFAILS,NUNSRT CHARACTER JOBSCL C .. C .. Local Arrays .. DOUBLE COMPLEX HP(LDH,LDH),UD(1),VT(1),ZWORK(LWORK) DOUBLE PRECISION A(LDA,LDA),H(LDH,LDH),QG(LDQG,LDQG+1), + RWORK(10*LD),S(LDH),U(1),WI(LD),WORK(LD* (LD+1)), + WR(LD) INTEGER ISEED(4) C .. C .. External Functions .. DOUBLE PRECISION DLAMCH,DLANGE,DLANSY,DLAPY2 LOGICAL LSAME EXTERNAL DLAMCH,DLANGE,DLANSY,DLAPY2,LSAME C .. C .. External Subroutines .. EXTERNAL DCOPY,DHAEVS,DHASRD,DLACPY,DLARNV,DLASCL,ZGESVD C .. C .. Intrinsic Functions .. INTRINSIC DCMPLX,MAX,SQRT C .. C .. Data statements .. * . (seed for random number generation) . * DATA ISEED/86,1967,2001,1995/ C .. * ****************************************************************** RTEPS = SQRT(DLAMCH('P')) * NFAILS = 0 NUNSRT = 0 MULT = 0 DO 100 ITER = 1,2 IF (ITER.EQ.1) JOBSCL = 'N' IF (ITER.EQ.2) JOBSCL = 'S' * ==== 2 test runs: with and without Hessenberg scaling ===== WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS IF (LSAME(JOBSCL,'N')) WRITE (NOUT, + FMT=*) 'TEST RUN WITHOUT SCALING' IF (LSAME(JOBSCL,'S')) WRITE (NOUT, + FMT=*) 'TEST RUN WITH SCALING' DO 90 N = NMIN,NMAX WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) 'DHAEVS TEST: ORDER N = ',N IF (N.LT.0) THEN WRITE (NOUT,FMT=*) + 'OF COURSE, THIS IS NOT A VALID TEST CASE.' WRITE (NOUT,FMT=*) + 'IT CHECKS THAT DHAEVS RETURNS WITH ERROR IERR = -2.' END IF * * ==== set up Hamiltonian ==== DO 10 I = 1,N CALL DLARNV(3,ISEED,N,A(1,I)) CALL DLARNV(3,ISEED,N,QG(1,I)) 10 CONTINUE IF (N.GE.0) THEN CALL DLARNV(3,ISEED,N,QG(1,N+1)) CALL DLACPY('A',N,N,A,LDA,H,LDH) CALL DLACPY('L',N,N,QG,LDQG,H(N+1,1),LDH) CALL DLACPY('U',N,N,QG(1,2),LDQG,H(1,N+1),LDH) END IF DO 20 J = 1,N CALL DCOPY(N,A(1,J),1,H(N+J,N+1),LDH) CALL DCOPY(J,QG(J,1),LDQG,H(N+1,J),1) CALL DCOPY(J,QG(1,J+1),1,H(J,N+1),LDH) 20 CONTINUE IF (N.GT.0) CALL DLASCL('G',0,0,ONE,-ONE,N,N,H(N+1,N+1), + LDH,INFO) * * ==== Frobenius norm of Hamiltonian Matrix ==== ANRM = DLANGE('F',N,N,A,LD,RWORK) GNRM = DLANSY('F','L',N,QG,LDQG,RWORK) QNRM = DLANSY('F','U',N,QG(1,2),LDQG,RWORK) HNRM = DLAPY2(DLAPY2(ANRM,GNRM),DLAPY2(ANRM,QNRM)) * * ==== square reduce ==== CALL DHASRD('N',N,A,LDA,QG,LDQG,U,1,RWORK,IERRS) * ==== find eigenvalues ==== CALL DHAEVS(JOBSCL,N,A,LDA,QG,LDQG,WR,WI,WORK,IERRH) * * ********************************************************** * ==== Test 0: Error exit from DHASRD or DHAEVS? ===== IF (IERRS.NE.0) WRITE (NOUT, + FMT=*) 'DHASRD RETURNS IERR = ',IERRS IF (IERRS.EQ.-1) WRITE (NOUT, + FMT=*) 'COMPU IS NOT ONE OF ''A'', ''F'', OR ''N''.' IF (IERRS.EQ.-2) WRITE (NOUT,FMT=*) 'N = ',N, + ' IS NEGATIVE.' IF (IERRS.EQ.-4) WRITE (NOUT,FMT=*) 'LDA = ',LDA, + ' IS LESS THAN N = ',N,'.' IF (IERRS.EQ.-6) WRITE (NOUT,FMT=*) 'LDQG = ',LDQG, + ' IS LESS THAN N = ',N,'.' IF (IERRS.EQ.-8) WRITE (NOUT,FMT=*) 'LDU = ',LDU, + ' IS LESS THAN N = ',N,'.' IF (IERRH.NE.0) WRITE (NOUT, + FMT=*) 'DHAEVS RETURNS IERR = ',IERRH IF (IERRH.GT.0) WRITE (NOUT,FMT=*) + 'DHSEQR FAILED TO CONVERGE WHILE COMPUTING THE',J, + 'TH EIGENVALUE. (THIS IS VERY RARE.)' IF (IERRH.EQ.-1) WRITE (NOUT,FMT=*) 'JOBSCL = '//JOBSCL// + ' IS NOT ONE OF ''S'' OR ''N''.' IF (IERRH.EQ.-2) WRITE (NOUT,FMT=*) 'N = ',N, + ' IS NEGATIVE.' IF (IERRH.EQ.-4) WRITE (NOUT,FMT=*) 'LDA = ',LDA, + ' IS LESS THAN N = ',N,'.' IF (IERRH.EQ.-6) WRITE (NOUT,FMT=*) 'LDQG = ',LDQG, + ' IS LESS THAN N = ',N,'.' IF (IERRS.EQ.0 .AND. IERRH.EQ.0) THEN WRITE (NOUT,FMT=*) 'NO ERRORS WERE DETECTED.' * ****************************************************** * ==== Test 1: Sort check ===== DO 30 I = 1,N - 1 IF ((WR(I).LT.WR(I+1)) .OR. + ((WR(I).EQ.ZERO).AND. (WR(I+1).EQ.ZERO).AND. + (WI(I).LT.WI(I+1)))) THEN WRITE (NOUT,FMT=*) 'EIGENVALUE OUT OF ORDER!' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ',WR(I+1) NUNSRT = NUNSRT + 1 GO TO 50 END IF 30 CONTINUE DO 40 I = 1,N - 1 IF (WR(I).EQ.WR(I+1) .AND. WR(I).NE.ZERO) THEN IF (WI(I).GT.ZERO .AND. + WI(I).NE.-WI(I+1)) THEN WRITE (NOUT,FMT=*) + 'COMPLEX CONJUGATES NOT ADJACENT OR ', + 'OUT OF ORDER.' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I), + ' WI(',I,') = ',WI(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ', + WR(I+1),' WI(',I + 1,') = ',WI(I+1) NUNSRT = NUNSRT + 1 GO TO 50 END IF IF (WR(I).EQ.WR(I+1) .AND. + WI(I).EQ.WI(I+1)) THEN WRITE (NOUT,FMT=*) + 'PECULIAR! THERE ARE MULTIPLE EIGENVALUES.' WRITE (NOUT,FMT=*) 'WR(',I,') = ',WR(I), + ' WI(',I,') = ',WI(I) WRITE (NOUT,FMT=*) 'WR(',I + 1,') = ', + WR(I+1),' WI(',I + 1,') = ',WI(I+1) MULT = MULT + 1 END IF END IF 40 CONTINUE 50 CONTINUE * ==== Test 2: eigenvalue check ==== RATMAX = ZERO DO 80 K = 1,N DO 70 J = 1,2*N DO 60 I = 1,2*N HP(I,J) = H(I,J) 60 CONTINUE HP(J,J) = HP(J,J) - DCMPLX(WR(K),WI(K)) 70 CONTINUE CALL ZGESVD('N','N',2*N,2*N,HP,LDH,S,UD,1,VT,1, + ZWORK,LWORK,RWORK,INFO) RATMAX = MAX(RATMAX,S(2*N)/HNRM) IF (S(2*N).GT.RTEPS*HNRM) THEN WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) 'SUSPICIOUS RATIO' WRITE (NOUT,FMT=*) 'COMPUTED EIGENVALUE = ', + WR(K),' + ',WI(K),'I' WRITE (NOUT,FMT=*) 'SIGMA-MIN(H-LAMBDA*I) = ', + S(2*N) WRITE (NOUT,FMT=*) END IF 80 CONTINUE WRITE (NOUT,FMT=*) 'MAXIMUM RATIO = ',RATMAX IF (RATMAX.GT.RTEPS) NFAILS = NFAILS + 1 END IF 90 CONTINUE 100 CONTINUE * ****************************************************************** WRITE (NOUT,FMT=*) WRITE (NOUT,FMT=*) EQULS WRITE (NOUT,FMT=*) IF (NFAILS.EQ.0 .AND. NUNSRT.EQ.0) WRITE (NOUT, + FMT=*) 'NO SUSPICIOUS TEST RESULTS.' IF (MULT.NE.0) WRITE (NOUT,FMT=*) 'PECULIAR! THERE WERE ',MULT, + ' MULTIPLE EIGENVALUES.' IF (NUNSRT.NE.0) WRITE (NOUT,FMT=*) 'SIGH..., THERE WERE ',NUNSRT, + ' EIGENVALUE SORTING FAILURES.' IF (NFAILS.NE.0) WRITE (NOUT,FMT=*) 'SIGH..., THERE WERE ',NFAILS, + ' SUSPICIOUS TEST RESULTS.' STOP END SHAR_EOF fi # end of overwriting check if test -f 'driver8.f' then echo shar: will not over-write existing file "'driver8.f'" else cat << SHAR_EOF > 'driver8.f' * DHASRD EXAMPLE PROGRAM TEXT. * C .. Parameters .. * NIN is the standard input. * NOUT is the standard output, e.g., a terminal INTEGER NIN,NOUT PARAMETER (NIN=5,NOUT=6) INTEGER NMAX PARAMETER (NMAX=20) INTEGER LDA,LDQG,LDU PARAMETER (LDA=NMAX,LDQG=NMAX,LDU=NMAX) INTEGER LRWORK PARAMETER (LRWORK= (NMAX+NMAX)* (NMAX+NMAX+1)) DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C .. C .. Local Scalars .. INTEGER I,IERR,IJ,J,JI,N,POS,WPOS CHARACTER COMPU C .. C .. Local Arrays .. DOUBLE PRECISION A(LDA,NMAX),QG(LDQG,NMAX+1),RWORK(LRWORK), + U(LDU,NMAX) C .. C .. External Subroutines .. EXTERNAL DCOPY,DGEMM,DHASRD,DSYMV C .. * WRITE (NOUT,FMT=9000) * Skip the heading in the data file and read the data. READ (NIN,FMT='()') READ (NIN,FMT=*) N,COMPU IF (N.LE.0 .OR. N.GT.NMAX) THEN WRITE (NOUT,FMT=9010) N ELSE READ (NIN,FMT=*) ((A(I,J),J=1,N),I=1,N) READ (NIN,FMT=*) ((QG(J,I+1),I=J,N),J=1,N) READ (NIN,FMT=*) ((QG(I,J),I=J,N),J=1,N) END IF * .. square-reduce by symplectic orthogonal similarity .. CALL DHASRD(COMPU,N,A,LDA,QG,LDQG,U,LDU,RWORK,IERR) IF (IERR.NE.0) THEN WRITE (NOUT,FMT=9020) IERR ELSE * .. show the square-reduced Hamiltonian .. WRITE (NOUT,FMT=9030) DO 10 I = 1,N WRITE (NOUT,FMT=9050) (A(I,J),J=1,N), + (QG(J,I+1),J=1,I-1), (QG(I,J+1),J=I,N) 10 CONTINUE DO 20 I = 1,N WRITE (NOUT,FMT=9050) (QG(I,J),J=1,I-1), (QG(J,I),J=I,N), + (-A(J,I),J=1,N) 20 CONTINUE * .. show the square of H .. WRITE (NOUT,FMT=9040) WPOS = (NMAX+NMAX)* (NMAX+NMAX) * T * .. compute N11 = A*A + G*Q and set N22 = N11 .. CALL DGEMM('N','N',N,N,N,ONE,A,LDA,A,LDA,ZERO,RWORK,N+N) DO 40 I = 1,N CALL DCOPY(N-I+1,QG(I,I),1,RWORK(WPOS+I),1) DO 30 J = 1,I - 1 RWORK(WPOS+J) = QG(I,J) 30 CONTINUE CALL DSYMV('U',N,ONE,QG(1,2),LDQG,RWORK(WPOS+1),1,ONE, + RWORK((I-1)* (N+N)+1),1) POS = N* (N+N) + N + I CALL DCOPY(N,RWORK((I-1)* (N+N)+1),1,RWORK(POS),N+N) 40 CONTINUE DO 50 I = 1,N CALL DSYMV('U',N,-ONE,QG(1,2),LDQG,A(I,1),LDA,ZERO, + RWORK((N+I-1)* (N+N)+1),1) CALL DSYMV('L',N,ONE,QG,LDQG,A(1,I),1,ZERO, + RWORK((I-1)* (N+N)+N+1),1) 50 CONTINUE DO 70 J = 1,N DO 60 I = J,N IJ = (N+J-1)* (N+N) + I JI = (N+I-1)* (N+N) + J RWORK(IJ) = RWORK(IJ) - RWORK(JI) RWORK(JI) = -RWORK(IJ) IJ = N + I + (J-1)* (N+N) JI = N + J + (I-1)* (N+N) RWORK(IJ) = RWORK(IJ) - RWORK(JI) RWORK(JI) = -RWORK(IJ) 60 CONTINUE 70 CONTINUE DO 80 I = 1,N + N WRITE (NOUT,FMT=9050) (RWORK(I+ (J-1)* (N+N)),J=1,N+N) 80 CONTINUE END IF STOP * 9000 FORMAT (' DHASRD EXAMPLE PROGRAM RESULTS',/,1X) 9010 FORMAT (/,' N is out of range.',/,' N = ',I5) 9020 FORMAT (' IERR on exit from DHASRD = ',I2) 9030 FORMAT (/,' The square-reduced Hamiltonian is ') 9040 FORMAT (/,' The square of the square-reduced Hamiltonian is ') 9050 FORMAT (1X,8 (F10.4)) END SHAR_EOF fi # end of overwriting check if test -f 'driver9.f' then echo shar: will not over-write existing file "'driver9.f'" else cat << SHAR_EOF > 'driver9.f' * test driver for DHASRD.f * ***************************************************************** * This program sets up several ad-hoc Hamiltonian matrices and * transforms it to square-reduced form by Van Loan's method. * * ***************************************************************** * There are three consistency checks. * Test 1: Similarity Check * Check that the output Hamiltonian matrix is similar to the * input Hamiltonian matrix ``modulo rounding errors.'' * Test 2: Symplectic Orthogonality * Check that the similarity transformation is symplectic * orthogonal ``modulo rounding errors.'' * Test 3: Square Reduced Check * Check that the output Hamiltonian matrix is square reduced * ``modulo rounding errors.'' * * ****************************************************************** C .. Parameters .. * NOUT is the standard output, e.g., a terminal INTEGER NOUT PARAMETER (NOUT=6) INTEGER NMAX,NMIN PARAMETER (NMAX=15,NMIN=-1) INTEGER LD PARAMETER (LD=NMAX) INTEGER LDA,LDQG,LDU PARAMETER (LDA=LD,LDQG=LD,LDU=LD) DOUBLE PRECISION ZERO,ONE,TWO PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0) CHARACTER*(*) EQULS PARAMETER (EQULS=' =============================================') C .. C .. Local Scalars .. * DOUBLE PRECISION ANRM,EPS,GNRM,HNRM,QNRM,SIMERR,SMPERR,SQRERR,Z11, + Z12,Z21 INTEGER I,IERR,IFUDGE,ITER,J,N,NFAILS CHARACTER COMPU C .. C .. Local Arrays .. * DOUBLE PRECISION A(LDA,LDA),AP(LDA,LDA),QG(LDQG,LDQG+1), + QPGP(LDQG,LDQG+1),RWORK(2*LD),U(LDU,2*LDU), + Y(LD,LD),Z(LD,LD) INTEGER ISEED(4) C .. C .. External Functions .. * DOUBLE PRECISION DLAMCH,DLANGE,DLANSY,DLAPY2,DNRM2 LOGICAL LSAME EXTERNAL DLAMCH,DLANGE,DLANSY,DLAPY2,DNRM2,LSAME C .. C .. External Subroutines .. * EXTERNAL DCOPY,DGEMM,DHASRD,DLACPY,DLARNV,DLASET,DSYMM,DSYRK C .. C .. Intrinsic Functions .. * INTRINSIC DBLE,SQRT C .. C .. Data statements .. * . (seed for random number generation) . * * ****************************************************************** DATA ISEED/86,1967,2001,1995/ C .. EPS = DLAMCH('P') * NFAILS = 0 DO 50 ITER = 1,2 IF (ITER.EQ.1) COMPU = 'F' IF (ITER.EQ.2) COMPU = 'A' * .. 2 test runs, 1 for generating transformation matrix, 1 for * accumulating transformations .. WRITE (NOUT,FMT=9020) EQULS WRITE (NOUT,FMT=9020) EQULS IF (LSAME(COMPU,'F')) THEN WRITE (NOUT,FMT=9020) + 'TEST RUN 1: GENERATE SYMPLECTIC ORTHOGONAL' WRITE (NOUT,FMT=9020) + ' TRANSFORMATION MATRIX.' ELSE WRITE (NOUT,FMT=9020) + 'TEST RUN 2: ACCUMULATE SYMPLECTIC ORTHOGONAL ' WRITE (NOUT,FMT=9020) ' TRANSFORMATIONS.' END IF DO 40 N = NMIN,NMAX WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) EQULS WRITE (NOUT,FMT=9030) 'DHASRD TESTS: ORDER N = ',N IF (N.LT.0) THEN WRITE (NOUT,FMT=9020) + 'OF COURSE, THIS IS NOT A VALID TEST CASE.' WRITE (NOUT,FMT=9020) + 'IT CHECKS THAT DHASRD RETURNS WITH ERROR IERR = -2.' END IF * * .. set up Hamiltonian .. DO 10 I = 1,N CALL DLARNV(3,ISEED,N,A(1,I)) CALL DLARNV(3,ISEED,N,QG(1,I)) 10 CONTINUE IF(N.GE.0) THEN CALL DLARNV(3,ISEED,N,QG(1,N+1)) ENDIF CALL DLACPY('A',N,N,A,LD,AP,LD) CALL DLACPY('A',N,N+1,QG,LDQG,QPGP,LDQG) * * .. set up symplectic orthogonal matrix for accumulating * transformations .. IF (COMPU.EQ.'A') CALL DLASET('A',N,2*N,ZERO,ONE,U,LDU) * * .. Frobenius norm of Hamiltonian Matrix .. ANRM = DLANGE('F',N,N,A,LD,RWORK) GNRM = DLANSY('F','L',N,QG,LDQG,RWORK) QNRM = DLANSY('F','U',N,QG(1,2),LDQG,RWORK) HNRM = DLAPY2(DLAPY2(ANRM,GNRM),DLAPY2(ANRM,QNRM)) * * .. square reduce .. CALL DHASRD(COMPU,N,A,LDA,QG,LDQG,U,LDU,RWORK,IERR) * ****************************************************************** * .. Test 0: Has DHASRD something to tell? ===== WRITE (NOUT,FMT=9010) 'IERR = ',IERR IF (IERR.EQ.-1) WRITE (NOUT,FMT=9020) 'COMPU = '//COMPU// + ' IS NOT ONE OF ''A'', ''F'', OR ''N''.' IF (IERR.EQ.-2) WRITE (NOUT,FMT=9040) 'N = ',N, + ' IS NEGATIVE.' IF (IERR.EQ.-4) WRITE (NOUT,FMT=9050) 'LDA = ',LDA, + ' IS LESS THAN N = ',N,'.' IF (IERR.EQ.-6) WRITE (NOUT,FMT=9050) 'LDQG = ',LDQG, + ' IS LESS THAN N = ',N,'.' IF (IERR.EQ.-8) WRITE (NOUT,FMT=9050) 'LDU = ',LDU, + ' IS LESS THAN N = ',N,'.' IF (IERR.EQ.0) THEN WRITE (NOUT,FMT=9020) 'NO ERRORS WERE DETECTED.' * ****************************************************************** * .. Test 1: similarity check (modulo rounding) .. * * * .. Z11 <- || (AP*S1 - GG*S2) - (S1*A + S2*Q)|| .. CALL DGEMM('N','N',N,N,N,ONE,AP,LDA,U,LDU,ZERO,Z,LD) CALL DSYMM('L','U',N,N,-ONE,QPGP(1,2),LDQG,U(1,N+1), + LDU,ONE,Z,LD) CALL DGEMM('N','N',N,N,N,-ONE,U,LDU,A,LDA,ONE,Z,LD) CALL DSYMM('R','L',N,N,-ONE,QG,LDQG,U(1,N+1),LDU,ONE, + Z,LD) Z11 = DLANGE('F',N,N,Z,LD,RWORK) * * .. Z12 <- ||(AP*S2 + GG*S1) - (S1*G - S2*A') || .. CALL DGEMM('N','N',N,N,N,ONE,AP,LDA,U(1,N+1),LDU,ZERO, + Z,LD) CALL DSYMM('L','U',N,N,ONE,QPGP(1,2),LDQG,U,LDU,ONE,Z, + LD) CALL DSYMM('R','U',N,N,-ONE,QG(1,2),LDQG,U,LDU,ONE,Z, + LD) CALL DGEMM('N','T',N,N,N,ONE,U(1,N+1),LDU,A,LDA,ONE,Z, + LD) Z12 = DLANGE('F',N,N,Z,LD,RWORK) * * .. Z21 <- || (QQ*S1 + AP'*S2) - (-S2*A + S1*Q) || .. CALL DSYMM('L','L',N,N,ONE,QPGP,LDQG,U,LDU,ZERO,Z,LD) CALL DGEMM('T','N',N,N,N,ONE,AP,LDA,U(1,N+1),LDU,ONE, + Z,LD) CALL DGEMM('N','N',N,N,N,ONE,U(1,N+1),LDU,A,LDA,ONE,Z, + LD) CALL DSYMM('R','L',N,N,-ONE,QG,LDQG,U,LDU,ONE,Z,LD) Z21 = DLANGE('F',N,N,Z,LD,RWORK) * * .. || HH*S - S*H || / || H || .. SIMERR = DLAPY2(DLAPY2(Z11,Z12),DLAPY2(Z11,Z21)) IF (HNRM.NE.ZERO) SIMERR = SIMERR/HNRM IFUDGE = 10*N * WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) 'SIMILARITY CHECK:' WRITE (NOUT,FMT=9000) 'RELATIVE ERROR = ',SIMERR * IF (SIMERR.GT.DBLE(IFUDGE)*EPS) THEN NFAILS = NFAILS + 1 WRITE (NOUT,FMT=9020) + 'SUSPICIOUS: THE RELATIVE ERROR OF SIMILARITY' WRITE (NOUT,FMT=9020) + ' IS LARGER THAN EXPECTED. ' ELSE WRITE (NOUT,FMT=9020) + 'OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE' WRITE (NOUT,FMT=9020) + ' INPUT ''MODULO ROUNDING ERRORS'' AS IT ' WRITE (NOUT,FMT=9020) ' SHOULD BE.' END IF * ****************************************************************** * .. Test 2: symplectic orthogonal check .. * * .. Z11 <- || S2'*S1 - S1'*S2 || .. CALL DGEMM('T','N',N,N,N,ONE,U(1,N+1),LDU,U,LDU,ZERO, + Z,LD) CALL DGEMM('T','N',N,N,N,-ONE,U,LDU,U(1,N+1),LDU,ONE, + Z,LD) Z11 = DLANGE('F',N,N,Z,LD,RWORK) * * .. Z12 <- || S2'S2 + S1'S1 - I || === CALL DLASET('U',N,N,ZERO,-ONE,Z,LD) CALL DSYRK('U','T',N,N,ONE,U(1,N+1),LDU,ONE,Z,LD) CALL DSYRK('U','T',N,N,ONE,U,LDU,ONE,Z,LD) Z12 = DLANSY('F','U',N,Z,LD,RWORK) * SMPERR = SQRT(TWO)*DLAPY2(Z11,Z12) IF (N.GT.0) SMPERR = SMPERR/SQRT(DBLE(N)) IFUDGE = 10*N * WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) + 'SYMPLECTIC ORTHOGONALITY CHECK:' WRITE (NOUT,FMT=9000) 'RELATIVE ERROR = ',SMPERR IF (SMPERR.GT.DBLE(IFUDGE)*EPS) THEN NFAILS = NFAILS + 1 WRITE (NOUT,FMT=9020) + 'SUSPICIOUS: THE RELATIVE DISTANCE TO ' WRITE (NOUT,FMT=9020) + ' THE SYMPLECTIC ORTHOGONAL CLASS IS' WRITE (NOUT,FMT=9020) + ' LARGER THAN EXPECTED.' ELSE WRITE (NOUT,FMT=9020) + 'OK: THE OUTPUT SIMILARITY TRANSFORMATION IS ' WRITE (NOUT,FMT=9020) + ' SYMPLECTIC ORTHOGONAL ''MODULO ROUNDING ' WRITE (NOUT,FMT=9020) + ' ERRORS'' AS IT SHOULD BE.' END IF * ****************************************************************** * .. Test 3: square reduced check .. * * .. Z <- A*A + G*Q .. CALL DLACPY('L',N,N,QG,LDQG,Y,LD) DO 20 I = 1,N - 1 CALL DCOPY(N-I,QG(I+1,I),1,Y(I,I+1),LD) 20 CONTINUE CALL DSYMM('L','U',N,N,ONE,QG(1,2),LDQG,Y,LD,ZERO,Z, + LD) CALL DGEMM('N','N',N,N,N,ONE,A,LDA,A,LDA,ONE,Z,LD) * * .. Z11 <- FROB NORM BELOW SUB DIAGONAL .. Z11 = ZERO DO 30 J = 1,N - 2 Z11 = DLAPY2(Z11,DNRM2(N-J-1,Z(J+2,J),1)) 30 CONTINUE * * .. R21 <- || Q*A - A'*Q || .. CALL DGEMM('T','N',N,N,N,-ONE,A,LDA,Y,LD,ZERO,Z,LD) CALL DSYMM('L','L',N,N,ONE,QG,LDQG,A,LDA,ONE,Z,LD) Z21 = DLANGE('F',N,N,Z,LD,RWORK) * * .. square reduced check .. SQRERR = DLAPY2(Z11,Z21) IF (HNRM.NE.ZERO) SQRERR = SQRERR/HNRM/HNRM IFUDGE = 10*N * WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) 'SQUARE REDUCED CHECK:' WRITE (NOUT,FMT=9000) 'RELATIVE ERROR = ',SQRERR IF (SQRERR.GT.DBLE(IFUDGE)*EPS) THEN NFAILS = NFAILS + 1 WRITE (NOUT,FMT=9020) + 'SUSPICIOUS: THE RELATIVE DIFFERENCE FROM' WRITE (NOUT,FMT=9020) + ' BEING SQUARE REDUCED IS ' WRITE (NOUT,FMT=9020) + ' LARGER THAN EXPECTED. ' ELSE WRITE (NOUT,FMT=9020) + 'OK: THE OUTPUT IS SQUARE REDUCED ''MODULO ' WRITE (NOUT,FMT=9020) + ' ROUNDING ERRORS'' AS IT SHOULD BE.' END IF END IF WRITE (NOUT,FMT=9020) 40 CONTINUE * .. repeat test for accumulated transformations .. 50 CONTINUE * ****************************************************************** WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) WRITE (NOUT,FMT=9020) EQULS IF (NFAILS.EQ.0) THEN WRITE (NOUT,FMT=9020) 'NO SUSPICIOUS TEST RESULTS.' ELSE WRITE (NOUT,FMT=9040) 'SIGH..., THERE WERE ',NFAILS, + ' SUSPICIOUS TEST RESULTS.' END IF STOP * 9000 FORMAT (A,E10.4) 9010 FORMAT (A,I5) 9020 FORMAT (A) 9030 FORMAT (A,I3) 9040 FORMAT (A,I3,A) 9050 FORMAT (A,I3,A,I3,A) END SHAR_EOF fi # end of overwriting check if test -f 'res1' then echo shar: will not over-write existing file "'res1'" else cat << SHAR_EOF > 'res1' DHABL/DHASRD/DHAEVS EXAMPLE PROGRAM RESULTS The eigenvalues are 11.5734 + ( 0.0000)i 0.9905 + ( 0.0000)i 0.0000 + ( 2.6316)i 0.0000 + ( -2.6316)i -0.9905 + ( 0.0000)i -11.5734 + ( 0.0000)i SHAR_EOF fi # end of overwriting check if test -f 'res2' then echo shar: will not over-write existing file "'res2'" else cat << SHAR_EOF > 'res2' DHABL EXAMPLE PROGRAM RESULTS The scaled Hamiltonian is -0.4000 0.5000 0.7000 816.5892 33.6243 0.1849 -0.4700 0.8000 2.5000 33.6243 -12.0087 0.0096 0.0810 0.2900 -0.9000 0.1849 0.0096 0.0007 -0.0001 -0.0005 0.1790 0.4000 0.4700 -0.0810 -0.0005 0.0412 17.4874 -0.5000 -0.8000 -0.2900 0.1790 17.4874 -832.7320 -0.7000 -2.5000 0.9000 The scaling factors are 0.0020 0.0204 2.0405 SHAR_EOF fi # end of overwriting check if test -f 'res3' then echo shar: will not over-write existing file "'res3'" else cat << SHAR_EOF > 'res3' ============================================= ============================================= TESTING SYMPLECTIC SCALING ============================================= DHAEV TESTS: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHAEV RETURNS WITH AN ERROR. SYMPLECTIC SCALING ERROR: DHABL RETURNS IERR = -2 ============================================= DHAEV TESTS: ORDER N = 0 SYMPLECTIC SCALING NO ERRORS DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 1 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 0.214824703705671 QPNRM = 0.214824703705671 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 2 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 2.39383677704275 QPNRM = 2.39383677704275 SIMILARITY CHECK: RELATIVE ERROR = 2.71358172647186D-016 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 3 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 15.2256806281845 QPNRM = 15.2256806281845 SIMILARITY CHECK: RELATIVE ERROR = 2.48715480617792D-016 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 4 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 18.6106681314024 QPNRM = 18.6106681314024 SIMILARITY CHECK: RELATIVE ERROR = 1.05324995919445D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 5 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 20.0403972386819 QPNRM = 20.0403972386819 SIMILARITY CHECK: RELATIVE ERROR = 8.98671786199432D-016 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 6 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 27.0871508893049 QPNRM = 27.0871508893049 SIMILARITY CHECK: RELATIVE ERROR = 2.02603694132078D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 7 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 31.9715120365886 QPNRM = 31.9715120365887 SIMILARITY CHECK: RELATIVE ERROR = 1.80862850391281D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 8 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 882.676073201900 QPNRM = 882.676073201900 SIMILARITY CHECK: RELATIVE ERROR = 1.65131705626505D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 9 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 439.327533212684 QPNRM = 439.327533212684 SIMILARITY CHECK: RELATIVE ERROR = 2.71075777868489D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 10 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 1309.25463730288 QPNRM = 1309.25463730288 SIMILARITY CHECK: RELATIVE ERROR = 1.88490594202794D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 11 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 1741.05138887116 QPNRM = 1741.05138887116 SIMILARITY CHECK: RELATIVE ERROR = 2.50411707039103D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 12 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 3176.84165062155 QPNRM = 3176.84165062155 SIMILARITY CHECK: RELATIVE ERROR = 2.28384097451001D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 13 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 15116.5791935657 QPNRM = 15116.5791935657 SIMILARITY CHECK: RELATIVE ERROR = 3.29590245254523D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 14 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 3473.23144926256 QPNRM = 3473.23144926256 SIMILARITY CHECK: RELATIVE ERROR = 3.90823874717891D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 15 SYMPLECTIC SCALING NO ERRORS DETECTED. SCALE OK: GPNRM = 23047.5564527910 QPNRM = 23047.5564527910 SIMILARITY CHECK: RELATIVE ERROR = 3.74980345485544D-015 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= ============================================= ============================================= DHAEV TESTS: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHAEV RETURNS WITH AN ERROR. NORM SCALING ERROR: DHABL RETURNS IERR = -2 ============================================= DHAEV TESTS: ORDER N = 0 NORM SCALING NO ERRORS DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 1 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 8.00000000000000 ANRM = 0.112480187217699 GNRM = 4.72836791261691D-002 QNRM = 7.15717213448476 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 2 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 32.0000000000000 ANRM = 2.08543567037179 GNRM = 7.35338166707300D-002 QNRM = 38.6716641328500 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 3 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 128.000000000000 ANRM = 6.92985441365365 GNRM = 8.27674077762709D-002 QNRM = 111.602234394339 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 4 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 128.000000000000 ANRM = 13.0185885521640 GNRM = 8.14982239347055D-002 QNRM = 96.6854420214134 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 5 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 256.000000000000 ANRM = 26.4153828934938 GNRM = 0.203519877652393 QNRM = 275.083577477601 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 6 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 512.000000000000 ANRM = 44.7518865401816 GNRM = 0.127302023790044 QNRM = 495.648625058800 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 7 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 512.000000000000 ANRM = 129.194114084613 GNRM = 8.86817468456751D-002 QNRM = 544.276163281934 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 8 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 512.000000000000 ANRM = 111.675056434900 GNRM = 0.102327196073017 QNRM = 701.419921152902 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 9 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 1024.00000000000 ANRM = 383.450623043022 GNRM = 9.51123274416865D-002 QNRM = 806.754300003111 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 10 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 1024.00000000000 ANRM = 453.189124764968 GNRM = 9.46373237162815D-002 QNRM = 1274.50458407210 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 11 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 2048.00000000000 ANRM = 3029.05979356483 GNRM = 9.41943554042665D-002 QNRM = 1462.80363077944 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 12 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 2048.00000000000 ANRM = 2436.91067122740 GNRM = 0.113971094531401 QNRM = 1399.15367777839 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 13 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 8192.00000000000 ANRM = 7090.34171108275 GNRM = 0.108161346046478 QNRM = 1905.00516147236 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 14 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 16384.0000000000 ANRM = 16078.4411578713 GNRM = 0.112534080495865 QNRM = 1966.31749866508 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHAEV TESTS: ORDER N = 15 NORM SCALING NO ERRORS DETECTED. SCALE OK: TAU = 32768.0000000000 ANRM = 32113.0244594608 GNRM = 0.112954504364062 QNRM = 2178.85440466470 SIMILARITY CHECK: RELATIVE ERROR = 0.00000000000000D+000 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= ============================================= ============================================= NO SUSPICIOUS TEST RESULTS. SHAR_EOF fi # end of overwriting check if test -f 'res4' then echo shar: will not over-write existing file "'res4'" else cat << SHAR_EOF > 'res4' DHAEVD EXAMPLE PROGRAM RESULTS The eigenvalues are 11.5734 + ( 0.0000)i 0.9905 + ( 0.0000)i 0.0000 + ( 2.6316)i 0.0000 + ( -2.6316)i -0.9905 + ( 0.0000)i -11.5734 + ( 0.0000)i SHAR_EOF fi # end of overwriting check if test -f 'res5' then echo shar: will not over-write existing file "'res5'" else cat << SHAR_EOF > 'res5' ============================================= ============================================= ============================================= DHAEVD TEST: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHAEVD RETURNS WITH ERROR IERR = -1. DHAEVD RETURNS IERR = -1 N = -1 IS NEGATIVE. ============================================= DHAEVD TEST: ORDER N = 0 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 0.00000000000000D+000 ============================================= DHAEVD TEST: ORDER N = 1 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 3.04914049851180D-017 ============================================= DHAEVD TEST: ORDER N = 2 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.32626163494360D-017 ============================================= DHAEVD TEST: ORDER N = 3 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.49766668244639D-016 ============================================= DHAEVD TEST: ORDER N = 4 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 3.20153267491709D-016 ============================================= DHAEVD TEST: ORDER N = 5 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.34089298326831D-016 ============================================= DHAEVD TEST: ORDER N = 6 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.26900142873816D-016 ============================================= DHAEVD TEST: ORDER N = 7 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.94648809454923D-017 ============================================= DHAEVD TEST: ORDER N = 8 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 2.04240075139212D-016 ============================================= DHAEVD TEST: ORDER N = 9 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.57456654601777D-017 ============================================= DHAEVD TEST: ORDER N = 10 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.46328999721847D-017 ============================================= DHAEVD TEST: ORDER N = 11 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.44489843930722D-016 ============================================= DHAEVD TEST: ORDER N = 12 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.30650717808721D-017 ============================================= DHAEVD TEST: ORDER N = 13 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.42408777704300D-017 ============================================= DHAEVD TEST: ORDER N = 14 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 5.44672709985278D-017 ============================================= DHAEVD TEST: ORDER N = 15 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.26263144193693D-017 ============================================= NO SUSPICIOUS TEST RESULTS. SHAR_EOF fi # end of overwriting check if test -f 'res6' then echo shar: will not over-write existing file "'res6'" else cat << SHAR_EOF > 'res6' DHAEVS EXAMPLE PROGRAM RESULTS The eigenvalues are 2.0000 + ( 1.0000)i 2.0000 + ( -1.0000)i 1.4142 + ( 0.0000)i -1.4142 + ( 0.0000)i -2.0000 + ( 1.0000)i -2.0000 + ( -1.0000)i SHAR_EOF fi # end of overwriting check if test -f 'res7' then echo shar: will not over-write existing file "'res7'" else cat << SHAR_EOF > 'res7' ============================================= ============================================= TEST RUN WITHOUT SCALING ============================================= DHAEVS TEST: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHAEVS RETURNS WITH ERROR IERR = -2. DHASRD RETURNS IERR = -2 N = -1 IS NEGATIVE. DHAEVS RETURNS IERR = -2 N = -1 IS NEGATIVE. ============================================= DHAEVS TEST: ORDER N = 0 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 0.00000000000000D+000 ============================================= DHAEVS TEST: ORDER N = 1 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 3.04914049851180D-017 ============================================= DHAEVS TEST: ORDER N = 2 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.32626163494360D-017 ============================================= DHAEVS TEST: ORDER N = 3 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.49766668244639D-016 ============================================= DHAEVS TEST: ORDER N = 4 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 3.20153267491709D-016 ============================================= DHAEVS TEST: ORDER N = 5 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.34089298326831D-016 ============================================= DHAEVS TEST: ORDER N = 6 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.29877393338627D-016 ============================================= DHAEVS TEST: ORDER N = 7 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.94648809454923D-017 ============================================= DHAEVS TEST: ORDER N = 8 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 2.04240075139212D-016 ============================================= DHAEVS TEST: ORDER N = 9 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.57456654601777D-017 ============================================= DHAEVS TEST: ORDER N = 10 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.82589212405420D-017 ============================================= DHAEVS TEST: ORDER N = 11 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.42381706309488D-017 ============================================= DHAEVS TEST: ORDER N = 12 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.30650717808721D-017 ============================================= DHAEVS TEST: ORDER N = 13 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 5.05376445522604D-017 ============================================= DHAEVS TEST: ORDER N = 14 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 8.81927893694238D-017 ============================================= DHAEVS TEST: ORDER N = 15 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 2.33923115848030D-017 ============================================= ============================================= TEST RUN WITH SCALING ============================================= DHAEVS TEST: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHAEVS RETURNS WITH ERROR IERR = -2. DHASRD RETURNS IERR = -2 N = -1 IS NEGATIVE. DHAEVS RETURNS IERR = -2 N = -1 IS NEGATIVE. ============================================= DHAEVS TEST: ORDER N = 0 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 0.00000000000000D+000 ============================================= DHAEVS TEST: ORDER N = 1 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.81321850595529D-016 ============================================= DHAEVS TEST: ORDER N = 2 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.27168745252820D-016 ============================================= DHAEVS TEST: ORDER N = 3 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 2.22750817839991D-016 ============================================= DHAEVS TEST: ORDER N = 4 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 9.62035219694833D-017 ============================================= DHAEVS TEST: ORDER N = 5 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.52542050682234D-017 ============================================= DHAEVS TEST: ORDER N = 6 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.26319855291347D-016 ============================================= DHAEVS TEST: ORDER N = 7 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 8.30976822568448D-017 ============================================= DHAEVS TEST: ORDER N = 8 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 1.21915819356464D-016 ============================================= DHAEVS TEST: ORDER N = 9 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 7.71183068963740D-017 ============================================= DHAEVS TEST: ORDER N = 10 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 3.48599858950901D-016 ============================================= DHAEVS TEST: ORDER N = 11 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 9.61460310186976D-017 ============================================= DHAEVS TEST: ORDER N = 12 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 9.37898461647893D-017 ============================================= DHAEVS TEST: ORDER N = 13 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 9.46424684079837D-017 ============================================= DHAEVS TEST: ORDER N = 14 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 6.03952805104437D-017 ============================================= DHAEVS TEST: ORDER N = 15 NO ERRORS WERE DETECTED. MAXIMUM RATIO = 4.25285054098332D-017 ============================================= NO SUSPICIOUS TEST RESULTS. SHAR_EOF fi # end of overwriting check if test -f 'res8' then echo shar: will not over-write existing file "'res8'" else cat << SHAR_EOF > 'res8' DHASRD EXAMPLE PROGRAM RESULTS The square-reduced Hamiltonian is 1.0000 3.3485 0.3436 1.0000 1.9126 -0.1072 6.7566 11.0750 -0.3014 1.9126 8.4479 -1.0790 2.3478 1.6899 -2.3868 -0.1072 -1.0790 -2.9871 7.0000 8.6275 -0.6352 -1.0000 -6.7566 -2.3478 8.6275 16.2238 -0.1403 -3.3485 -11.0750 -1.6899 -0.6352 -0.1403 1.2371 -0.3436 0.3014 2.3868 The square of the square-reduced Hamiltonian is 48.0000 80.6858 -2.5217 0.0000 1.8590 -10.5824 167.8362 298.4815 -4.0310 -1.8590 0.0000 -33.1160 0.0000 4.5325 2.5185 10.5824 33.1160 0.0000 0.0000 0.0000 0.0000 48.0000 167.8362 0.0000 0.0000 0.0000 0.0000 80.6858 298.4815 4.5325 0.0000 0.0000 0.0000 -2.5217 -4.0310 2.5185 SHAR_EOF fi # end of overwriting check if test -f 'res9' then echo shar: will not over-write existing file "'res9'" else cat << SHAR_EOF > 'res9' ============================================= ============================================= TEST RUN 1: GENERATE SYMPLECTIC ORTHOGONAL TRANSFORMATION MATRIX. ============================================= DHASRD TESTS: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHASRD RETURNS WITH ERROR IERR = -2. IERR = -2 N = -1 IS NEGATIVE. ============================================= DHASRD TESTS: ORDER N = 0 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 1 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 2 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.9309E-16 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.6852E-16 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.7065E-18 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 3 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.5776E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.7567E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5012E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 4 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.4201E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8232E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6567E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 5 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.4592E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.4234E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5775E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 6 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.5469E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.6165E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6139E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 7 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7371E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.7638E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6747E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 8 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.6731E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8331E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.4146E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 9 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.5766E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8345E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.3707E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 10 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.8684E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8681E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6717E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 11 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7600E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.1039E-14 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6265E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 12 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.6507E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.6864E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.4831E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 13 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7996E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.9113E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.7316E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 14 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.8700E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.1001E-14 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5954E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 15 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7674E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8019E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5804E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= ============================================= TEST RUN 2: ACCUMULATE SYMPLECTIC ORTHOGONAL TRANSFORMATIONS. ============================================= DHASRD TESTS: ORDER N = -1 OF COURSE, THIS IS NOT A VALID TEST CASE. IT CHECKS THAT DHASRD RETURNS WITH ERROR IERR = -2. IERR = -2 N = -1 IS NEGATIVE. ============================================= DHASRD TESTS: ORDER N = 0 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 1 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 2 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.3836E-16 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.0000E+00 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.2687E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 3 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.5912E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.5425E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5281E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 4 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.4659E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.7732E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.7809E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 5 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7273E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8631E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.2484E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 6 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.5056E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.5897E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.3751E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 7 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.6284E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8475E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.8228E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 8 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.6185E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.8405E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.3792E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 9 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.6004E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.7230E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.7241E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 10 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7288E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.1014E-14 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.4731E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 11 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7396E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.9724E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5115E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 12 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7823E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.9687E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5684E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 13 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.9421E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.1279E-14 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5695E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 14 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.7472E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.9712E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.5340E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= DHASRD TESTS: ORDER N = 15 IERR = 0 NO ERRORS WERE DETECTED. SIMILARITY CHECK: RELATIVE ERROR = 0.8658E-15 OK: THE OUTPUT APPEARS TO BE SIMILAR TO THE INPUT 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SYMPLECTIC ORTHOGONALITY CHECK: RELATIVE ERROR = 0.9001E-15 OK: THE OUTPUT SIMILARITY TRANSFORMATION IS SYMPLECTIC ORTHOGONAL 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. SQUARE REDUCED CHECK: RELATIVE ERROR = 0.6217E-16 OK: THE OUTPUT IS SQUARE REDUCED 'MODULO ROUNDING ERRORS' AS IT SHOULD BE. ============================================= NO SUSPICIOUS TEST RESULTS. SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' SUBROUTINE DCROOT(XR,XI,YR,YI) C C PURPOSE C C To compute the complex square root YR + i*YI of a complex number C XR + i*XI in real arithmetic. The branch is chosen so that C YR .GE. 0.0 and SIGN(YI) .EQ. SIGN(XI). C C ARGUMENTS C C XR (input) DOUBLE PRECISION C XI (input) DOUBLE PRECISION C These scalars define the real and imaginary part of the C complex number of which the root is sought. C C YR (input) DOUBLE PRECISION C YI (input) DOUBLE PRECISION C These scalars define the real and imaginary part of the C complex root. C C REFERENCE C C Adapted from EISPACK subroutine CSROOT. C C CONTRIBUTOR C C P. Benner, Universitaet Bremen, Germany, and C R. Byers, University of Kansas, Lawrence, USA. C C REVISIONS C C 1993, November 15. C 1998, August 26. C C*********************************************************************** C C C .. Scalar Arguments .. DOUBLE PRECISION XI,XR,YI,YR C .. C .. Parameters .. C DOUBLE PRECISION ZERO,HALF PARAMETER (ZERO=0.0D0,HALF=1.0D0/2.0D0) C .. C .. Local Scalars .. C DOUBLE PRECISION S C .. C .. External Subroutines .. C . LAPACK . C C .. C .. Intrinsic Functions .. C C*********************************************************************** INTRINSIC ABS,SQRT C .. C .. External Functions .. DOUBLE PRECISION DLAPY2 EXTERNAL DLAPY2 C .. C S = SQRT(HALF* (DLAPY2(XR,XI)+ABS(XR))) IF (XR.GE.ZERO) YR = S IF (XI.LT.ZERO) S = -S IF (XR.LE.ZERO) YI = S IF (XR.LT.ZERO) YR = HALF* (XI/YI) IF (XR.GT.ZERO) YI = HALF* (XI/YR) C RETURN END SUBROUTINE DHABL(JOBSCL,N,A,LDA,QG,LDQG,D,RWORK,IERR) C C .. Scalar Arguments .. INTEGER IERR,LDA,LDQG,N CHARACTER JOBSCL C .. C .. Array Arguments .. C C PURPOSE C C Perform a symplectic scaling on the Hamiltonian matrix C C ( A G ) C (*) H = ( T ), C ( Q -A ) C C i.e., perform either the symplectic scaling transformation C C -1 C ( A' G' ) ( D 0 ) ( A G ) ( D 0 ) C (**) H' <-- ( T ) = ( ) ( T) ( -1 ) C ( Q' -A' ) ( 0 D ) ( Q -A ) ( 0 D ) C C where D is a diagonal scaling matrix, or the symplectic norm C scaling transformation C C ( A'' G'' ) 1 ( A G/tau ) C (***) H'' <-- ( T ) = --- ( T ) C ( Q'' -A'' ) tau (tau Q -A ) C C where tau is a real scalar. Note that if tau is not equal to 1 C then (***) is NOT a similarity transformation. The eigenvalues C of H are then tau times the eigenvalues of H''. C C For symplectic scaling (**), D is chosen to give the rows and C columns of A' approximately equal 1-norms and to give Q' and G' C approximately equal norms. (See METHOD below for details.) For C norm scaling, tau = MAX(1, ||A||, ||G||, ||Q||) where ||.|| C denotes the 1-norm (column sum norm). C C ARGUMENTS C C Mode Parameters C C JOBSCL CHARACTER*1 C Indicates which scaling strategy is used. C = 'A' or 'a': do the symplectic scaling (**); C = 'B' or 'b': do the norm scaling (***); C = 'N' or 'n': do nothing; set IERR = 0 and return. C C Input/Output Parameters C C N (input) INTEGER C The order of the matrices A, G, and Q. Normally, N.GE.1. C If N .EQ. 0, then IERR is set to zero and DHABL returns C without referencing any other argument. C C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) C On input the leading N-by-N part of this array contains C the upper left block A of the Hamiltonian matrix H in C (*). C On output A has been overwritten by the leading N-by-N C part of the scaled Hamiltonian matrix H' in (**) or H'' in C (***) depending on the setting of JOBSCL. C C LDA (input) INTEGER C The leading dimension of array A as declared in the C calling program. LDA .GE. N. C C QG (input/output) DOUBLE PRECISION array, C dimension (LDQG,N+1) C On input QG contains the upper right symmetric block G and C the lower left symmetric block Q of the Hamiltonian matrix C in (*). On output these are overwritten by the C corresponding blocks of the scaled Hamiltonian matrix H' C in (**) or H'' in (***) depending on the setting of C JOBSCL. The diagonal and lower triangle of the submatrix C in columns 1 to N contain the lower triangular part of Q. C The diagonal and upper triangle of the submatrix in C columns 2 to N+1 contain the upper triangle of G. C So, if I. GE. J, then Q(i,j) is stored in QG(i,j) C and G(i,j) is in QG(i,j+1). C C LDQG (input) INTEGER C The leading dimension of array QG just as declared in the C calling procedure. LDQG .GE. N. C C D (output) DOUBLE PRECISION array, dimension nd C If JOBSCL = 'A' or 'a', then nd = N and on output, D C contains the diagonal elements of the diagonal scaling C matrix in (**). C If JOBSCL = 'B' or 'b', then nd = 1 and D(1) is set to tau C from (***). In this case, no other elements of D are C referenced. C C Workspace C C RWORK DOUBLE PRECISION array, dimension at least N C C Error Indicator C C IERR INTEGER C = 0: successful exit; C < 0: if IERR = -j, then the j-th argument had an illegal C value; C C METHOD C C 1. Symplectic scaling (JOBSCL = 'A' or 'a'): C First, LAPACK subroutine DGEBAL is used to equilibrate the C 1-norms of the rows and columns of A using a diagonal scaling C matrix D_A. Then, H is similarily transformed by the symplectic C diagonal matrix D1 = diag(D_A,D_A**(-1)). Next, the C off-diagonal blocks of the resulting Hamiltonian matrix are C equilibrated in the 1-norm using the symplectic diagonal matrix C D2 of the form C C ( I/rho 0 ) C D2 = ( ) C ( I rho*I ) C C where rho is a real scalar. Thus, in (**), D = D1*D2. C C 2. Norm scaling (JOBSCL = 'B' or 'b'): C The norm of the matrices A and G of (*) is reduced by setting C A := A/tau and G := G/(tau**2) where tau is the power of the C base of the arithmetic closest to MAX(1, ||A||, ||G||, ||Q||) C and ||.|| denotes the 1-norm. C C NUMERICAL ASPECTS C C For symplectic scaling, the complexity of the used algorithms is C hard to estimate and depends upon how well the rows and columns of C A in (*) are equilibrated. In one sweep, each row/column of A is C scaled once, i.e., the cost of one sweep is N**2 multiplications. C Usually, 3-6 sweeps are enough to equilibrate the norms of the C rows and columns of a matrix. Roundoff errors are possible as C LAPACK routine DGEBAL does NOT use powers of the machine base for C scaling. The second stage (equilibrating ||G|| and ||Q||) requires C N**2 multiplications. C For norm scaling, 3*N**2 + O(N) multiplications are required and C NO rounding errors occur as all multiplications are performed with C powers of the machine base. C C CONTRIBUTOR C C P. Benner, Universitaet Bremen, Germany, and C R. Byers, University of Kansas, Lawrence, USA. C C REVISIONS C C 1997, December 5. C 1998, August 25. C C KEYWORDS C C Hamiltonian matrix, symplectic similarity transformation, C scaling. C C ************************************************************* * DOUBLE PRECISION A(LDA,N),D(*),QG(LDQG,N+1),RWORK(N) C .. C .. Parameters .. * DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C .. C .. Local Scalars .. * DOUBLE PRECISION ANRM,BASE,EPS,GNRM,OFL,QNRM,RHO,SFMAX,SFMIN,TAU, + UFL,Y INTEGER I,IHI,ILO,INFO,J LOGICAL NONE,NORM,SYMP C .. C .. External Functions .. * . BLAS, LAPACK . * DOUBLE PRECISION DLAMCH,DLANGE,DLANSY LOGICAL LSAME EXTERNAL DLAMCH,DLANGE,DLANSY,LSAME C .. C .. External Subroutines .. * . BLAS, LAPACK . * . other . * EXTERNAL DGEBAL,DLABAD,DLASCL,DRSCL C .. C .. Intrinsic Functions .. * INTRINSIC ABS,MAX,SQRT C .. * IERR = 0 IF (N.NE.0) THEN SYMP = LSAME(JOBSCL,'A') NORM = LSAME(JOBSCL,'B') NONE = LSAME(JOBSCL,'N') * IF (.NOT.SYMP .AND. .NOT.NORM .AND. .NOT.NONE) IERR = -1 IF (N.LT.0) IERR = -2 IF (LDA.LT.N) IERR = -4 IF (LDQG.LT.N) IERR = -6 * IF (IERR.EQ.0) THEN * * .. some machine dependant constants .. BASE = DLAMCH('B') EPS = DLAMCH('P') UFL = DLAMCH('S') OFL = ONE/UFL CALL DLABAD(UFL,OFL) SFMAX = (EPS/BASE)/UFL SFMIN = ONE/SFMAX * .. set tolerance TOL .. IF (NORM) THEN ANRM = DLANGE('1',N,N,A,LDA,RWORK) GNRM = DLANSY('1','U',N,QG(1,2),LDQG,RWORK) QNRM = DLANSY('1','L',N,QG,LDQG,RWORK) Y = MAX(ONE,ANRM,GNRM,QNRM) TAU = ONE * 10 CONTINUE * .. while TAU < Y .. IF ((TAU.LT.Y) .AND. (TAU.LT.SFMAX)) THEN TAU = TAU*BASE GO TO 10 END IF * .. end while .. IF (TAU.GT.ONE) THEN IF (ABS(TAU/BASE-Y).LT.ABS(TAU-Y)) TAU = TAU/BASE CALL DLASCL('G',N,N,TAU,ONE,N,N,A,LDA,INFO) CALL DLASCL('U',N,N,TAU,ONE,N,N,QG(1,2),LDQG,INFO) CALL DLASCL('U',N,N,TAU,ONE,N,N,QG(1,2),LDQG,INFO) END IF * D(1) = TAU * ELSE IF (SYMP) THEN CALL DGEBAL('S',N,A,LDA,ILO,IHI,D,INFO) DO 30 J = 1,N DO 20 I = J,N QG(I,J) = QG(I,J)*D(J)*D(I) QG(J,I+1) = QG(J,I+1)/D(J)/D(I) 20 CONTINUE 30 CONTINUE * GNRM = DLANSY('1','U',N,QG(1,2),LDQG,RWORK) QNRM = DLANSY('1','L',N,QG,LDQG,RWORK) IF (GNRM.EQ.ZERO) THEN IF (QNRM.EQ.ZERO) RHO = ONE IF (QNRM.NE.ZERO) RHO = SFMAX ELSE IF (QNRM.EQ.ZERO) THEN IF (GNRM.EQ.ZERO) RHO = ONE IF (GNRM.NE.ZERO) RHO = SFMIN ELSE RHO = SQRT(QNRM)/SQRT(GNRM) END IF * CALL DLASCL('L',0,0,RHO,ONE,N,N,QG,LDQG,INFO) CALL DLASCL('U',0,0,ONE,RHO,N,N,QG(1,2),LDQG,INFO) CALL DRSCL(N,SQRT(RHO),D,1) * END IF * END IF * END IF * *** Last Line of DHABL *** END SUBROUTINE DHAEVD(N,A,LDA,QG,LDQG,WR,WI,RWORK,IERR) C C PURPOSE C C This is an ``easy-to-use'' driver that uses the method of [1] to C calculate the eigenvalues of a Hamiltonian matrix. A matrix H C is Hamiltonian if it has the form of C C ( A G ) T T C (*) H = ( T), where G = G , Q = Q . C ( Q -A ) C C C ARGUMENTS C C Input/Output Parameters C C N (input) INTEGER C The order of the matrices A, G, and Q in (*). Ordinarily, C N must be positive. If N equals zero, then DHAEVD sets C IERR = 0 and returns immediately without accessing the C other arguments. C C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) C On input, the leading N-by-N part of this array contains C the upper left block A of the Hamiltonian matrix H in (*). C On output A is overwritten by the leading N-by-N upper C left block of the square reduced Hamiltonian H' in (**). C (See Section METHOD below.) C C LDA (input) INTEGER C The leading dimension of array A as declared in the C calling program. LDA .GE. N. C C QG (input/output) DOUBLE PRECISION array, C dimension (LDQG,N+1) C On input QG contains the upper right symmetric block G and C the lower left symmetric block Q of the Hamiltonian H in C (*). C On output QG contains the upper right symmetric block G' C and lower left symmetric block Q' of the square reduced C Hamiltonian matrix H' (**) as produced by DHASRD; see C Section METHOD below. C The diagonal and lower triangle of the submatrix in C columns 1 to N contain the lower triangular part of Q. C The diagonal and upper triangle of the submatrix in C columns 2 to N+1 contain the upper triangle of G. C So, if I. GE. J, then Q(i,j) = Q(j,i) is stored in QG(i,j) C and G(i,j) = G(j,i) is in QG(j,i+1). C C LDQG (input) INTEGER C The leading dimension of array QG as declared in the C calling procedure. LDQG .GE. N. C C WR (output) DOUBLE PRECISION array, dimension at least N C WI (output) DOUBLE PRECISION array, dimension at least N. C WR and WI are overwritten by N eigenvalues of H with C non-negative real part. The remaining N eigenvalues are C the negatives of these eigenvalues. The real parts C overwrite WR. The imaginary parts overwrite WI. C Eigenvalues are stored in WR and WI in decreasing order of C magnitude of the real parts, i.e., WR(I) .GE. WR(I+1). C (In particular, an eigenvalue closest to the imaginary C axis is WR(N)+WI(N)i.) C In addition, eigenvalues with zero real part are sorted in C decreasing order of magnitude of imaginary parts. Note C that non-real eigenvalues with non-zero real part appear C in complex conjugate pairs, but eigenvalues with zero real C part do not, in general, appear in complex conjugate C pairs. C C Workspace C C RWORK DOUBLE PRECISION array, dimension at least N*(N+1) C C Error Indicator C C IERR INTEGER C = 0: successful exit; C < 0: if IERR = -j, then the j-th argument had an illegal C value; C > 0: if IERR = j, then LAPACK subroutine DHSEQR failed C to converge while computing the j-th eigenvalue. C C METHOD C C DHAEVD calls DHASRD to transform the Hamiltonian matrix H in (*) C into a square-reduced Hamiltonian matrix C C (A' G' ) C (**) H' = ( T) C (Q' -A' ) C T C by an orthogonal symplectic similarity transformation H' = U H U C where C ( U1 U2 ) C (***) U = ( ). C ( -U2 U1 ) C T C The square reduced Hamiltonian matrix satisfies Q'A' - A' Q' = 0, C and C C 2 T 2 ( A'' G'') C H' := (U H U) = ( T ). C ( 0 A'') C C In addition, A'' is upper Hessenberg and G'' is skew symmetric. C The square roots of the eigenvalues of A'' = A'A' + G'Q' are the C eigenvalues of H. C DHAEVD calls DHAEVS to form A'', to compute its eigenvalues, and C to recover the eigenvalues of H from those of A''. C C REFERENCES C C [1] Van Loan, C. F. C A symplectic method for approximating all the eigenvalues of C a Hamiltonian matrix. C Linear Algebra and its Applications 61 (1984), pp. 233-251. C C [2] Byers, R. C Hamiltonian and symplectic algorithms for the algebraic C Riccati equation. C Ph. D. Thesis, Cornell University, Ithaca, NY, January 1983. C C [3] Benner, P., Byers, R., and Barth, E. C Fortran 77 Subroutines for Computing the Eigenvalues of C Hamiltonian Matrices {I}: The Square-Reduced Method. C Submitted to ACM Trans. Math. Software, 1998. C C NUMERICAL ASPECTS C C The algorithm requires (86/3)*N**3 + O(N**2) floating point C operations. C Eigenvalues computed by this subroutine are exact eigenvalues C of a perturbed Hamiltonian matrix H + E where C C || E || <= c sqrt(eps) || H ||, C C c is a modest constant depending on the dimension N and eps is the C machine precision. Moreover, if the norm of H and an eigenvalue C are of roughly the same magnitude, the computed eigenvalue is C essentially as accurate as the computed eigenvalue obtained by C traditional methods. See [1] or [2]. C C CONTRIBUTOR C C P. Benner, Universitaet Bremen, Germany, and C R. Byers, University of Kansas, Lawrence, USA. C C REVISIONS C C 1997, December 5. C 1998, August 25. C C KEYWORDS C C Eigenvalues, (square-reduced) Hamiltonian matrix. C C ****************************************************************** C C .. Scalar Arguments .. INTEGER IERR,LDA,LDQG,N C .. C .. Array Arguments .. DOUBLE PRECISION A(LDA,N),QG(LDQG,N+1),RWORK(N* (N+1)),WI(N),WR(N) C .. C .. Local Scalars .. C .. C .. External Subroutines .. EXTERNAL DHAEVS,DHASRD C .. C .. Local Arrays .. DOUBLE PRECISION DUMMY(1) C .. * IERR = 0 IF (N.NE.0) THEN * IF (N.LT.0) IERR = -1 IF (LDA.LT.N) IERR = -3 IF (LDQG.LT.N) IERR = -5 * IF (IERR.EQ.0) CALL DHASRD('N',N,A,LDA,QG,LDQG,DUMMY,1,RWORK, + IERR) IF (IERR.EQ.0) CALL DHAEVS('S',N,A,LDA,QG,LDQG,WR,WI,RWORK, + IERR) * END IF * RETURN * *** Last Line of DHAEVD *** END SUBROUTINE DHAEVS(JOBSCL,N,A,LDA,QG,LDQG,WR,WI,RWORK,IERR) C C .. Scalar Arguments .. INTEGER IERR,LDA,LDQG,N CHARACTER JOBSCL C .. C .. Array Arguments .. C C PURPOSE C C This subroutine computes the eigenvalues of an N-by-N square- C reduced Hamiltonian matrix C C ( A G ) C (*) H = ( T ) C ( Q -A ) C C Here, A is an N-by-N matrix, and G and Q are symmetric N-by-N C matrices. It is assumed without a check that H is square-reduced, C i.e., that C C 2 ( A' G' ) C (**) H = ( T ) with A' upper Hessenberg. C ( 0 A' ) C C T 2 C (Equivalently, Q'A'- A' Q' = 0, A'' = A' + G'Q', and for i > j+1, C A''(i,j) = 0.) Ordinarily, H is the output from DHASRD. C The eigenvalues of H are computed as the square roots of the C eigenvalues of A''. C C ARGUMENTS C C Mode Parameters C C JOBSCL CHARACTER*1 C JOBSCL is the mode parameter passed to LAPACK subroutine C DGEBAL to optionally permute and balance the Hessenberg C matrix A' in (**). See LAPACK subroutine DGEBAL and C Section METHOD below. C = 'N' or 'n': do neither scaling nor balancing; C = 'P' or 'p': attempt to permute the A' in (**) to block C triangular form; C = 'S' or 's': do scaling in order to equilibrate the rows C and columns of A'; C = 'B' or 'b': perform both, 'P' and 'S'. C C Input/Output Parameters C C N (input) INTEGER C The order of the matrices A, G, and Q. Ordinarily N C must be positive. If N equals zero, then DHAEVS sets C IERR = 0 and returns immediately without accessing the C other arguments. C C A (input) DOUBLE PRECISION array, dimension (LDA,N) C The leading N-by-N part of this array contains the C upper left block A of the square-reduced Hamiltonian C matrix H in (*) as produced by DHASRD. C C LDA (input) INTEGER C The leading dimension of array A as declared in the C calling program. LDA .GE. N. C C QG (input) DOUBLE PRECISION array, dimension (LDQG,N+1) C QG contains the upper right symmetric block G and the C lower left symmetric block Q of the square-reduced C Hamiltonian matrix H in (*) as produced by DHASRD. C The diagonal and lower triangle of the submatrix in C columns 1 to N contain the lower triangular part of Q. C The diagonal and upper triangle of the submatrix in C columns 2 to N+1 contain the upper triangle of G. C So, if I. GE. J, then Q(i,j) = Q(j,i) is stored in QG(i,j) C and G(i,j) = G(j,i) is in QG(j,i+1). C C LDQG (input) INTEGER C The leading dimension of array QG just as declared in the C calling procedure. LDQG. GE. N. C C WR (output) DOUBLE PRECISION array, dimension at least N. C WI (output) DOUBLE PRECISION array, dimension at least N. C WR and WI are overwritten by N eigenvalues of H with C non-negative real part. The remaining N eigenvalues are C the negatives of these eigenvalues. The real parts C overwrite WR. The imaginary parts overwrite WI. C Eigenvalues are stored in WR and WI in decreasing order of C magnitude of the real parts, i.e., WR(I) .GE. WR(I+1). C (In particular, an eigenvalue closest to the imaginary C axis is WR(N)+WI(N)i.) C In addition, eigenvalues with zero real part are sorted in C decreasing order of magnitude of imaginary parts. Note C that non-real eigenvalues with non-zero real part appear C in complex conjugate pairs, but eigenvalues with zero real C part do not, in general, appear in complex conjugate C pairs. C C Workspace C C RWORK DOUBLE PRECISION array, dimension at least N*(N+1) C C Error Indicator C C IERR INTEGER C = 0: successful exit; C < 0: if IERR = -j, then the j-th argument had an illegal C value; C > 0: if IERR = j, then LAPACK subroutine DHSEQR failed C to converge while computing the j-th eigenvalue. C C METHOD C C DHAEVS forms the upper Hessenberg matrix A' in (**) and calls C LAPACK subroutines to calculate its eigenvalues. The eigenvalues C of H are the square roots of the eigenvalues of A'. C C REFERENCES C C [1] Van Loan, C. F. C A symplectic method for approximating all the eigenvalues of C a Hamiltonian matrix. C Linear Algebra and its Applications 61 (1984), pp. 233-251. C C [2] Byers, R. C Hamiltonian and symplectic algorithms for the algebraic C Riccati equation. C Ph. D. Thesis, Cornell University, Ithaca, NY, January 1983. C C [3] Benner, P., Byers, R., and Barth, E. C Fortran 77 Subroutines for Computing the Eigenvalues of C Hamiltonian Matrices {I}: The Square-Reduced Method. C Submitted to ACM Trans. Math. Software, 1998. C C NUMERICAL ASPECTS C C The algorithm requires (32/3)*N**3 + O(N**2) floating point C operations. C Eigenvalues computed by this subroutine are exact eigenvalues C of a perturbed Hamiltonian matrix H + E where C C || E || <= c sqrt(eps) || H ||, C C c is a modest constant depending on the dimension N and eps is the C machine precision. Moreover, if the norm of H and an eigenvalue C are of roughly the same magnitude, the computed eigenvalue is C essentially as accurate as the computed eigenvalue obtained by C traditional methods. See [1] or [2]. C C CONTRIBUTOR C C P. Benner, Universitaet Bremen, Germany, and C R. Byers, University of Kansas, Lawrence, USA. C C REVISIONS C C 1994, February 1. C 1994, October 4. C 1995, March 22. C 1996, April 10. C 1997, December 15. C 1998, August 24. C C KEYWORDS C C Eigenvalues, (square-reduced) Hamiltonian matrix. C C ****************************************************************** * DOUBLE PRECISION A(LDA,N),QG(LDQG,N+1),RWORK(N* (N+1)),WI(N),WR(N) C .. C .. Parameters .. DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C .. C .. Local Scalars .. DOUBLE PRECISION SWAP,X,Y INTEGER I,IGNORE,IHI,ILO,J,M LOGICAL SORTED C .. C .. External Functions .. * . BLAS, LAPACK . LOGICAL LSAME EXTERNAL LSAME C .. C .. External Subroutines .. * . BLAS, LAPACK . * . others . EXTERNAL DCOPY,DCROOT,DGEBAL,DGEMM,DHSEQR,DLACPY,DSYMV C .. C .. Local Arrays .. DOUBLE PRECISION DUMMY(1) C .. * IERR = 0 IF (N.NE.0) THEN IF (.NOT. (LSAME(JOBSCL,'S').OR.LSAME(JOBSCL, + 'P').OR.LSAME(JOBSCL,'B').OR.LSAME(JOBSCL,'N'))) IERR = -1 IF (N.LT.0) IERR = -2 IF (LDA.LT.N) IERR = -4 IF (LDQG.LT.N) IERR = -6 * IF (IERR.EQ.0) THEN * .. 2 * .. form the eigenvalues of A' = A + GQ .. CALL DLACPY('L',N,N,QG,LDQG,RWORK,N) DO 10 I = 1,N CALL DCOPY(N-I+1,RWORK(I+N* (I-1)),1, + RWORK(I+N* (I-1)),N) 10 CONTINUE DO 20 J = 1,N CALL DSYMV('U',N,ONE,QG(1,2),LDQG,RWORK(1+N* (J-1)),1, + ZERO,WR,1) CALL DCOPY(N,WR,1,RWORK(1+N* (J-1)),1) 20 CONTINUE CALL DGEMM('N','N',N,N,N,ONE,A,LDA,A,LDA,ONE,RWORK,N) * .. 2 * .. Find the eigenvalues of A + GQ .. CALL DGEBAL(JOBSCL,N,RWORK,N,ILO,IHI,RWORK(1+N*N),IGNORE) CALL DHSEQR('E','N',N,ILO,IHI,RWORK,N,WR,WI,DUMMY,1, + RWORK(1+N*N),N,IERR) IF (IERR.EQ.0) THEN * .. * .. eigenvalues of H are the square roots .. DO 30 I = 1,N X = WR(I) Y = WI(I) CALL DCROOT(X,Y,WR(I),WI(I)) 30 CONTINUE * .. * .. Sort eigenvalues into decreasing order by real part * and, for eigenvalues with zero real part only, * decreasing order of imaginary part. (This simple * bubble sort preserves the relative order of * eigenvalues with equal but nonzero real part. * This ensures that complex conjugate pairs remain * together. .. SORTED = .FALSE. DO 50 M = N,1,-1 IF (SORTED) GO TO 60 SORTED = .TRUE. DO 40 I = 1,M - 1 IF (((WR(I).LT.WR(I+1)).OR. + ((WR(I).EQ.ZERO).AND. + (WR(I+1).EQ.ZERO).AND. + (WI(I).LT.WI(I+1))))) THEN SWAP = WR(I) WR(I) = WR(I+1) WR(I+1) = SWAP SWAP = WI(I) WI(I) = WI(I+1) WI(I+1) = SWAP * SORTED = .FALSE. * END IF * 40 CONTINUE 50 CONTINUE 60 CONTINUE * END IF * END IF * END IF * RETURN * *** Last Line of DHAEVS *** END SUBROUTINE DHASRD(COMPU,N,A,LDA,QG,LDQG,U,LDU,RWORK,IERR) C C .. Scalar Arguments .. INTEGER IERR,LDA,LDQG,LDU,N CHARACTER COMPU C .. C .. Array Arguments .. C C C PURPOSE C C To transform a Hamiltonian matrix C C ( A G ) C (*) H = ( T ) C ( Q -A ) C C into a square-reduced Hamiltonian matrix C C (A' G' ) C (**) H' = ( T) C (Q' -A' ) C T C by an orthogonal symplectic similarity transformation H' = U H U C where C ( U1 U2 ) C (***) U = ( ). C ( -U2 U1 ) C T C The square-reduced Hamiltonian matrix satisfies Q'A' - A' Q' = 0, C and C C 2 T 2 ( A'' G'') C H' := (U H U) = ( T ). C ( 0 A'') C C In addition, A'' is upper Hessenberg and G'' is skew symmetric. C The square roots of the eigenvalues of A'' = A'*A'+G'*Q' are the C eigenvalues of H. C C ARGUMENTS C C Mode Parameters C C COMPU CHARACTER*1 C Indicates whether the orthogonal symplectic similarity C transformation matrix U in (***) is returned or C accumulated into an orthogonal symplectic matrix or if the C transformation matrix is not required. C = 'N' or 'n': U is not referenced; C = 'V', 'v', 'A', or 'a': the orthogonal symplectic C similarity transformations are C accumulated into U. On input, U C must contain an orthogonal C symplectic matrix S; on exit , U C contains S*U with U from (**); C = 'I', 'i', 'F', or 'f': on entry, U need not be set, and C on exit, U contains the C orthogonal symplectic matrix U C from (**). C See the description of U below for details. C C N (input) INTEGER C The order of the matrices A, G, and Q. Normally, N.GE.1. C If N = 0, then IERR is set to 0, DHASRD returns C immediately without referencing any other argument. C C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) C On input, the leading N-by-N part of this array contains C the upper left block A of the Hamiltonian matrix H in (*). C On output, the leading N-by-N part of this array is C overwritten by the upper right block A' of the C square-reduced Hamiltonian matrix H' in (**). C C LDA (input) INTEGER C The leading dimension of array A as declared in the C calling program. LDA .GE. N. C C QG (input/output) DOUBLE PRECISION array, C dimension (LDQG,N+1) C On input QG contains the upper right symmetric block G and C the lower left symmetric block Q of the Hamiltonian matrix C in (*). On output these are overwritten by the C corresponding blocks of the square reduced Hamiltonian C matrix H' in (**). C The diagonal and lower triangle of the submatrix in C columns 1 to N contain the lower triangular part of Q. C The diagonal and upper triangle of the submatrix in C columns 2 to N+1 contains the upper triangle of G. C So, if I. GE. J, then Q(i,j) = Q(j,i) is stored in QG(i,j) C and G(i,j) = G(j,i) is in QG(j,i+1). C C LDQG (input) INTEGER C The leading dimension of array QG just as declared in the C calling procedure. LDQG. GE. N. C C U (input/output) DOUBLE PRECISION array, dimension (LDU,2*N) C If COMPU = 'N' or 'n', then U is not referenced. C If COMPU = 'I', 'i', 'F' or 'f', then the input contents C of U are not specified. On output the leading (N-by-2*N) C segment of array U is overwritten by the the first N rows C of the orthogonal symplectic matrix U in (**). C If COMPU = 'V', 'v', 'A' or 'a', then, on input, the C leading N-by-(2*N) segment of this array is assumed to C contain the first N rows of an orthogonal symplectic C matrix S. On output this array is overwritten by the first C N rows of the product S*U where U is the orthogonal C symplectic matrix from (**). C The storage scheme implied by (***) is used for orthogonal C symplectic matrices, i.e., only the first N rows are C stored as they contain all relevant information. C C LDU (input) INTEGER C The leading dimension of array U as declared in the C calling program. If COMPU is not 'N' or 'n', then C LDU .GE. N. C If COMPU = 'N' or 'n', then the array U is not referenced C and LDU may be set to 1. C C Workspace C C RWORK DOUBLE PRECISION array, dimension at least 2*N C C Error Indicator C C IERR INTEGER C = 0: successful exit; C < 0: if IERR = -j, then the j-th argument had an illegal C value. C C METHOD C C DHASRD transforms a Hamiltonian matrix H into a square-reduced C Hamiltonian matrix H' using the implicit version of Van Loan's C method as proposed in [1,2,3]. C C REFERENCES C C [1] Van Loan, C. F. C A symplectic method for approximating all the eigenvalues of C a Hamiltonian matrix. C Linear Algebra and its Applications 61 (1984), pp. 233-251. C C [2] Byers, R. C Hamiltonian and symplectic algorithms for the algebraic C Riccati equation. C Ph. D. Thesis, Cornell University, Ithaca, NY, January 1983. C C [3] Benner, P., Byers, R., and Barth, E. C Fortran 77 Subroutines for Computing the Eigenvalues of C Hamiltonian Matrices {I}: The Square-Reduced Method. C Submitted to ACM Trans. Math. Software, 1998. C C NUMERICAL ASPECTS C C This algorithm requires approximately 20*N**3 flops for C transforming H into square-reduced form. If the transformations C are required, this adds another 8*N**3 flops. The method is C strongly backward stable in the sense that if H' and U are the C computed square-reduced Hamiltonian and computed orthogonal C symplectic similarity transformation, then there is an orthogonal C symplectic matrix T and a Hamiltonian matrix M such that C C H T = T M C C || T - U || <= c1 * eps C C || H' - M || <= c2 * eps * || H || C C where c1, c2 are modest constants depending on the dimension N and C eps is the machine precision. C C Eigenvalues computed by explicitly forming the upper Hessenberg C matrix W = A'A' + G'Q' with A', G', Q' as in (**) and applying C the Hessenberg QR iteration to W are exactly eigenvalues of a C perturbed Hamiltonian matrix H + E, where C C || E || <= c3 * sqrt(eps) * || H ||, C C and c3 is a modest constant depending on the dimension N and eps C is the machine precision. Moreover, if the norm of H and an C eigenvalue lambda are of roughly the same magnitude, the computed C eigenvalue is essentially as accurate as the computed eigenvalue C from traditional methods. See [1] or [2]. C C CONTRIBUTOR C C P. Benner, Universitaet Bremen, Germany, and C R. Byers, University of Kansas, Lawrence, USA. C C REVISIONS C C 1982, July. C 1992, June. C 1995, March. C 1996, April. C 1997, December. C 1998, April. C 1998, August 26. C C KEYWORDS C C (Square-reduced) Hamiltonian matrix, symplectic orthogonal C similarity transformation. C C ****************************************************************** C DOUBLE PRECISION A(LDA,N),QG(LDQG,N+1),RWORK(2*N),U(LDU,*) C .. C .. Parameters .. C DOUBLE PRECISION ZERO,ONE,TWO PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0) C .. C .. Local Scalars .. DOUBLE PRECISION COSINE,SINE,TAU,TEMP,X,Y INTEGER I,J LOGICAL ACCUM,FORGET,FORM C .. C .. Local Arrays .. DOUBLE PRECISION T(2,2) C .. C .. External Functions .. DOUBLE PRECISION DDOT LOGICAL LSAME EXTERNAL DDOT,LSAME C .. C .. External Subroutines .. EXTERNAL DAXPY,DCOPY,DGEMV,DLARFG,DLARFX,DLARTG,DROT,DSYMV,DSYR2 C .. C IERR = 0 IF (N.NE.0) THEN ACCUM = LSAME(COMPU,'A') .OR. LSAME(COMPU,'V') FORM = LSAME(COMPU,'F') .OR. LSAME(COMPU,'I') FORGET = LSAME(COMPU,'N') IF (.NOT.ACCUM .AND. .NOT.FORM .AND. .NOT.FORGET) IERR = -1 IF (N.LT.0) IERR = -2 IF (LDA.LT.N) IERR = -4 IF (LDQG.LT.N) IERR = -6 IF (LDU.LT.N .AND. .NOT.FORGET) IERR = -8 IF (IERR.EQ.0) THEN C C .. TRANSFORM TO SQUARE REDUCED FORM .. DO 10 J = 1,N - 1 C T C .. RWORK <- (Q A - A Q)(J+1:N,J) .. CALL DCOPY(J-1,QG(J,1),LDQG,RWORK(N+1),1) CALL DCOPY(N-J+1,QG(J,J),1,RWORK(N+J),1) CALL DGEMV('T',N,N-J,-ONE,A(1,J+1),LDA,RWORK(N+1),1, + ZERO,RWORK(J+1),1) CALL DGEMV('N',N-J,J,ONE,QG(J+1,1),LDQG,A(1,J),1,ONE, + RWORK(J+1),1) CALL DSYMV('L',N-J,ONE,QG(J+1,J+1),LDQG,A(J+1,J),1, + ONE,RWORK(J+1),1) C C .. SYMPLECTIC REFLECTION TO ZERO (H*H)((N+J+2):2N,J) CALL DLARFG(N-J,RWORK(J+1),RWORK(J+1+1),1,TAU) Y = RWORK(J+1) RWORK(J+1) = ONE C CALL DLARFX('L',N-J,N,RWORK(J+1),TAU,A(J+1,1),LDA, + RWORK(N+1)) CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,A(1,J+1),LDA, + RWORK(N+1)) C CALL DLARFX('L',N-J,J,RWORK(J+1),TAU,QG(J+1,1),LDQG, + RWORK(N+1)) CALL DSYMV('L',N-J,TAU,QG(J+1,J+1),LDQG,RWORK(J+1),1, + ZERO,RWORK(N+J+1),1) CALL DAXPY(N-J,-TAU*DDOT(N-J,RWORK(N+J+1),1, + RWORK(J+1),1)/TWO,RWORK(J+1),1, + RWORK(N+J+1),1) CALL DSYR2('L',N-J,-ONE,RWORK(J+1),1,RWORK(N+J+1),1, + QG(J+1,J+1),LDQG) C CALL DLARFX('R',J,N-J,RWORK(J+1),TAU,QG(1,J+2),LDQG, + RWORK(N+1)) CALL DSYMV('U',N-J,TAU,QG(J+1,J+2),LDQG,RWORK(J+1),1, + ZERO,RWORK(N+J+1),1) CALL DAXPY(N-J,-TAU*DDOT(N-J,RWORK(N+J+1),1, + RWORK(J+1),1)/TWO,RWORK(J+1),1, + RWORK(N+J+1),1) CALL DSYR2('U',N-J,-ONE,RWORK(J+1),1,RWORK(N+J+1),1, + QG(J+1,J+2),LDQG) C IF (FORM) THEN C .. SAVE REFLECTION .. CALL DCOPY(N-J,RWORK(J+1),1,U(J+1,J),1) U(J+1,J) = TAU C ELSE IF (ACCUM) THEN C .. ACCUMULATE REFLECTION .. CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,U(1,J+1),LDU, + RWORK(N+1)) CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,U(1,N+J+1), + LDU,RWORK(N+1)) END IF C C .. (X,Y) := ((J+1,J),(N+J+1,J)) COMPONENT OF H*H .. X = DDOT(J,QG(1,J+2),1,QG(J,1),LDQG) + + DDOT(N-J,QG(J+1,J+2),LDQG,QG(J+1,J),1) + + DDOT(N,A(J+1,1),LDA,A(1,J),1) C C .. SYMPLECTIC ROTATION TO ZERO (H*H)(N+J+1,J) .. CALL DLARTG(X,Y,COSINE,SINE,TEMP) C CALL DROT(J,A(J+1,1),LDA,QG(J+1,1),LDQG,COSINE,SINE) CALL DROT(J,A(1,J+1),1,QG(1,J+2),1,COSINE,SINE) IF (N-J-1.GT.0) THEN CALL DROT(N-J-1,A(J+1,J+2),LDA,QG(J+2,J+1),1, + COSINE,SINE) CALL DROT(N-J-1,A(J+2,J+1),1,QG(J+1,J+3),LDQG, + COSINE,SINE) END IF C T(1,1) = A(J+1,J+1) T(1,2) = QG(J+1,J+2) T(2,1) = QG(J+1,J+1) T(2,2) = -T(1,1) CALL DROT(2,T(1,1),1,T(1,2),1,COSINE,SINE) CALL DROT(2,T(1,1),2,T(2,1),2,COSINE,SINE) A(J+1,J+1) = T(1,1) QG(J+1,J+2) = T(1,2) QG(J+1,J+1) = T(2,1) C IF (FORM) THEN C .. SAVE ROTATION .. U(J,J) = COSINE U(J,N+J) = SINE C ELSE IF (ACCUM) THEN C .. ACCUMULATE ROTATION .. CALL DROT(N,U(1,J+1),1,U(1,N+J+1),1,COSINE,SINE) END IF C C .. RWORK := (A + GQ)(J+1:N,J) .. CALL DGEMV('N',N-J,N,ONE,A(J+1,1),LDA,A(1,J),1,ZERO, + RWORK(J+1),1) CALL DCOPY(J-1,QG(J,1),LDQG,RWORK(N+1),1) CALL DCOPY(N-J+1,QG(J,J),1,RWORK(N+J),1) CALL DGEMV('T',J,N-J,ONE,QG(1,J+2),LDQG,RWORK(N+1),1, + ONE,RWORK(J+1),1) CALL DSYMV('U',N-J,ONE,QG(J+1,J+2),LDQG,RWORK(N+J+1), + 1,ONE,RWORK(J+1),1) C C .. SYMPLECTIC REFLECTION TO ZERO (H*H)(J+2:N,J) CALL DLARFG(N-J,RWORK(J+1),RWORK(J+1+1),1,TAU) RWORK(J+1) = ONE CALL DLARFX('L',N-J,N,RWORK(J+1),TAU,A(J+1,1),LDA, + RWORK(N+1)) CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,A(1,J+1),LDA, + RWORK(N+1)) C CALL DLARFX('L',N-J,J,RWORK(J+1),TAU,QG(J+1,1),LDQG, + RWORK(N+1)) CALL DSYMV('L',N-J,TAU,QG(J+1,J+1),LDQG,RWORK(J+1),1, + ZERO,RWORK(N+J+1),1) CALL DAXPY(N-J,-TAU*DDOT(N-J,RWORK(N+J+1),1, + RWORK(J+1),1)/TWO,RWORK(J+1),1, + RWORK(N+J+1),1) CALL DSYR2('L',N-J,-ONE,RWORK(J+1),1,RWORK(N+J+1),1, + QG(J+1,J+1),LDQG) C CALL DLARFX('R',J,N-J,RWORK(J+1),TAU,QG(1,J+2),LDQG, + RWORK(N+1)) CALL DSYMV('U',N-J,TAU,QG(J+1,J+2),LDQG,RWORK(J+1),1, + ZERO,RWORK(N+J+1),1) CALL DAXPY(N-J,-TAU*DDOT(N-J,RWORK(N+J+1),1, + RWORK(J+1),1)/TWO,RWORK(J+1),1, + RWORK(N+J+1),1) CALL DSYR2('U',N-J,-ONE,RWORK(J+1),1,RWORK(N+J+1),1, + QG(J+1,J+2),LDQG) C IF (FORM) THEN C .. SAVE REFLECTION .. CALL DCOPY(N-J,RWORK(J+1),1,U(J+1,N+J),1) U(J+1,N+J) = TAU C ELSE IF (ACCUM) THEN C .. ACCUMULATE REFLECTION .. CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,U(1,J+1),LDU, + RWORK(N+1)) CALL DLARFX('R',N,N-J,RWORK(J+1),TAU,U(1,N+J+1), + LDU,RWORK(N+1)) END IF C 10 CONTINUE C C IF (FORM) THEN C .. FORM S BY ACCUMULATING TRANSFORMATIONS .. DO 40 J = N - 1,1,-1 C .. INITIALIZE (J+1)ST COLUMN OF S .. DO 20 I = 1,N U(I,J+1) = ZERO 20 CONTINUE U(J+1,J+1) = ONE DO 30 I = 1,N U(I,N+J+1) = ZERO 30 CONTINUE C C .. SECOND REFLECTION .. TAU = U(J+1,N+J) U(J+1,N+J) = ONE CALL DLARFX('L',N-J,N-J,U(J+1,N+J),TAU,U(J+1,J+1), + LDU,RWORK(N+1)) CALL DLARFX('L',N-J,N-J,U(J+1,N+J),TAU, + U(J+1,N+J+1),LDU,RWORK(N+1)) C .. ROTATION .. CALL DROT(N-J,U(J+1,J+1),LDU,U(J+1,N+J+1),LDU, + U(J,J),U(J,N+J)) C .. FIRST REFLECTION .. TAU = U(J+1,J) U(J+1,J) = ONE CALL DLARFX('L',N-J,N-J,U(J+1,J),TAU,U(J+1,J+1), + LDU,RWORK(N+1)) CALL DLARFX('L',N-J,N-J,U(J+1,J),TAU,U(J+1,N+J+1), + LDU,RWORK(N+1)) 40 CONTINUE C .. FIRST COLUMN IS FIRST COLUMN OF IDENTITY .. DO 50 I = 1,N U(I,1) = ZERO U(I,N+1) = ZERO 50 CONTINUE U(1,1) = ONE END IF C END IF C END IF C .. LAST LINE OF DHASRD .. END SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0