C ALGORITHM 718, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 2, JUNE, 1993, PP. 224-232. CXXXXXXXXXXXXXXXXXXXXX START OF FILE: README XXXXXXXXXXXXXXXXXXXXX ************************ DSEVAS/SSEVAS ************************ This diskette contains the source code for SEVAS, a FORTRAN subroutine that solves the single-input eigenvalue allocation problem. Files beginning with 'D' are for the double precision version, while those beginning with 'S' are for single precision. Files used: DSEVAS/SSEVAS -eigenvalue allocation subroutine DTEST/STEST -test driver program DFLRD/SFLRD -test driver that reads user input from a file DEIGVL/SEIGVL -subroutine to find the eigenvalues of a general matrix DBLAS/SBLAS -Level 1 BLAS routines used by DSEVAS/SSEVAS Makefile -constructs DTEST/STEST/DFLRD/SFLRD DDAT01.DAT -sample data file for DFLRD/SFLRD README -this file BUILDING PROGRAMS: ------------------ (The discussion here works for the single-precision version as well. Replace DTEST by STEST, and so on.) There are two top-level programs: DTEST and DFLRD. DTEST is the main test program that generates its own data, while DFLRD takes user-generated data from a file. To build: 1) compile DTEST (and/or DFLRD), DSEVAS, DBLAS, DEIGVL. 2) link the above routines, using either DTEST or DFLRD as the main (executable). On most machines, listing DTEST(DFLRD) first on the linker command line is necessary. EXAMPLES: --------- 1) For UNIX systems: ==================== Type "make dtest" or "make dflrd" in the directory containing the SEVAS subroutines. To compile DTEST explicitly, use the following: f77 -c dtest.f f77 -c dsevas.f f77 -c deigvl.f f77 -c dblas.f f77 -o dtest dsevas.o deigvl.o dblas.o (You may need to use "fc" or some other compiler name instead of "f77".) To run, type: dtest 2) For VMS systems: =================== Type the following: FORTRAN DTEST.F FORTRAN DSEVAS.F FORTRAN DEIGVL.F FORTRAN DBLAS.F LINK DTEST,DSEVAS,DEIGVL,DBLAS To run, type: RUN DTEST #XXXXXXXXXXXXXXXXXXXXX START OF FILE: MAKEFILE XXXXXXXXXXXXXXXXXXXXX #=========================== # Makefile for dsevas/ssevas #=========================== # # The compiler name may be different on # your machine. # COMP=f77 # The following files are used: # # dtest/stest -test program # dflrd/sflrd -user input test program # dsevas/ssevas -eigenvalue allocation routine # deigvl/eigvl -eigenvalue problem routine (for test drivers) # dblas/sblas -BLAS1 routines used by DSEVAS/SSEVAS # DOBJ = dtest.o dsevas.o deigvl.o dblas.o dtest: ${DOBJ} ${COMP} -o $@ ${DOBJ} DRDOBJ = dflrd.o dsevas.o deigvl.o dblas.o dflrd: ${DRDOBJ} ${COMP} -o $@ ${DRDOBJ} SOBJ = stest.o ssevas.o seigvl.o sblas.o stest: ${SOBJ} ${COMP} -o $@ ${SOBJ} SRDOBJ = sflrd.o ssevas.o seigvl.o sblas.o sflrd: ${SRDOBJ} ${COMP} -o $@ ${SRDOBJ} clean: rm *.o 5 T 0.0 1.000000000000000 0.5000000000000000 0.3333333333333333 0.2500000000000000 0.2000000000000000 2.000000000000000 1.000000000000000 0.6666666666666666 0.5000000000000000 0.4000000000000000 0.0000000000000000E+00 1.500000000000000 1.000000000000000 0.7500000000000000 0.6000000000000000 0.0000000000000000E+00 0.0000000000000000E+00 1.333333333333333 1.000000000000000 0.8000000000000000 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 1.250000000000000 1.000000000000000 1.000000000000000 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 1.000000000000000 0.0000000000000000E+00 2.000000000000000 0.0000000000000000E+00 3.000000000000000 0.0000000000000000E+00 4.000000000000000 0.0000000000000000E+00 5.000000000000000 0.0000000000000000E+00 CXXXXXXXXXXXXXXXXXXXXX START OF FILE: SSEVAS.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C********************************************************************* C C Title: SSEVAS C C Authors: G. Miminis and M. Reid C C Purpose: Calculates f such that (A - b*f) has C eigenvalues eigvl. A and b may be C general, or in "controllability" form, C where A is unreduced upper hessenberg C and b has one non-zero element in b(1). C The eigenvalues stored in eigvl may be C real or complex conjugate pairs. C C C Algorithm: C ---------- C C This subroutine implements an algorithm by G. Miminis and C C. Paige, submitted to the SIAM Journal on Matrix Analysis C and Applications. For details please see the paper. C C C Input Arguments: C ---------------- C C [ (+) means argument is altered by the subroutine. ] C C (+) A real(lda,n) Matrix A C C n integer Order of matrix A C C lda integer First dimension of matrix A as C defined in the calling program C C (+) binp real(n) Input vector b. C C eigvl real(2,n) Desired eigenvalues. EIGVL(1,*) C holds the real part; EIGVL(2,*) C holds the imaginary. C C (+) rstor real(4,*) Working storage for reflectors. C Must be dimensioned as C RSTOR(4,k), where: C C k >= (n**2 - 2*n + 1)/4 C C C (+) irpos integer(*) Position of reflectors. Must C be dimensioned as IRPOS(k), C where k is defined as for C RSTOR. C C (+) cstor real(*) Working storage for C conversion to controllability C form. Must be dimensioned as C CSTOR(w), where: C C w >= (n**2 + 3*n - 4)/2 C C C iform logical Input form. C .TRUE. = controllability C form C .FALSE. = general C C C (+) rmu real*4 The machine unit. C If positive, the subroutine C will use the given value C as the machine unit. C If the given value of rmu is C zero or negative the subroutine C will automatically calculate the C machine unit. (See note 3 C below.) The calculated C machine unit is returned via C this parameter. C C The tolerance used in testing C for controllability is: C C tol = ||(b,A)|| * u , C C where u is the machine unit C and ||.|| is the 1-norm. C C Output Arguments: C -------------- C C (+) f real(n) Solution vector C C (+) ieigal integer Number of eigenvalues actually C allocated. C C Local Variables: C ---------------- C C isize integer Current size of matrix A. C C le integer Position of leading edge of C matrix A in array A. C C irefc integer Index of last elementary C reflector stored in rstor. C C eps real Tolerance used in testing C for controllability. C C C C Subroutines Used: C ----------------- C C SAPPSH, SCRH, SPREM, SGTEPS, SEIGCK, SCNTMK, SPOST, STORE C C C Notes on Usage: C --------------- C C 1) SSEVAS only takes eigenvalues that are real or in C complex conjugate pairs. These are passed to the C subroutine via EIGVL, dimensioned as EIGVL(2,N). C EIGVL(1,x) holds the real part, while EIGVL(2,x) holds C the imaginary part. Real valued eigenvalues can be in C any order. Conjugate pairs of eigenvalues must be C stored contiguously, with the first of the pair on an C odd-numbered index in the array EIGVL. For example, C the eigenvalues 1,2,3, and 1 +/- 4 can be stored: C C [1.0,0.0] C [2.0,0.0] C [1.0,+4.0] C [1.0,-4.0] C [3.0,0.0] C C but the following would generate an error: C C [1.0,0.0] C [1.0,+4.0] C [1.0,-4.0] C [2.0,0.0] C [3.0,0.0] C C since the subroutine working with two eigenvalues C at a time, will work with eigenvalues 1.0 and C 1.0 + i*4.0 , which is not allowed. An easy way C to list the eigenvalues could be to first list C all the complex eigenvalues, followed by the real C eigenvalues. C C 2) Before allocating a pair of eigenvalues, the subroutine C checks to see if the current matrix is controllable. C If this is not the case, no further allocations will C take place, and the subroutine will print a warning C indicating the number of successful allocations. C The resulting matrix (A-b*f) will have those C eigenvalues allocated; the remainder will be undefined. C C 3) The SSEVAS routine uses the machine unit as part of the C tolerance for controllability determination. If C the input parameter RMU is zero or negative, the C internal routine SGTEPS automatically calculates C the machine unit. The user should note, however, C that such a calculation can fail and give the C double precision machine unit, even in the single C precision version. (See Gentleman and Marovich C [CACM 17, p. 276, 1974]). For this reason, the user C has the option of supplying the machine unit, C possibly from the PORTS machine constants sub- C routines. (See TOMS 4 (1978), p. 177). The C parameter RMU returns the machine unit actually C used in the calculations. C C 4) This subroutine makes extensive use of the level-1 BLAS C (Basic Linear Algebra Subroutines) package. The C BLAS-1 subroutines used are: C C SASUM, SAXPY, SCOPY, SDOT, SNRM2, SSCAL, ISAMAX C C These routines are available from various sources, C such as the IMSL subroutine library. C C 5) A named common block (MACEPS) is used in this subroutine. C C C********************************************************************* C********************************************************************* C********************************************************************* SUBROUTINE SSEVAS(A,N,LDA,BINP,EIGVL,RSTOR,IRPOS, + CSTOR,IFORM,RMU,IEIGAL,F) C C........input arguments C INTEGER N,LDA,IRPOS(*),IEIGAL REAL A(LDA,*),BINP(N),F(N),RSTOR(4,*),CSTOR(*) REAL EIGVL(2,N),RMU LOGICAL IFORM C C.......local variables C INTEGER ISIZE,LE,IREFC,IOLDRF REAL S1,S2,LM,EPS,SGTEPS LOGICAL SEIGCK C C.......convert general form to controllability form C IF(.NOT. IFORM) CALL SCNTMK(A,N,LDA,BINP,CSTOR) C C.......calculate epsilon for controllability determination C EPS = SGTEPS(A,N,LDA,BINP,RMU) C C.......clear zero elements of system C DO 50 J=1,N-2 F(J) = 0.0 DO 50 I=J+2,N A(I,J) = 0.0 50 CONTINUE F(N-1) = 0.0 F(N) = 0.0 C C.......start deflation: continue by double-step until problem is C smaller than size 3 C IEIGAL = 0 IREFC = 0 ISIZE = N 100 IF ( ISIZE .LE. 2) GOTO 200 LE = N - ISIZE + 1 IOLDRF = IREFC C C.......check eigenvalues for complex conjugates C IF(.NOT. SEIGCK(EIGVL,LE,EPS)) THEN WRITE(*,5010) WRITE(*,5005) N-ISIZE GOTO 900 ENDIF C C.......apply reflector containing eigenvalues to be allocated. C CALL SAPPSH(A,LE,N,LDA,EIGVL,RSTOR,IRPOS,IREFC) C C.......transform A to block upper hessenberg. C CALL SCRH(A,LE,N,LDA,RSTOR,IRPOS,IREFC) C C.......apply latest reflector to binp (vector of inputs) C BINP(LE+1) = 0.0 BINP(LE+2) = 0.0 CALL SPREM(BINP,LDA,LE,1,1,RSTOR(1,IREFC),3) C C......check system for controllability C IF(ABS(BINP(LE+2)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS WRITE(*,5005) LE-1 IREFC = IOLDRF F(LE) = 0.0 F(LE+1) = 0.0 GOTO 900 C C.......specify le,le+1 elements of f C ELSE F(LE) = A(LE+2,LE) / BINP(LE+2) F(LE+1) = A(LE+2,LE+1) / BINP(LE+2) ENDIF C C.......reduce problem size for next step C ISIZE = ISIZE - 2 IEIGAL = IEIGAL + 2 GOTO 100 200 CONTINUE C C.......allocate eigenvalues for matrices of size one or two. C IF (ISIZE .EQ. 1) THEN C C...........check for controllability C IF(ABS(BINP(N)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS WRITE(*,5005) N-1 F(N) = 0.0 GOTO 900 ELSE F(N) = (A(N,N) - EIGVL(1,N))/BINP(N) IEIGAL = IEIGAL + 1 ENDIF ELSE C C...........check for controllability and conjugate C eigenvalues (size 2) C IF(.NOT. SEIGCK(EIGVL,LE,EPS) .OR. + ABS(BINP(N-1)) .LT. EPS .OR. + ABS(A(N,N-1)) .LT. EPS) THEN IF(ABS(BINP(N-1)) .LT. EPS .OR. + ABS(A(N,N-1)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS ELSE WRITE(*,5010) ENDIF WRITE(*,5005) N-2 F(N-1) = 0.0 F(N) = 0.0 GOTO 900 ELSE LM = EIGVL(1,N-1)*EIGVL(1,N) - EIGVL(2,N-1)*EIGVL(2,N) S1 = EIGVL(1,N-1)+EIGVL(1,N) - A(N,N) S2 = (S1*A(N,N) - LM)/A(N,N-1) F(N-1) = (A(N-1,N-1) - S1)/BINP(N-1) F(N) = (A(N-1,N) - S2)/BINP(N-1) IEIGAL = IEIGAL + 2 ENDIF ENDIF 900 CONTINUE C C.......apply accumulated reflectors to f vector C DO 300 K=IREFC,1,-1 CALL SPREM(F,LDA,IRPOS(K),1,1,RSTOR(1,K),3) 300 CONTINUE C C.......reverse transformation to controllability form for C vector f C IF(.NOT. IFORM) THEN II=N*(N+1)/2 + N - 4 DO 350 K=2,N CALL SPOST(F,1,N-K+1,N-K+1,1,CSTOR(II),K) II = II - K - 2 350 CONTINUE ENDIF 5000 FORMAT(1X,'SYSTEM IS CLOSE TO UNCONTROLLABLE') 5002 FORMAT(1X,'USING TOLERANCE ',E20.12) 5005 FORMAT(1X,I4,' EIGENVALUES ALLOCATED.') 5010 FORMAT(1X,'COMPLEX EIGENVALUES ARE NOT COMPLEX CONJUGATES.') RETURN END C********************************************************************* C C Title: SCNTMK C C Purpose: Converts the system (b,A) into C controllability form. C C C Input Arguments: C ---------------- C C A real(lda,n) Control matrix A C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C binp real(n) Input vector C C cstor real(*) Working storage for C conversion to controllability C form. C C C C Local Variables: C ---------------- C C ipnt integer Pointer to current reflector C in CSTOR() C C C Subroutines Used: C ----------------- C C SCOPY, SNRM2, SVHOUS, SPREM, SPOST C C********************************************************************* SUBROUTINE SCNTMK(A,N,LDA,BINP,CSTOR) INTEGER N,LDA,IPNT REAL A(LDA,*),BINP(N),CSTOR(*),SNRM2 IPNT = 1 DO 100 I=N,2,-1 IF(I .EQ. N) THEN CALL SCOPY(I, BINP(1), 1, CSTOR(IPNT+1), 1) BINP(1) = -SIGN(SNRM2(I, BINP(1), 1),BINP(1)) DO 150 J=2,N BINP(J) = 0.0 150 CONTINUE ELSE CALL SCOPY(I, A(N-I+1,N-I), 1, CSTOR(IPNT+1), 1) A(N-I+1,N-I) = -SIGN(SNRM2(I, A(N-I+1,N-I), 1), + A(N-I+1,N-I)) DO 175 J=N-I+2,N A(J,N-I) = 0.0 175 CONTINUE ENDIF CALL SVHOUS(CSTOR(IPNT), I) CALL SPREM(A,LDA,N-I+1,N-I+1,N,CSTOR(IPNT),I) CALL SPOST(A,LDA,1,N,N-I+1,CSTOR(IPNT),I) IPNT = IPNT + I + 1 100 CONTINUE RETURN END C********************************************************************* C C Title: SVHOUS C C Purpose: Given a n-element vector stored in v(2)- C v(n+1), SVHOUS calculates pi and u such that C I - u * uT / pi will introduce zeroes into C first two elements of v. Pi is returned in C v(1), while u overwrites the input vector C v(2)-v(n+1). C C Input Arguments: C ---------------- C C v real(*) The last n elements of v contain C the input vector to be zeroed. C C C Output Values: C -------------- C C v v(1) = pi C v(2)-v(n+1) = u C C C Subroutines used: C ----------------- C C SSCAL, SNRM2, ISAMAX C C********************************************************************* C C.......overwrites v with u of elem. reflector; returns p C SUBROUTINE SVHOUS(V,N) INTEGER N,I REAL V(*),SIGMA,ALFA,SNRM2 I = ISAMAX( N, V(2), 1) ALFA = 1.0/ABS(V(I+1)) CALL SSCAL( N, ALFA, V(2), 1) SIGMA = SIGN(SNRM2(N, V(2), 1),V(2)) V(2) = V(2) + SIGMA V(1) = SIGMA*V(2) RETURN END C********************************************************************* C C Title: SEIGCK C C Purpose: Checks that complex eigenvalues are C complex conjugates. The function C returns .TRUE. if the eigenvalues are C real or complex conjugates, .FALSE. C otherwise. C C C C Input Arguments: C ---------------- C C eigvl real(2,n) Array of eigenvalues to C be allocated. C C le integer Position of leading edge of C matrix A in array A. C C eps real Epsilon value for C comparison to zero. C C eps = ||(b,A)||u C C where ||.|| is the C 1-norm of the system, C and u is the machine unit. C C Local Variables: C ---------------- C C rd1 real Relative difference in real C parts. C C rd2 real Relative difference in complex C parts. C C norm,tnorm,eps C C C Subroutines Used: C ----------------- C C SASUM C C********************************************************************* LOGICAL FUNCTION SEIGCK(EIGVL,LE,EPS) REAL EIGVL(2,*),EPS,RD1,RD2 INTEGER LE SEIGCK = .TRUE. C C.......check if complex eigenvalues are complex conjugates C IF(ABS(EIGVL(2,LE)) .GT. EPS .OR. + ABS(EIGVL(2,LE+1)) .GT. EPS) THEN RD1 = ABS((EIGVL(1,LE) - EIGVL(1,LE+1))/EIGVL(1,LE+1)) RD2 = ABS((EIGVL(2,LE) + EIGVL(2,LE+1))/EIGVL(2,LE+1)) IF(RD1 .GT. EPS .OR. RD2 .GT. EPS) THEN SEIGCK = .FALSE. ENDIF ENDIF RETURN END C********************************************************************* C C Title: SGTEPS C C C Purpose: Calculates the tolerance used in C determining the uncontrollability of C a system. To compute the tolerance, the C routine computes the 1-norm of the system C (b,A) and multiplies this value by the C machine unit: C C eps = ||(b,A)||*u C C where ||.|| is the 1-norm, and u is the C machine unit. C C If the input parameter RMU is positive, C that value is used as the machine unit in C calculating EPS. If RMU is zero or negative, C the routine calculates the machine unit. C C C C Input Arguments: C ---------------- C C A real(lda,n) Control matrix A C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C b real(n) Input vector b C C rmu real User input value for machine C unit C C Local Variables: C ---------------- C C norm,tnorm,eps C C C Subroutines Used: C ----------------- C C SASUM,SSTORE C C********************************************************************* REAL FUNCTION SGTEPS(A,N,LDA,B,RMU) INTEGER N,LDA REAL A(LDA,*),B(N),RMU REAL NORM,TNORM,EEPS,EPS,SASUM COMMON/MACEPS/EEPS C C.......use user value for machine unit if non-zero C IF(RMU .GT. 0.0) THEN EPS = RMU ELSE C.......calculate machine unit EPS = 1.0 30 EPS = EPS / 2.0 EEPS=EPS + 1.0 CALL SSTORE(EEPS) IF(EEPS .NE. 1.0) GOTO 30 EPS = EPS * 2.0 RMU = EPS ENDIF C C.......calculate 1 norm of system (b,A) C NORM = SASUM(N, B(1), 1) DO 50 I=1,N TNORM = SASUM(N, A(1,I), 1) IF(TNORM .GT. NORM) NORM = TNORM 50 CONTINUE SGTEPS = EPS * NORM RETURN END C C.......This subroutine ensures that the value E is actually C stored in memory, with appropriate rounding errors. C (From Gentleman and Marovich, Communications of the ACM, C vol. 17, no. 5, May 1974. p. 277) C C SUBROUTINE SSTORE(E) REAL E,V COMMON/MACEPS/V V = E RETURN END C********************************************************************* C C Title: SCRH C C Purpose: Uses a series of elementary reflectors C to restore the matrix A to a "near" upper C Hessenberg form. C C C Input Arguments: C ---------------- C C A real(lda,n) Control matrix A C C le integer Position of leading edge of C matrix A in array A. C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C rstor real(4,*) Working storage for C reflectors C C irpos integer(*) Position of reflectors C C irefc integer Index of last elementary C reflector stored in rstor. C C C Output Values: C -------------- C C A A is returned, transformed into near C upper Hessenberg. C C rstor,irpos,irefc Contain reflectors used to transform C A to near upper Hessenberg. C C C C Subroutines Used: C ----------------- C C SCOPY, SHOUSE, SNRM2, SPOST, SPREM C C********************************************************************* SUBROUTINE SCRH(A,LE,N,LDA,RSTOR,IRPOS,IREFC) INTEGER LE,N,LDA,I,IREFC,IRPOS(*) REAL A(LDA,*),P,RSTOR(4,*),SNRM2 DO 100 I=N,LE+3,-1 IREFC = IREFC + 1 IRPOS(IREFC) = I-3 CALL SCOPY(3, A(I,I-3), LDA, RSTOR(2,IREFC), 1) CALL SHOUSE(RSTOR(1,IREFC)) P = -SIGN(SNRM2(3, A(I,I-3), LDA),A(I,I-1)) CALL SPOST(A,LDA,LE,I-1,I-3,RSTOR(1,IREFC),3) A(I,I-1) = P A(I,I-2) = 0.0 A(I,I-3) = 0.0 CALL SPREM(A,LDA,I-3,LE,N,RSTOR(1,IREFC),3) 100 CONTINUE RETURN END C********************************************************************* C C Title: SAPPSH C C Purpose: Calculates an elementary reflector C containing the implicit shift and pre- C and post multiplies A. C C Input Arguments: C ---------------- C C A real(lda,n) Control matrix A. C C le integer Position of leading edge of C matrix A in array A. C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C eigvl real(2,n) Vector of eigenvalues to be C allocated. C C rstor real(4,*) Working storage for reflectors. C C irpos integer(*) Position of reflectors. C C irefc integer Index of last reflector in C irpos and rstor. C C C Output Values: C -------------- C C A Returns the original matrix A pre- and C post-multiplied by the elementary reflector. C C rstor,irpos Stores the value and position of the C elementary reflector. C C la,lm C C Subroutines used: C ----------------- C C SHOUSE, SPREM, SPOST C C********************************************************************* SUBROUTINE SAPPSH(A,LE,N,LDA,EIGVL,RSTOR,IRPOS,IREFC) INTEGER LE,N,LDA,IRPOS(*),IREFC REAL A(LDA,*),RSTOR(4,*) REAL EIGVL(2,*) REAL LA,LM IREFC = IREFC + 1 IRPOS(IREFC) = N-2 C C.......construct bottom row of A~ C LA = EIGVL(1,LE) + EIGVL(1,LE+1) LM = EIGVL(1,LE)*EIGVL(1,LE+1) - EIGVL(2,LE)*EIGVL(2,LE+1) RSTOR(2,IREFC) = A(N,N-1)*A(N-1,N-2) RSTOR(3,IREFC) = A(N,N-1)*A(N-1,N-1) + A(N,N)*A(N,N-1) - + LA*A(N,N-1) RSTOR(4,IREFC) = A(N,N-1)*A(N-1,N) + A(N,N)*A(N,N) - + LA*A(N,N) + LM C C.......convert A to PtAP C CALL SHOUSE(RSTOR(1,IREFC)) CALL SPOST(A,LDA,LE,N,N-2,RSTOR(1,IREFC),3) CALL SPREM(A,LDA,N-2,LE,N,RSTOR(1,IREFC),3) RETURN END C********************************************************************* C C Title: SPREM C C Purpose: Premultiplies the matrix contained in array A C by the 3-element elementary reflector stored C in vector v. C C Input Arguments: C ---------------- C C A real(lda,nv) Control matrix A. C C lda integer First dimension of array A as C defined in the calling program C C irow integer First row of A affected by C the reflector. (Reflector will C be applied to rows *irow* to C *irow+nv-1*.) C C j,k integer Range of columns of A on C which to apply reflector. C C v real(4) Contains the elementary C reflector. v(1) = pi and C v(2)->v(nv+1) = u, where C I - u*uT/pi is the reflector. C C nv integer Dimension of vector v. C C C Output Values: C -------------- C C A Returns the original matrix A pre- C multiplied by the elementary reflector. C C C C Subroutines used: C ----------------- C C SDOT, SAXPY C C********************************************************************* SUBROUTINE SPREM(A,LDA,IROW,J,K,V,NV) INTEGER IROW,J,K,LDA,I,NV REAL A(LDA,*),V(*),ALPHA,SDOT DO 100 I=J,K ALPHA = SDOT(NV, A(IROW,I), 1, V(2), 1)/V(1) CALL SAXPY(NV, -ALPHA, V(2), 1, A(IROW,I), 1) 100 CONTINUE RETURN END C********************************************************************* C C Title: SPOST C C Purpose: post-multiplies the matrix contained in C array A by the 3-element elementary reflector C stored in vector v. C C Input Arguments: C ---------------- C C A real(lda,n) Control matrix A. C C lda integer First dimension of array A as C defined in the calling program C C j,k integer Range of rows of A on C which to apply reflector. C C icol integer First column of A affected by C the reflector. (Reflector will C be applied to columns *icol* C to *icol+nv-1*.) C C v real(4) Contains the elementary C reflector. v(1) = pi and C v(2)->v(nv+1) = u, where C I - u*uT/pi is the reflector. C C nv integer Dimension of vector v. C C C Output Values: C -------------- C C A Returns the original matrix A C post-multiplied by the elementary reflector. C C C C Subroutines used: C ----------------- C C SDOT, SAXPY C C********************************************************************* SUBROUTINE SPOST(A,LDA,J,K,ICOL,V,NV) INTEGER J,K,ICOL,LDA,I,NV REAL A(LDA,*),V(*),ALPHA,SDOT DO 100 I=J,K ALPHA = SDOT(NV, A(I,ICOL), LDA, V(2), 1)/V(1) CALL SAXPY(NV, -ALPHA, V(2), 1, A(I,ICOL), LDA) 100 CONTINUE RETURN END C********************************************************************* C C Title: SHOUSE C C Purpose: Given a 3-element vector stored in v(2)- C v(4), SHOUSE calculates pi and u such that C I - u * uT / pi will introduce zeroes into C first two elements of v. Pi is returned in C v(1), while u overwrites the input vector C v(2)-v(4). C C Input Arguments: C ---------------- C C v real(4) The last three elements of v contain C the input vector to be zeroed. C C C Output Values: C -------------- C C v v(1) = pi C v(2)-v(4) = u C C C Subroutines used: C ----------------- C C SSCAL, SNRM2, ISAMAX C C********************************************************************* C C.......overwrites v with u of elem. reflector; returns p C SUBROUTINE SHOUSE(V) INTEGER I REAL V(4),SIGMA,ALFA,SNRM2 I = ISAMAX( 3, V(2), 1) ALFA = 1.0/ABS(V(I+1)) CALL SSCAL( 3, ALFA, V(2), 1) SIGMA = SIGN(SNRM2(3, V(2), 1),V(4)) V(4) = V(4) + SIGMA V(1) = SIGMA*V(4) RETURN END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: DTEST.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C C TEST DRIVER: DSEVAS C C These test segments use DEIGVL, a public domain routine C for computing the eigenvalues of a matrix. The code C for this subroutine can be found in DEIGVL.F, and was C originally found in the EISPACK subroutine library. C (Any other similar subroutine can be substituted.) C The DSEVAS subroutine requires the level 1 BLAS routines, C which are included in the file DBLAS.F, and can also be C found in the IMSL library. C C********************************************************************* C********************************************************************* PROGRAM SAMPLE PARAMETER (ISL = 25, JSL = ISL*ISL+2*ISL) DOUBLE PRECISION A(ISL,ISL),AA(ISL,ISL),B(ISL),BB(ISL) DOUBLE PRECISION F(ISL),RST(4,JSL),CST(JSL) INTEGER IST(JSL),N,IEAL COMPLEX*16 EIG(ISL) DOUBLE PRECISION L(2,ISL),RMU LOGICAL FORM C********************************************************************* C C Example 1: system in controllability form. C C This program segment computes a system in controllability C form and computes the vector f such that the matrix (A - b*f) C has the eigenvalues (1,2,3,4,5). C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 5 DO 100 I=1,N L(1,I) = DBLE(I) L(2,I) = 0.0 B(I) = 0.0 BB(I) = 0.0 DO 100 J=I-1,N IF(J .GT. 0) THEN A(I,J) = DBLE(I)/DBLE(J) AA(I,J) = A(I,J) ENDIF 100 CONTINUE B(1) = 1.0 BB(1) = 1.0 WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 1' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 110 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 110 CONTINUE C C.......call DSEVAS routine C FORM = .TRUE. RMU = 0.0 CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 120 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 120 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 140 I=1,N WRITE(*,*) EIG(I) 140 CONTINUE WRITE(*,*) WRITE(*,*) ' (THE MATRIX (A - B*F) SHOULD HAVE EIGENVALUES' WRITE(*,*) ' 1.0, 2.0, 3.0, 4.0 and 5.0, BUT NOT NECESSARILY' WRITE(*,*) ' IN THIS ORDER. THIS EXAMPLE USES A SYSTEM IN' WRITE(*,*) ' CONTROLLABILITY FORM.)' WRITE(*,*) C********************************************************************* C C Example 2: near-uncontrollable system in controllability C form. C C This program segment computes a system in controllability C form and attempts to compute the vector f such that the C matrix (A - b*f) has the eigenvalues (1,2,3,4,5). C The subroutine will detect that the system is near to C an uncontrollable one, and stop allocations. C C********************************************************************* C C.......create system (in controllability form) C N = 5 DO 200 I=1,N L(1,I) = DBLE(I) L(2,I) = 0.0 B(I) = 0.0 BB(I) = 0.0 F(I) = 0.0 DO 200 J=1,N IF(J .GE. I-1) THEN A(I,J) = DBLE(I)/DBLE(J) AA(I,J) = A(I,J) ELSE A(I,J) = 0.0 AA(I,J) = 0.0 ENDIF 200 CONTINUE B(1) = 1.0 BB(1) = 1.0 C C.......make sub-diagonal element zero. C A(3,2) = 0.0 AA(3,2) = A(4,3) WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 2' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 210 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 210 CONTINUE C C.......call DSEVAS routine C FORM = .TRUE. RMU = 0.0 CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 220 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 220 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F):' WRITE(*,*) DO 240 I=1,N WRITE(*,*) EIG(I) 240 CONTINUE WRITE(*,*) WRITE(*,*) ' (SINCE THE SYSTEM (b,A) IS UNCONTROLLABLE,' WRITE(*,*) ' NO EIGENVALUES ARE ALLOCATED. THE VECTOR F' WRITE(*,*) ' SHOULD BE ZERO, AND THE EIGENVALUES OF' WRITE(*,*) ' (A - B*F) SHOULD BE THOSE OF THE MATRIX A)' WRITE(*,*) C********************************************************************* C C Example 3: near-uncontrollable system in general form. C C This program segment computes a system in general C form and attempts to compute the vector f such that the C matrix (A - b*f) has the eigenvalues (1,2,3,4,5). C The subroutine will detect that the system is near to C an uncontrollable one after allocating two eigenvalues, C and stop further allocations. The resulting matrix (A - b*f) C will thus have eigenvalues 1, 2, and three other undetermined C values. C C********************************************************************* C C.......create general system C N = 5 DO 300 I=1,N L(1,I) = DBLE(I) L(2,I) = 0.0 B(I) = 1.0/DBLE(I) BB(I) = B(I) F(I) = 0.0 DO 300 J=1,N A(I,J) = (DBLE(I) - 3.0) / DBLE(J) IF(I .EQ. J+1 .AND. J .LT. 3) + A(I,J) = A(I,J)*20.0 AA(I,J) = A(I,J) 300 CONTINUE WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 3 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 310 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 310 CONTINUE C C.......call DSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 320 I=1,N WRITE(*,*) F(I) DO 320 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 320 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F):' WRITE(*,*) DO 340 I=1,N WRITE(*,*) EIG(I) 340 CONTINUE WRITE(*,*) WRITE(*,*) 'MACHINE UNIT: ',RMU WRITE(*,*) WRITE(*,*) ' (THE SUBROUTINE SHOULD ALLOCATE TWO' WRITE(*,*) ' EIGENVALUES AND THEN IF THE CORRECT MACHINE' WRITE(*,*) ' UNIT IS COMPUTED THE SUBROUTINE SHOULD' WRITE(*,*) ' DETECT AN UNCONTROLLABLE SYSTEM.' WRITE(*,*) ' NOTE THAT SOME MACHINES/COMPILERS MAY NOT COMPUTE' WRITE(*,*) ' THE CORRECT MACHINE UNIT. INSTEAD THEY MAY' WRITE(*,*) ' COMPUTE A MACHINE UNIT AFFECTED BY THE LENGTH OF' WRITE(*,*) ' THE REGISTER. IN SUCH A CASE THE SUBROUTINE MAY' WRITE(*,*) ' NOT DETECT UNCONTROLLABILITY AND IT MAY CONTINUE' WRITE(*,*) ' THE ALLOCATIONS WITH AN ILL-CONDITIONED' WRITE(*,*) ' EIGENVALUE ALLOCATION PROBLEM. THIS WILL RESULT' WRITE(*,*) ' TO AN INACCURATE F AND THUS TO INACCURATE' WRITE(*,*) ' EIGENVALUES FOR (A - B*F).' WRITE(*,*) ' FOR MORE DETAILS SEE REFERENCES IN THE PAPER.' WRITE(*,*) ' IF UNCONTROLLABILITY IS DETECTED THE MATRIX' WRITE(*,*) ' (A - B*F) SHOULD HAVE EIGENVALUES 1.0, 2.0,' WRITE(*,*) ' AND THREE OTHER UNSPECIFIED EIGENVALUES.)' WRITE(*,*) C********************************************************************* C C Example 4: general system C C This program segment computes a general system C and computes the vector f such that the matrix (A - b*f) C has the eigenvalues 7.0 +- 28.0i, -15.0 +- 20.0, 63.0, and C -24.0 C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 6 DO 400 I=1,N B(I) = 1.0/DBLE(I) BB(I) = B(I) DO 400 J=1,N A(I,J) = DBLE(I) + 1.0/DBLE(J) IF(J .EQ. I+1) A(I,J) = A(I,J)*15.0 AA(I,J) = A(I,J) 400 CONTINUE L(1,1) = 7.0 L(2,1) = 28.0 L(1,2) = 7.0 L(2,2) = -28.0 L(1,3) = -15.0 L(2,3) = 20.0 L(1,4) = -15.0 L(2,4) = -20.0 L(1,5) = 63.0 L(2,5) = 0.0 L(1,6) = -24.0 L(2,6) = 0.0 WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 4 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 410 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 410 CONTINUE C C.......call DSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 420 I=1,N WRITE(*,*) F(I) DO 420 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 420 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 440 I=1,N WRITE(*,*) EIG(I) 440 CONTINUE WRITE(*,*) WRITE(*,*) ' (THE MATRIX (A - B*F) SHOULD HAVE' WRITE(*,*) ' EIGENVALUES 63.0, -24.0, 7.0 + 28.0i,' WRITE(*,*) ' 7.0 - 28.0i, -15.0 + 20.0i, and -15.0 - 20.0i,' WRITE(*,*) ' IN ANY ORDER.) THIS EXAMPLE ALLOCATES' WRITE(*,*) ' EIGENVALUES TO A GENERAL SYSTEM.' WRITE(*,*) C********************************************************************* C C Example 5: general system with C ill-conditioned eigenvalue problem C C This program segment creates a general system of order 8 C but then attempts to allocate 7 identical eigenvalues. C The resulting F vector is accurate. (This can be seen by C comparing the single-precision output of this test driver C to the double-precision version.) The eigenvalue problem C for the matrix (A - f*b) is ill-conditioned, however, C because of the identical eigenvalues. This results C in extremely inaccurate results from the DEIGVL subroutine. C The example produces a an upper-Hessenberg matrix C with non-zero subdiagonal elements and multiple C eigenvalues. This matrix is defective, and C has a Jordan canonical form with a Jordan block of the C order of the multiplicity of the eigenvalues. (In this case, C seven). A Jordan block of this size indicates a matrix C that is ill-conditioned with respect to the eigenvalue C problem. C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 8 DO 500 I=1,N IF(N-I+1 .LE. 7) THEN L(1,I) = 8.0 ELSE L(1,I) = DBLE(I) ENDIF L(2,I) = 0.0 B(I) = 1.0/DBLE(I) BB(I) = B(I) DO 500 J=1,N A(I,J) = DBLE(I) + 1.0/DBLE(J) IF(J .EQ. I+1) A(I,J) = A(I,J)*15.0 AA(I,J) = A(I,J) 500 CONTINUE WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 5 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 510 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 510 CONTINUE C C.......call DSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 520 I=1,N WRITE(*,*) F(I) DO 520 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 520 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 540 I=1,N WRITE(*,*) EIG(I) 540 CONTINUE WRITE(*,*) WRITE(*,*) ' (THIS EXAMPLE ALLOCATES THE EIGENVALUE 8.0' WRITE(*,*) ' SEVEN TIMES. HOWEVER, THE COMPUTED MULTIPLE' WRITE(*,*) ' EIGENVALUE IS NOT ACCURATE TO MACHINE PRECISION.' WRITE(*,*) ' THIS IS DUE TO THE FACT THAT UPPER HESSENBERG' WRITE(*,*) ' MATRICES WITH NON-ZERO SUBDIAGONAL ELEMENTS ARE' WRITE(*,*) ' DEFECTIVE IF THEY HAVE MULTIPLE EIGENVALUES. THIS' WRITE(*,*) ' RESULTS IN A JORDAN CANONICAL FORM WITH' WRITE(*,*) ' JORDAN BLOCKS AS LARGE AS THE MULTIPLICITY' WRITE(*,*) ' OF THE EIGENVALUES, INDICATING A MATRIX THAT' WRITE(*,*) ' IS ILL-CONDITIONED WITH RESPECT TO THE' WRITE(*,*) ' EIGENVALUE PROBLEM. THE VECTOR F IS ACCURATE,' WRITE(*,*) ' BUT THE COMPUTED EIGENVALUES FROM DEIGVL ARE' WRITE(*,*) ' NOT. THIS IS A RESULT OF THE EIGENPROBLEM' WRITE(*,*) ' OF (A - B*F), AND NOT THE POLE PLACEMENT' WRITE(*,*) ' ALGORITHM ITSELF. A COMPARISON OF THE SINGLE ' WRITE(*,*) ' AND DOUBLE PRECISION VERSIONS OF THIS PROGRAM ' WRITE(*,*) ' WILL DEMONSTRATE THIS DISTINCTION.)' WRITE(*,*) 5000 FORMAT(1X,8E15.6) STOP END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: STEST.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C C TEST DRIVER: SSEVAS C C These test segments use SEIGVL, a public domain routine C for computing the eigenvalues of a matrix. The code C for this subroutine can be found in SEIGVL.F, and was C originally found in the EISPACK subroutine library. C (Any other similar subroutine can be substituted.) C The SSEVAS subroutine requires the level 1 BLAS routines, C which are included in the file SBLAS.F, and can also be C found in the IMSL library. C C********************************************************************* C********************************************************************* PROGRAM SAMPLE PARAMETER (ISL = 25, JSL = ISL*ISL+2*ISL) REAL A(ISL,ISL),AA(ISL,ISL),B(ISL),BB(ISL) REAL F(ISL),RST(4,JSL),CST(JSL) INTEGER IST(JSL),N,IEAL COMPLEX EIG(ISL) REAL L(2,ISL),RMU LOGICAL FORM C********************************************************************* C C Example 1: system in controllability form. C C This program segment computes a system in controllability C form and computes the vector f such that the matrix (A - b*f) C has the eigenvalues (1,2,3,4,5). C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 5 DO 100 I=1,N L(1,I) = REAL(I) L(2,I) = 0.0 B(I) = 0.0 BB(I) = 0.0 DO 100 J=I-1,N IF(J .GT. 0) THEN A(I,J) = REAL(I)/REAL(J) AA(I,J) = A(I,J) ENDIF 100 CONTINUE B(1) = 1.0 BB(1) = 1.0 WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 1' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 110 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 110 CONTINUE C C.......call SSEVAS routine C FORM = .TRUE. RMU = 0.0 CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 120 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 120 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 140 I=1,N WRITE(*,*) EIG(I) 140 CONTINUE WRITE(*,*) WRITE(*,*) ' (THE MATRIX (A - B*F) SHOULD HAVE EIGENVALUES' WRITE(*,*) ' 1.0, 2.0, 3.0, 4.0 and 5.0, BUT NOT NECESSARILY' WRITE(*,*) ' IN THIS ORDER. THIS EXAMPLE USES A SYSTEM IN' WRITE(*,*) ' CONTROLLABILITY FORM.)' WRITE(*,*) C********************************************************************* C C Example 2: near-uncontrollable system in controllability C form. C C This program segment computes a system in controllability C form and attempts to compute the vector f such that the C matrix (A - b*f) has the eigenvalues (1,2,3,4,5). C The subroutine will detect that the system is near to C an uncontrollable one, and stop allocations. C C********************************************************************* C C.......create system (in controllability form) C N = 5 DO 200 I=1,N L(1,I) = REAL(I) L(2,I) = 0.0 B(I) = 0.0 BB(I) = 0.0 F(I) = 0.0 DO 200 J=1,N IF(J .GE. I-1) THEN A(I,J) = REAL(I)/REAL(J) AA(I,J) = A(I,J) ELSE A(I,J) = 0.0 AA(I,J) = 0.0 ENDIF 200 CONTINUE B(1) = 1.0 BB(1) = 1.0 C C.......make sub-diagonal element zero. C A(3,2) = 0.0 AA(3,2) = A(4,3) WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 2' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 210 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 210 CONTINUE C C.......call SSEVAS routine C FORM = .TRUE. RMU = 0.0 CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 220 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 220 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F):' WRITE(*,*) DO 240 I=1,N WRITE(*,*) EIG(I) 240 CONTINUE WRITE(*,*) WRITE(*,*) ' (SINCE THE SYSTEM (b,A) IS UNCONTROLLABLE,' WRITE(*,*) ' NO EIGENVALUES ARE ALLOCATED. THE VECTOR F' WRITE(*,*) ' SHOULD BE ZERO, AND THE EIGENVALUES OF' WRITE(*,*) ' (A - B*F) SHOULD BE THOSE OF THE MATRIX A)' WRITE(*,*) C********************************************************************* C C Example 3: near-uncontrollable system in general form. C C This program segment computes a system in general C form and attempts to compute the vector f such that the C matrix (A - b*f) has the eigenvalues (1,2,3,4,5). C The subroutine will detect that the system is near to C an uncontrollable one after allocating two eigenvalues, C and stop further allocations. The resulting matrix (A - b*f) C will thus have eigenvalues 1, 2, and three other undetermined C values. C C********************************************************************* C C.......create general system C N = 5 DO 300 I=1,N L(1,I) = REAL(I) L(2,I) = 0.0 B(I) = 1.0/REAL(I) BB(I) = B(I) F(I) = 0.0 DO 300 J=1,N A(I,J) = (REAL(I) - 3.0) / REAL(J) IF(I .EQ. J+1 .AND. J .LT. 3) + A(I,J) = A(I,J)*20.0 AA(I,J) = A(I,J) 300 CONTINUE WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 3 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 310 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 310 CONTINUE C C.......call SSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 320 I=1,N WRITE(*,*) F(I) DO 320 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 320 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F):' WRITE(*,*) DO 340 I=1,N WRITE(*,*) EIG(I) 340 CONTINUE WRITE(*,*) WRITE(*,*) 'MACHINE UNIT: ',RMU WRITE(*,*) WRITE(*,*) ' (THE SUBROUTINE SHOULD ALLOCATE TWO' WRITE(*,*) ' EIGENVALUES AND THEN IF THE CORRECT MACHINE' WRITE(*,*) ' UNIT IS COMPUTED THE SUBROUTINE SHOULD' WRITE(*,*) ' DETECT AN UNCONTROLLABLE SYSTEM.' WRITE(*,*) ' NOTE THAT SOME MACHINES/COMPILERS MAY NOT COMPUTE' WRITE(*,*) ' THE CORRECT MACHINE UNIT. INSTEAD THEY MAY' WRITE(*,*) ' COMPUTE A MACHINE UNIT AFFECTED BY THE LENGTH OF' WRITE(*,*) ' THE REGISTER. IN SUCH A CASE THE SUBROUTINE MAY' WRITE(*,*) ' NOT DETECT UNCONTROLLABILITY AND IT MAY CONTINUE' WRITE(*,*) ' THE ALLOCATIONS WITH AN ILL-CONDITIONED' WRITE(*,*) ' EIGENVALUE ALLOCATION PROBLEM. THIS WILL RESULT' WRITE(*,*) ' TO AN INACCURATE F AND THUS TO INACCURATE' WRITE(*,*) ' EIGENVALUES FOR (A - B*F).' WRITE(*,*) ' FOR MORE DETAILS SEE REFERENCES IN THE PAPER.' WRITE(*,*) ' IF UNCONTROLLABILITY IS DETECTED THE MATRIX' WRITE(*,*) ' (A - B*F) SHOULD HAVE EIGENVALUES 1.0, 2.0,' WRITE(*,*) ' AND THREE OTHER UNSPECIFIED EIGENVALUES.)' WRITE(*,*) C********************************************************************* C C Example 4: general system C C This program segment computes a general system C and computes the vector f such that the matrix (A - b*f) C has the eigenvalues 7.0 +- 28.0i, -15.0 +- 20.0, 63.0, and C -24.0 C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 6 DO 400 I=1,N B(I) = 1.0/REAL(I) BB(I) = B(I) DO 400 J=1,N A(I,J) = REAL(I) + 1.0/REAL(J) IF(J .EQ. I+1) A(I,J) = A(I,J)*15.0 AA(I,J) = A(I,J) 400 CONTINUE L(1,1) = 7.0 L(2,1) = 28.0 L(1,2) = 7.0 L(2,2) = -28.0 L(1,3) = -15.0 L(2,3) = 20.0 L(1,4) = -15.0 L(2,4) = -20.0 L(1,5) = 63.0 L(2,5) = 0.0 L(1,6) = -24.0 L(2,6) = 0.0 WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 4 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 410 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 410 CONTINUE C C.......call SSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 420 I=1,N WRITE(*,*) F(I) DO 420 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 420 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 440 I=1,N WRITE(*,*) EIG(I) 440 CONTINUE WRITE(*,*) WRITE(*,*) ' (THE MATRIX (A - B*F) SHOULD HAVE' WRITE(*,*) ' EIGENVALUES 63.0, -24.0, 7.0 + 28.0i,' WRITE(*,*) ' 7.0 - 28.0i, -15.0 + 20.0i, and -15.0 - 20.0i,' WRITE(*,*) ' IN ANY ORDER.) THIS EXAMPLE ALLOCATES' WRITE(*,*) ' EIGENVALUES TO A GENERAL SYSTEM.' WRITE(*,*) C********************************************************************* C C Example 5: general system with C ill-conditioned eigenvalue problem C C This program segment creates a general system of order 10 C but then attempts to allocate 7 identical eigenvalues. C The resulting F vector is accurate. (This can be seen by C comparing the single-precision output of this test driver C to the double-precision version.) The eigenvalue problem C for the matrix (A - f*b) is ill-conditioned, however, C because of the identical eigenvalues. This results C in extremely inaccurate results from the SEIGVL subroutine. C The example produces a an upper-Hessenberg matrix C with non-zero subdiagonal elements and multiple C eigenvalues. This matrix is defective, and C has a Jordan canonical form with a Jordan block of the C order of the multiplicity of the eigenvalues. (In this case, C seven). A Jordan block of this size indicates a matrix C that is ill-conditioned with respect to the eigenvalue C problem. C C********************************************************************* C C.......construct matrix A, vector b, and eigenvalue matrix L C N = 8 DO 500 I=1,N IF(N-I+1 .LE. 7) THEN L(1,I) = 8.0 ELSE L(1,I) = REAL(I) ENDIF L(2,I) = 0.0 B(I) = 1.0/REAL(I) BB(I) = B(I) DO 500 J=1,N A(I,J) = REAL(I) + 1.0/REAL(J) IF(J .EQ. I+1) A(I,J) = A(I,J)*15.0 AA(I,J) = A(I,J) 500 CONTINUE WRITE(*,*) WRITE(*,*) '*********' WRITE(*,*) 'EXAMPLE 5 ' WRITE(*,*) '*********' WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX A:' WRITE(*,*) DO 510 I=1,N WRITE(*,5000) (A(I,J),J=1,N) 510 CONTINUE C C.......call SSEVAS routine C FORM = .FALSE. RMU = 0.0 CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 520 I=1,N WRITE(*,*) F(I) DO 520 J=1,N AA(I,J) = AA(I,J) - BB(I)*F(J) 520 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 540 I=1,N WRITE(*,*) EIG(I) 540 CONTINUE WRITE(*,*) WRITE(*,*) ' (THIS EXAMPLE ALLOCATES THE EIGENVALUE 8.0' WRITE(*,*) ' SEVEN TIMES. HOWEVER, THE COMPUTED MULTIPLE' WRITE(*,*) ' EIGENVALUE IS NOT ACCURATE TO MACHINE PRECISION.' WRITE(*,*) ' THIS IS DUE TO THE FACT THAT UPPER HESSENBERG' WRITE(*,*) ' MATRICES WITH NON-ZERO SUBDIAGONAL ELEMENTS ARE' WRITE(*,*) ' DEFECTIVE IF THEY HAVE MULTIPLE EIGENVALUES. THIS' WRITE(*,*) ' RESULTS IN A JORDAN CANONICAL FORM WITH' WRITE(*,*) ' JORDAN BLOCKS AS LARGE AS THE MULTIPLICITY' WRITE(*,*) ' OF THE EIGENVALUES, INDICATING A MATRIX THAT' WRITE(*,*) ' IS ILL-CONDITIONED WITH RESPECT TO THE' WRITE(*,*) ' EIGENVALUE PROBLEM. THE VECTOR F IS ACCURATE,' WRITE(*,*) ' BUT THE COMPUTED EIGENVALUES FROM DEIGVL ARE' WRITE(*,*) ' NOT. THIS IS A RESULT OF THE EIGENPROBLEM' WRITE(*,*) ' OF (A - B*F), AND NOT THE POLE PLACEMENT' WRITE(*,*) ' ALGORITHM ITSELF. A COMPARISON OF THE SINGLE ' WRITE(*,*) ' AND DOUBLE PRECISION VERSIONS OF THIS PROGRAM ' WRITE(*,*) ' WILL DEMONSTRATE THIS DISTINCTION.)' WRITE(*,*) 5000 FORMAT(1X,8E15.6) STOP END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: SBLAS.F XXXXXXXXXXXXXXXXXXXXX C******************************************************** C C Level 1 BLAS routines used in SSEVAS C C******************************************************** INTEGER FUNCTION ISAMAX(N,SX,INCX) C C finds the index of element having max. absolute value. C jack dongarra, linpack, 3/11/78. C modified to correct problem with negative increments, 9/29/88. C REAL SX(1),SMAX INTEGER I,INCX,IX,N C ISAMAX = 0 IF( N .LT. 1 ) RETURN ISAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C IX = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 SMAX = ABS(SX(IX)) IX = IX + INCX DO 10 I = 2,N IF(ABS(SX(IX)).LE.SMAX) GO TO 5 ISAMAX = I SMAX = ABS(SX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C code for increment equal to 1 C 20 SMAX = ABS(SX(1)) DO 30 I = 2,N IF(ABS(SX(I)).LE.SMAX) GO TO 30 ISAMAX = I SMAX = ABS(SX(I)) 30 CONTINUE RETURN END REAL FUNCTION SASUM(N,SX,INCX) C C takes the sum of the absolute values. C uses unrolled loops for increment equal to one. C jack dongarra, linpack, 3/11/78. C modified to correct problem with negative increments, 9/29/88. C REAL SX(1),STEMP INTEGER I,IX,INCX,M,MP1,N C SASUM = 0.0E0 STEMP = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C IX = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 DO 10 I = 1,N STEMP = STEMP + ABS(SX(IX)) IX = IX + INCX 10 CONTINUE SASUM = STEMP RETURN C C code for increment equal to 1 C C C clean-up loop C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M STEMP = STEMP + ABS(SX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 STEMP = STEMP + ABS(SX(I)) + ABS(SX(I + 1)) + ABS(SX(I + 2)) * + ABS(SX(I + 3)) + ABS(SX(I + 4)) + ABS(SX(I + 5)) 50 CONTINUE 60 SASUM = STEMP RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C C constant times a vector plus a vector. C uses unrolled loop for increments equal to one. C jack dongarra, linpack, 3/11/78. C REAL SX(1),SY(1),SA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (SA .EQ. 0.0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SY(I) + SA*SX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 SY(I) = SY(I) + SA*SX(I) SY(I + 1) = SY(I + 1) + SA*SX(I + 1) SY(I + 2) = SY(I + 2) + SA*SX(I + 2) SY(I + 3) = SY(I + 3) + SA*SX(I + 3) 50 CONTINUE RETURN END SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) C C copies a vector, x, to a vector, y. C uses unrolled loops for increments equal to 1. C jack dongarra, linpack, 3/11/78. C REAL SX(1),SY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SY(I) = SX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 SY(I) = SX(I) SY(I + 1) = SX(I + 1) SY(I + 2) = SX(I + 2) SY(I + 3) = SX(I + 3) SY(I + 4) = SX(I + 4) SY(I + 5) = SX(I + 5) SY(I + 6) = SX(I + 6) 50 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C C forms the dot product of two vectors. C uses unrolled loops for increments equal to one. C jack dongarra, linpack, 3/11/78. C REAL SX(1),SY(1),STEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C STEMP = 0.0E0 SDOT = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE SDOT = STEMP RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M STEMP = STEMP + SX(I)*SY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 STEMP = STEMP + SX(I)*SY(I) + SX(I + 1)*SY(I + 1) + * SX(I+2)*SY(I+2) + SX(I + 3)*SY(I + 3) + SX(I + 4)*SY(I + 4) 50 CONTINUE 60 SDOT = STEMP RETURN END REAL FUNCTION SNRM2 ( N, SX, INCX) INTEGER NEXT REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE DATA ZERO, ONE /0.0E0, 1.0E0/ C C euclidean norm of the n-vector stored in sx() with storage C increment incx . C if n .le. 0 return with result = 0. C if n .ge. 1 then incx must be .ge. 1 C C c.l.lawson, 1978 jan 08 C C four phase method using two built-in constants that are C hopefully applicable to all machines. C cutlo = maximum of sqrt(u/eps) over all known machines. C cuthi = minimum of sqrt(v) over all known machines. C where C eps = smallest no. such that eps + 1. .gt. 1. C u = smallest positive no. (underflow limit) C v = largest no. (overflow limit) C C brief outline of algorithm.. C C phase 1 scans zero components. C move to phase 2 when a component is nonzero and .le. cutlo C move to phase 3 when a component is .gt. cutlo C move to phase 4 when a component is .ge. cuthi/m C where m = n for x() real and m = 2*n for complex. C C values for cutlo and cuthi.. C from the environmental parameters listed in the imsl converter C document the limiting values are as follows.. C cutlo, s.p. u/eps = 2**(-102) for honeywell.close seconds are C univac and dec at 2**(-103) C thus cutlo = 2**(-51) = 4.44089e-16 C cuthi, s.p. v = 2**127 for univac, honeywell, and dec. C thus cuthi = 2**(63.5) = 1.30438e19 C cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. C thus cutlo = 2**(-33.5) = 8.23181d-11 C cuthi, d.p. same as s.p. cuthi = 1.30438d19 C data cutlo, cuthi / 8.232d-11, 1.304d19 / C data cutlo, cuthi / 4.441e-16, 1.304e19 / DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C IF(N .GT. 0) GO TO 10 SNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C begin main loop I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C phase 1. sum is zero C 50 IF( SX(I) .EQ. ZERO) GO TO 200 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 C C prepare for phase 2. ASSIGN 70 TO NEXT GO TO 105 C C prepare for phase 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / SX(I)) / SX(I) 105 XMAX = ABS(SX(I)) GO TO 115 C C phase 2. sum is small. C scale to avoid destructive underflow. C 70 IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75 C C common code for phases 2 and 4. C in phase 4 sum is large. scale to avoid overflow. C 110 IF( ABS(SX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / SX(I))**2 XMAX = ABS(SX(I)) GO TO 200 C 115 SUM = SUM + (SX(I)/XMAX)**2 GO TO 200 C C C prepare for phase 3. C 75 SUM = (SUM * XMAX) * XMAX C C C for real or d.p. set hitest = cuthi/n C for complex set hitest = cuthi/(2*n) C 85 HITEST = CUTHI/FLOAT( N ) C C phase 3. sum is mid-range. no scaling. C DO 95 J =I,NN,INCX IF(ABS(SX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + SX(J)**2 SNRM2 = SQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C end of main loop. C C compute square root and adjust for scaling. C SNRM2 = XMAX * SQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE SSCAL(N,SA,SX,INCX) C C scales a vector by a constant. C uses unrolled loops for increment equal to 1. C jack dongarra, linpack, 3/11/78. C modified to correct problem with negative increments, 9/29/88. C REAL SA,SX(1) INTEGER I,IX,INCX,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C IX = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 DO 10 I = 1,N SX(IX) = SA*SX(IX) IX = IX + INCX 10 CONTINUE RETURN C C code for increment equal to 1 C C C clean-up loop C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M SX(I) = SA*SX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 SX(I) = SA*SX(I) SX(I + 1) = SA*SX(I + 1) SX(I + 2) = SA*SX(I + 2) SX(I + 3) = SA*SX(I + 3) SX(I + 4) = SA*SX(I + 4) 50 CONTINUE RETURN END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: SEIGVL.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C********************************************************************* C C Title: SEIGVL C C Purpose: Finds the eigenvalues of a real-valued C matrix AA. Uses BALANC, HQR and ELMHES, from C the EISPACK subroutine library. (These C routines are included in this file.) C C C Input Arguments: C ---------------- C C n integer Row dimension of matrix AA C C AA real*4(adim,n) Input matrix AA C C llda integer Row dimension of array AA C C C C Output Arguments: C -------------- C C eig complex(n) Desired eigenvalues C C C C Subroutines Used: C ----------------- C C BALANC, HQR, ELMHES C C********************************************************************* C********************************************************************* C********************************************************************* SUBROUTINE SEIGVL(N,AA,LLDA,EIG) PARAMETER (MAXN = 25) C C........input arguments C INTEGER N,LLDA REAL AA(LLDA,*) COMPLEX EIG(MAXN) C C.......local variables C REAL A(MAXN,MAXN) REAL REIG(MAXN), IMEIG(MAXN), SCALE(MAXN) INTEGER INT(MAXN), LOW, IGH, IERR C C.......check size of input system C IF(N .GT. MAXN) THEN WRITE(*,*) 'SEIGVL: INPUT SYSTEM TOO LARGE.' STOP ENDIF C C.......copy input matrix to working matrix C DO 40 I=1,N DO 40 J=1,N A(I,J) = AA(I,J) 40 CONTINUE LDA = MAXN C C.......convert general form to Upper Hessenberg form C and find eigenvalues. C CALL BALANC( LDA, N, A, LOW, IGH, SCALE) CALL ELMHES( LDA, N, LOW, IGH, A, INT) CALL HQR( LDA, N, LOW, IGH, A, REIG, IMEIG, IERR) DO 100 I=1,N EIG(I) = CMPLX(REIG(I),IMEIG(I)) 100 CONTINUE RETURN END C The following subroutines from EISPACK C Ccccccccccccccccccccccccccccccccccccccccccccccccccc Ccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Mon Apr 1 12:47:07 EST 1991 *** C SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL A(NM,N),SCALE(N) REAL C,F,G,R,S,B2,RADIX LOGICAL NOCONV C C this subroutine is a translation of the algol procedure balance, C num. math. 13, 293-304(1969) by parlett and reinsch. C handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). C C this subroutine balances a real matrix and isolates C eigenvalues whenever possible. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C a contains the input matrix to be balanced. C C on output C C a contains the balanced matrix. C C low and igh are two integers such that a(i,j) C is equal to zero if C (1) i is greater than j and C (2) j=1,...,low-1 or i=igh+1,...,n. C C scale contains information determining the C permutations and scaling factors used. C C suppose that the principal submatrix in rows low through igh C has been balanced, that p(j) denotes the index interchanged C with j during the permutation step, and that the elements C of the diagonal matrix used are denoted by d(i,j). then C scale(j) = p(j), for j = 1,...,low-1 C = d(j,j), j = low,...,igh C = p(j) j = igh+1,...,n. C the order in which the interchanges are made is n to igh+1, C then 1 to low-1. C C note that 1 is returned for igh if igh is zero formally. C C the algol procedure exc contained in balance appears in C balanc in line. (note that the algol roles of identifiers C k,l have been reversed.) C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated august 1983. C C ------------------------------------------------------------------ C RADIX = 16.0E0 C B2 = RADIX * RADIX K = 1 L = N GO TO 100 C .......... in-line procedure for row and C column exchange .......... 20 SCALE(M) = J IF (J .EQ. M) GO TO 50 C DO 30 I = 1, L F = A(I,J) A(I,J) = A(I,M) A(I,M) = F 30 CONTINUE C DO 40 I = K, N F = A(J,I) A(J,I) = A(M,I) A(M,I) = F 40 CONTINUE C 50 GO TO (80,130), IEXC C .......... search for rows isolating an eigenvalue C and push them down .......... 80 IF (L .EQ. 1) GO TO 280 L = L - 1 C .......... for j=l step -1 until 1 do -- .......... 100 DO 120 JJ = 1, L J = L + 1 - JJ C DO 110 I = 1, L IF (I .EQ. J) GO TO 110 IF (A(J,I) .NE. 0.0E0) GO TO 120 110 CONTINUE C M = L IEXC = 1 GO TO 20 120 CONTINUE C GO TO 140 C .......... search for columns isolating an eigenvalue C and push them left .......... 130 K = K + 1 C 140 DO 170 J = K, L C DO 150 I = K, L IF (I .EQ. J) GO TO 150 IF (A(I,J) .NE. 0.0E0) GO TO 170 150 CONTINUE C M = K IEXC = 2 GO TO 20 170 CONTINUE C .......... now balance the submatrix in rows k to l .......... DO 180 I = K, L 180 SCALE(I) = 1.0E0 C .......... iterative loop for norm reduction .......... 190 NOCONV = .FALSE. C DO 270 I = K, L C = 0.0E0 R = 0.0E0 C DO 200 J = K, L IF (J .EQ. I) GO TO 200 C = C + ABS(A(J,I)) R = R + ABS(A(I,J)) 200 CONTINUE C .......... guard against zero c or r due to underflow .......... IF (C .EQ. 0.0E0 .OR. R .EQ. 0.0E0) GO TO 270 G = R / RADIX F = 1.0E0 S = C + R 210 IF (C .GE. G) GO TO 220 F = F * RADIX C = C * B2 GO TO 210 220 G = R * RADIX 230 IF (C .LT. G) GO TO 240 F = F / RADIX C = C / B2 GO TO 230 C .......... now balance .......... 240 IF ((C + R) / F .GE. 0.95E0 * S) GO TO 270 G = 1.0E0 / F SCALE(I) = SCALE(I) * F NOCONV = .TRUE. C DO 250 J = K, N 250 A(I,J) = A(I,J) * G C DO 260 J = 1, L 260 A(J,I) = A(J,I) * F C 270 CONTINUE C IF (NOCONV) GO TO 190 C 280 LOW = K IGH = L RETURN END C Ccccccccccccccccccccccccccccccccccccccccccccccccccc Ccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Mon Apr 1 12:47:08 EST 1991 *** C SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 REAL A(NM,N) REAL X,Y INTEGER INT(IGH) C C this subroutine is a translation of the algol procedure elmhes, C num. math. 12, 349-368(1968) by martin and wilkinson. C handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). C C given a real general matrix, this subroutine C reduces a submatrix situated in rows and columns C low through igh to upper hessenberg form by C stabilized elementary similarity transformations. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C low and igh are integers determined by the balancing C subroutine balanc. if balanc has not been used, C set low=1, igh=n. C C a contains the input matrix. C C on output C C a contains the hessenberg matrix. the multipliers C which were used in the reduction are stored in the C remaining triangle under the hessenberg matrix. C C int contains information on the rows and columns C interchanged in the reduction. C only elements low through igh are used. C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated august 1983. C C ------------------------------------------------------------------ C LA = IGH - 1 KP1 = LOW + 1 IF (LA .LT. KP1) GO TO 200 C DO 180 M = KP1, LA MM1 = M - 1 X = 0.0E0 I = M C DO 100 J = M, IGH IF (ABS(A(J,MM1)) .LE. ABS(X)) GO TO 100 X = A(J,MM1) I = J 100 CONTINUE C INT(M) = I IF (I .EQ. M) GO TO 130 C .......... interchange rows and columns of a .......... DO 110 J = MM1, N Y = A(I,J) A(I,J) = A(M,J) A(M,J) = Y 110 CONTINUE C DO 120 J = 1, IGH Y = A(J,I) A(J,I) = A(J,M) A(J,M) = Y 120 CONTINUE C .......... end interchange .......... 130 IF (X .EQ. 0.0E0) GO TO 180 MP1 = M + 1 C DO 160 I = MP1, IGH Y = A(I,MM1) IF (Y .EQ. 0.0E0) GO TO 160 Y = Y / X A(I,MM1) = Y C DO 140 J = M, N 140 A(I,J) = A(I,J) - Y * A(M,J) C DO 150 J = 1, IGH 150 A(J,M) = A(J,M) + Y * A(J,I) C 160 CONTINUE C 180 CONTINUE C 200 RETURN END C Ccccccccccccccccccccccccccccccccccccccccccccccccccc Ccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Mon Apr 1 12:47:10 EST 1991 *** C SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) C RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG) C INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR REAL H(NM,N),WR(N),WI(N) REAL P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2 LOGICAL NOTLAS C C this subroutine is a translation of the algol procedure hqr, C num. math. 14, 219-231(1970) by martin, peters, and wilkinson. C handbook for auto. comp., vol.ii-linear algebra, 359-371(1971). C C this subroutine finds the eigenvalues of a real C upper hessenberg matrix by the qr method. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C low and igh are integers determined by the balancing C subroutine balanc. if balanc has not been used, C set low=1, igh=n. C C h contains the upper hessenberg matrix. information about C the transformations used in the reduction to hessenberg C form by elmhes or orthes, if performed, is stored C in the remaining triangle under the hessenberg matrix. C C on output C C h has been destroyed. therefore, it must be saved C before calling hqr if subsequent calculation and C back transformation of eigenvectors is to be performed. C C wr and wi contain the real and imaginary parts, C respectively, of the eigenvalues. the eigenvalues C are unordered except that complex conjugate pairs C of values appear consecutively with the eigenvalue C having the positive imaginary part first. if an C error exit is made, the eigenvalues should be correct C for indices ierr+1,...,n. C C ierr is set to C zero for normal return, C j if the limit of 30*n iterations is exhausted C while the j-th eigenvalue is being sought. C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated september 1989. C C ------------------------------------------------------------------ C IERR = 0 NORM = 0.0E0 K = 1 C .......... store roots isolated by balanc C and compute matrix norm .......... DO 50 I = 1, N C DO 40 J = K, N 40 NORM = NORM + ABS(H(I,J)) C K = I IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 WR(I) = H(I,I) WI(I) = 0.0E0 50 CONTINUE C EN = IGH T = 0.0E0 ITN = 30*N C .......... search for next eigenvalues .......... 60 IF (EN .LT. LOW) GO TO 1001 ITS = 0 NA = EN - 1 ENM2 = NA - 1 C .......... look for single small sub-diagonal element C for l=en step -1 until low do -- .......... 70 DO 80 LL = LOW, EN L = EN + LOW - LL IF (L .EQ. LOW) GO TO 100 S = ABS(H(L-1,L-1)) + ABS(H(L,L)) IF (S .EQ. 0.0E0) S = NORM TST1 = S TST2 = TST1 + ABS(H(L,L-1)) IF (TST2 .EQ. TST1) GO TO 100 80 CONTINUE C .......... form shift .......... 100 X = H(EN,EN) IF (L .EQ. EN) GO TO 270 Y = H(NA,NA) W = H(EN,NA) * H(NA,EN) IF (L .EQ. NA) GO TO 280 IF (ITN .EQ. 0) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 C .......... form exceptional shift .......... T = T + X C DO 120 I = LOW, EN 120 H(I,I) = H(I,I) - X C S = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75E0 * S Y = X W = -0.4375E0 * S * S 130 ITS = ITS + 1 ITN = ITN - 1 C .......... look for two consecutive small C sub-diagonal elements. C for m=en-2 step -1 until l do -- .......... DO 140 MM = L, ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R * S - W) / H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = ABS(P) + ABS(Q) + ABS(R) P = P / S Q = Q / S R = R / S IF (M .EQ. L) GO TO 150 TST1 = ABS(P)*(ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1))) TST2 = TST1 + ABS(H(M,M-1))*(ABS(Q) + ABS(R)) IF (TST2 .EQ. TST1) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.0E0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0E0 160 CONTINUE C .......... double qr step involving rows l to en and C columns m to en .......... DO 260 K = M, NA NOTLAS = K .NE. NA IF (K .EQ. M) GO TO 170 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0E0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X .EQ. 0.0E0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = SIGN(SQRT(P*P+Q*Q+R*R),P) IF (K .EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 190 P = P + S X = P / S Y = Q / S ZZ = R / S Q = Q / P R = R / P IF (NOTLAS) GO TO 225 C .......... row modification .......... DO 200 J = K, EN P = H(K,J) + Q * H(K+1,J) H(K,J) = H(K,J) - P * X H(K+1,J) = H(K+1,J) - P * Y 200 CONTINUE C J = MIN0(EN,K+3) C .......... column modification .......... DO 210 I = L, J P = X * H(I,K) + Y * H(I,K+1) H(I,K) = H(I,K) - P H(I,K+1) = H(I,K+1) - P * Q 210 CONTINUE GO TO 255 225 CONTINUE C .......... row modification .......... DO 230 J = K, EN P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) H(K,J) = H(K,J) - P * X H(K+1,J) = H(K+1,J) - P * Y H(K+2,J) = H(K+2,J) - P * ZZ 230 CONTINUE C J = MIN0(EN,K+3) C .......... column modification .......... DO 240 I = L, J P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) H(I,K) = H(I,K) - P H(I,K+1) = H(I,K+1) - P * Q H(I,K+2) = H(I,K+2) - P * R 240 CONTINUE 255 CONTINUE C 260 CONTINUE C GO TO 70 C .......... one root found .......... 270 WR(EN) = X + T WI(EN) = 0.0E0 EN = NA GO TO 60 C .......... two roots found .......... 280 P = (Y - X) / 2.0E0 Q = P * P + W ZZ = SQRT(ABS(Q)) X = X + T IF (Q .LT. 0.0E0) GO TO 320 C .......... real pair .......... ZZ = P + SIGN(ZZ,P) WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ .NE. 0.0E0) WR(EN) = X - W / ZZ WI(NA) = 0.0E0 WI(EN) = 0.0E0 GO TO 330 C .......... complex pair .......... 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 C .......... set error -- all eigenvalues have not C converged after 30*n iterations .......... 1000 IERR = EN 1001 RETURN END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: SFLRD.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C C DRIVER PROGRAM FOR SSEVAS C C This program prompts the user for a file name that C contains the input, reads the file, allocates the desired C eigenvalues, and finds the eigenvalues of the C matrix (A - b*f). The input values, output vector f, C and the resulting eigenvalues are printed. C C INPUT FORMAT: C ------------- C C 1) Line 1 of the file contains an integer in free format C showing the size of the system. C C 2) Line 2 of the file contains a logical value in free format C showing the form of the system. A value of "T" C indicates a system in controllability form, whereas C a value of "F" indicates a system in general form. C C 3) Line 3 of the file contains a real value in free format C giving the machine unit used for testing for C controllability. If the value 0.0 is given C the subroutine calculates the machine unit. C C 4) The matrix A is read next, one element per line, C in free format. The matrix is read by row. C C 5) The vector b is read next, one element per line, C in free format. C C 6) The desired eigenvalues are read next. Each line C contains an ordered pair of two real, free format C values which represent the real and the C imaginary part of a complex eigenvalue. C Real eigenvalues should obviously be represented C by an ordered pair with a zero imaginary part. C Complex conjugates should appear both, and C contiguously. Real eigenvalues should appear in C groups of even numbers, unless they are at the end C of the list. A good way to list the eigenvalues C would be, to first list all the complex eigenvalues C followed be the real eigenvalues. C C A typical example could be the following: C C 4 C T C 0.0 C 1 For A= 1 2 3 4 C 2 5 6 7 8 C 3 0 9 1 2 C 4 0 0 3 4 C 5 C 6 C 7 C 8 C 0 C 9 C 1 C 2 C 0 C 0 C 3 C 4 C 10 For b= 10 C 0 0 C 0 0 C 0 0 C 1 2 C 1 -2 C 3 0 C 4 0 C C Note that in this test driver program we have set a limit C of 25 for the size of the matrix A. This obviously C can be altered by the user. C C C********************************************************************* C********************************************************************* PROGRAM SAMPLE PARAMETER (ISL = 25, JSL = ISL*ISL+2*ISL) REAL A(ISL,ISL),AA(ISL,ISL),B(ISL),BB(ISL) REAL F(ISL),RST(4,JSL),CST(JSL) INTEGER IST(JSL),N,IEAL COMPLEX EIG(ISL) REAL L(2,ISL), RMU LOGICAL FORM CHARACTER*20 FNAME C C.......read matrix A, vector b, and eigenvalue matrix L C WRITE(*,*) WRITE(*,*) '*****************************' WRITE(*,*) ' SSEVAS TEST PROGRAM' WRITE(*,*) '*****************************' WRITE(*,*) WRITE(*,4000) 4000 FORMAT(1X,'DATA FILE: ',/) READ(*,7000) FNAME 7000 FORMAT(A) OPEN(5, FILE=FNAME, STATUS='OLD') READ(5,*) N IF(N .GT. ISL) THEN WRITE(*,*) 'INPUT SYSTEM TOO LARGE.' STOP ENDIF READ(5,*) FORM READ(5,*) RMU DO 100 I=1,N DO 100 J=1,N READ(5,*) A(I,J) AA(I,J) = A(I,J) 100 CONTINUE DO 110 I=1,N READ(5,*) B(I) BB(I) = B(I) 110 CONTINUE DO 120 I=1,N READ(5,*) L(1,I), L(2,I) 120 CONTINUE CLOSE(5) C C.......Echo input C WRITE(*,*) WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX (b,A):' WRITE(*,*) DO 130 I=1,N WRITE(*,5000) B(I),(A(I,J),J=1,N) 130 CONTINUE WRITE(*,*) WRITE(*,*) 'EIGENVALUES TO BE ALLOCATED:' WRITE(*,*) DO 135 I=1,N WRITE(*,5000) L(1,I),L(2,I) 135 CONTINUE C C.......call SSEVAS routine C CALL SSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 140 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 140 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL SEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 150 I=1,N WRITE(*,*) EIG(I) 150 CONTINUE 5000 FORMAT(1X,10E15.6) STOP END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: DBLAS.F XXXXXXXXXXXXXXXXXXXXX C******************************************************** C C Level 1 BLAS routines used by DSEVAS C C******************************************************** DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) C C takes the sum of the absolute values. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DX(1),DTEMP INTEGER I,INCX,M,MP1,N,NINCX C DASUM = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DTEMP = DTEMP + DABS(DX(I)) 10 CONTINUE DASUM = DTEMP RETURN C C code for increment equal to 1 C C C clean-up loop C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 DTEMP = DTEMP+DABS(DX(I))+DABS(DX(I + 1))+DABS(DX(I + 2)) * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 50 CONTINUE 60 DASUM = DTEMP RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C constant times a vector plus a vector. C uses unrolled loops for increments equal to one. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C copies a vector, x, to a vector, y. C uses unrolled loops for increments equal to one. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C forms the dot product of two vectors. C uses unrolled loops for increments equal to one. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C code for unequal increments or equal increments C not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C code for both increments equal to 1 C C C clean-up loop C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER NEXT DOUBLE PRECISION DX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C euclidean norm of the n-vector stored in dx() with storage C increment incx . C if n .le. 0 return with result = 0. C if n .ge. 1 then incx must be .ge. 1 C C c.l.lawson, 1978 jan 08 C C four phase method using two built-in constants that are C hopefully applicable to all machines. C cutlo = maximum of dsqrt(u/eps) over all known machines. C cuthi = minimum of dsqrt(v) over all known machines. C where C eps = smallest no. such that eps + 1. .gt. 1. C u = smallest positive no. (underflow limit) C v = largest no. (overflow limit) C C brief outline of algorithm.. C C phase 1 scans zero components. C move to phase 2 when a component is nonzero and .le. cutlo C move to phase 3 when a component is .gt. cutlo C move to phase 4 when a component is .ge. cuthi/m C where m = n for x() real and m = 2*n for complex. C C values for cutlo and cuthi.. C from the environmental parameters listed in the imsl converter C document the limiting values are as follows.. C cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are C univac and dec at 2**(-103) C thus cutlo = 2**(-51) = 4.44089e-16 C cuthi, s.p. v = 2**127 for univac, honeywell, and dec. C thus cuthi = 2**(63.5) = 1.30438e19 C cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. C thus cutlo = 2**(-33.5) = 8.23181d-11 C cuthi, d.p. same as s.p. cuthi = 1.30438d19 C data cutlo, cuthi / 8.232d-11, 1.304d19 / C data cutlo, cuthi / 4.441e-16, 1.304e19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C begin main loop I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C phase 1. sum is zero C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C prepare for phase 2. ASSIGN 70 TO NEXT GO TO 105 C C prepare for phase 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C phase 2. sum is small. C scale to avoid destructive underflow. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C common code for phases 2 and 4. C in phase 4 sum is large. scale to avoid overflow. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C prepare for phase 3. C 75 SUM = (SUM * XMAX) * XMAX C C C for real or d.p. set hitest = cuthi/n C for complex set hitest = cuthi/(2*n) C 85 HITEST = CUTHI/FLOAT( N ) C C phase 3. sum is mid-range. no scaling. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C end of main loop. C C compute square root and adjust for scaling. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C C scales a vector by a constant. C uses unrolled loops for increment equal to one. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C code for increment equal to 1 C C C clean-up loop C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END INTEGER FUNCTION IDAMAX(N,DX,INCX) C C finds the index of element having max. absolute value. C jack dongarra, linpack, 3/11/78. C DOUBLE PRECISION DX(1),DMAX INTEGER I,INCX,IX,N C IDAMAX = 0 IF( N .LT. 1 ) RETURN IDAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C code for increment not equal to 1 C IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(DABS(DX(IX)).LE.DMAX) GO TO 5 IDAMAX = I DMAX = DABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C code for increment equal to 1 C 20 DMAX = DABS(DX(1)) DO 30 I = 2,N IF(DABS(DX(I)).LE.DMAX) GO TO 30 IDAMAX = I DMAX = DABS(DX(I)) 30 CONTINUE RETURN END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: DEIGVL.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C********************************************************************* C C Title: DEIGVL C C Purpose: Finds the eigenvalues of a real-valued C matrix AA. Uses BALANC, HQR, and ELMHES, C from the EISPACK subroutine library. C (These subroutines are included in this C file.) C C C Input Arguments: C ---------------- C C n integer Row dimension of matrix AA C C AA real*8(adim,n) Input matrix AA C C llda integer Row dimension of array AA C C C C Output Arguments: C -------------- C C eig complex*16(n) Desired eigenvalues C C C C Subroutines Used: C ----------------- C C BALANC, HQR, ELMHES C C********************************************************************* C********************************************************************* C********************************************************************* SUBROUTINE DEIGVL(N,AA,LLDA,EIG) PARAMETER (MAXN = 25) C C........input arguments C INTEGER N,LLDA DOUBLE PRECISION AA(LLDA,*) COMPLEX*16 EIG(MAXN) C C.......local variables C DOUBLE PRECISION A(MAXN,MAXN) DOUBLE PRECISION REIG(MAXN), IMEIG(MAXN), SCALE(MAXN) INTEGER INT(MAXN),LOW,IGH,IERR C C.......check for size of input system C IF(N .GT. MAXN) THEN WRITE(*,*) 'DEIGVL: INPUT SYSTEM TOO LARGE.' STOP ENDIF C C.......copy input matrix to working matrix C DO 40 I=1,N DO 40 J=1,N A(I,J) = AA(I,J) 40 CONTINUE LDA = MAXN C C.......convert general form to Upper Hessenberg form C and find eigenvalues. C CALL BALANC( LDA, N, A, LOW, IGH, SCALE) CALL ELMHES( LDA, N, LOW, IGH, A, INT) CALL HQR( LDA, N, LOW, IGH, A, REIG, IMEIG, IERR) DO 100 I=1,N EIG(I) = DCMPLX(REIG(I),IMEIG(I)) 100 CONTINUE RETURN END C The following subroutines from EISPACK C Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Wed Mar 20 07:57:52 EST 1991 *** C SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC DOUBLE PRECISION A(NM,N),SCALE(N) DOUBLE PRECISION C,F,G,R,S,B2,RADIX LOGICAL NOCONV C C this subroutine is a translation of the algol procedure balance, C num. math. 13, 293-304(1969) by parlett and reinsch. C handbook for auto. comp., vol.ii-linear algebra, 315-326(1971). C C this subroutine balances a real matrix and isolates C eigenvalues whenever possible. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C a contains the input matrix to be balanced. C C on output C C a contains the balanced matrix. C C low and igh are two integers such that a(i,j) C is equal to zero if C (1) i is greater than j and C (2) j=1,...,low-1 or i=igh+1,...,n. C C scale contains information determining the C permutations and scaling factors used. C C suppose that the principal submatrix in rows low through igh C has been balanced, that p(j) denotes the index interchanged C with j during the permutation step, and that the elements C of the diagonal matrix used are denoted by d(i,j). then C scale(j) = p(j), for j = 1,...,low-1 C = d(j,j), j = low,...,igh C = p(j) j = igh+1,...,n. C the order in which the interchanges are made is n to igh+1, C then 1 to low-1. C C note that 1 is returned for igh if igh is zero formally. C C the algol procedure exc contained in balance appears in C balanc in line. (note that the algol roles of identifiers C k,l have been reversed.) C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated august 1983. C C ------------------------------------------------------------------ C RADIX = 16.0D0 C B2 = RADIX * RADIX K = 1 L = N GO TO 100 C .......... in-line procedure for row and C column exchange .......... 20 SCALE(M) = J IF (J .EQ. M) GO TO 50 C DO 30 I = 1, L F = A(I,J) A(I,J) = A(I,M) A(I,M) = F 30 CONTINUE C DO 40 I = K, N F = A(J,I) A(J,I) = A(M,I) A(M,I) = F 40 CONTINUE C 50 GO TO (80,130), IEXC C .......... search for rows isolating an eigenvalue C and push them down .......... 80 IF (L .EQ. 1) GO TO 280 L = L - 1 C .......... for j=l step -1 until 1 do -- .......... 100 DO 120 JJ = 1, L J = L + 1 - JJ C DO 110 I = 1, L IF (I .EQ. J) GO TO 110 IF (A(J,I) .NE. 0.0D0) GO TO 120 110 CONTINUE C M = L IEXC = 1 GO TO 20 120 CONTINUE C GO TO 140 C .......... search for columns isolating an eigenvalue C and push them left .......... 130 K = K + 1 C 140 DO 170 J = K, L C DO 150 I = K, L IF (I .EQ. J) GO TO 150 IF (A(I,J) .NE. 0.0D0) GO TO 170 150 CONTINUE C M = K IEXC = 2 GO TO 20 170 CONTINUE C .......... now balance the submatrix in rows k to l .......... DO 180 I = K, L 180 SCALE(I) = 1.0D0 C .......... iterative loop for norm reduction .......... 190 NOCONV = .FALSE. C DO 270 I = K, L C = 0.0D0 R = 0.0D0 C DO 200 J = K, L IF (J .EQ. I) GO TO 200 C = C + DABS(A(J,I)) R = R + DABS(A(I,J)) 200 CONTINUE C .......... guard against zero c or r due to underflow .......... IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270 G = R / RADIX F = 1.0D0 S = C + R 210 IF (C .GE. G) GO TO 220 F = F * RADIX C = C * B2 GO TO 210 220 G = R * RADIX 230 IF (C .LT. G) GO TO 240 F = F / RADIX C = C / B2 GO TO 230 C .......... now balance .......... 240 IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270 G = 1.0D0 / F SCALE(I) = SCALE(I) * F NOCONV = .TRUE. C DO 250 J = K, N 250 A(I,J) = A(I,J) * G C DO 260 J = 1, L 260 A(J,I) = A(J,I) * F C 270 CONTINUE C IF (NOCONV) GO TO 190 C 280 LOW = K IGH = L RETURN END C Cccccccccccccccccccccccccccccccccccccccccccccccccccccccc Cccccccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Tue Mar 19 12:19:40 EST 1991 *** C SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT) C INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1 DOUBLE PRECISION A(NM,N) DOUBLE PRECISION X,Y INTEGER INT(IGH) C C this subroutine is a translation of the algol procedure elmhes, C num. math. 12, 349-368(1968) by martin and wilkinson. C handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). C C given a real general matrix, this subroutine C reduces a submatrix situated in rows and columns C low through igh to upper hessenberg form by C stabilized elementary similarity transformations. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C low and igh are integers determined by the balancing C subroutine balanc. if balanc has not been used, C set low=1, igh=n. C C a contains the input matrix. C C on output C C a contains the hessenberg matrix. the multipliers C which were used in the reduction are stored in the C remaining triangle under the hessenberg matrix. C C int contains information on the rows and columns C interchanged in the reduction. C only elements low through igh are used. C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated august 1983. C C ------------------------------------------------------------------ C LA = IGH - 1 KP1 = LOW + 1 IF (LA .LT. KP1) GO TO 200 C DO 180 M = KP1, LA MM1 = M - 1 X = 0.0D0 I = M C DO 100 J = M, IGH IF (DABS(A(J,MM1)) .LE. DABS(X)) GO TO 100 X = A(J,MM1) I = J 100 CONTINUE C INT(M) = I IF (I .EQ. M) GO TO 130 C .......... interchange rows and columns of a .......... DO 110 J = MM1, N Y = A(I,J) A(I,J) = A(M,J) A(M,J) = Y 110 CONTINUE C DO 120 J = 1, IGH Y = A(J,I) A(J,I) = A(J,M) A(J,M) = Y 120 CONTINUE C .......... end interchange .......... 130 IF (X .EQ. 0.0D0) GO TO 180 MP1 = M + 1 C DO 160 I = MP1, IGH Y = A(I,MM1) IF (Y .EQ. 0.0D0) GO TO 160 Y = Y / X A(I,MM1) = Y C DO 140 J = M, N 140 A(I,J) = A(I,J) - Y * A(M,J) C DO 150 J = 1, IGH 150 A(J,M) = A(J,M) + Y * A(J,I) C 160 CONTINUE C 180 CONTINUE C 200 RETURN END C Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C C Careful! Anything free comes with no guarantee. C *** from netlib, Tue Mar 19 12:19:41 EST 1991 *** C SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) C RESTORED CORRECT INDICES OF LOOPS (200,210,230,240). (9/29/89 BSG) C INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR DOUBLE PRECISION H(NM,N),WR(N),WI(N) DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2 LOGICAL NOTLAS C C this subroutine is a translation of the algol procedure hqr, C num. math. 14, 219-231(1970) by martin, peters, and wilkinson. C handbook for auto. comp., vol.ii-linear algebra, 359-371(1971). C C this subroutine finds the eigenvalues of a real C upper hessenberg matrix by the qr method. C C on input C C nm must be set to the row dimension of two-dimensional C array parameters as declared in the calling program C dimension statement. C C n is the order of the matrix. C C low and igh are integers determined by the balancing C subroutine balanc. if balanc has not been used, C set low=1, igh=n. C C h contains the upper hessenberg matrix. information about C the transformations used in the reduction to hessenberg C form by elmhes or orthes, if performed, is stored C in the remaining triangle under the hessenberg matrix. C C on output C C h has been destroyed. therefore, it must be saved C before calling hqr if subsequent calculation and C back transformation of eigenvectors is to be performed. C C wr and wi contain the real and imaginary parts, C respectively, of the eigenvalues. the eigenvalues C are unordered except that complex conjugate pairs C of values appear consecutively with the eigenvalue C having the positive imaginary part first. if an C error exit is made, the eigenvalues should be correct C for indices ierr+1,...,n. C C ierr is set to C zero for normal return, C j if the limit of 30*n iterations is exhausted C while the j-th eigenvalue is being sought. C C questions and comments should be directed to burton s. garbow, C mathematics and computer science div, argonne national laboratory C C this version dated september 1989. C C ------------------------------------------------------------------ C IERR = 0 NORM = 0.0D0 K = 1 C .......... store roots isolated by balanc C and compute matrix norm .......... DO 50 I = 1, N C DO 40 J = K, N 40 NORM = NORM + DABS(H(I,J)) C K = I IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 WR(I) = H(I,I) WI(I) = 0.0D0 50 CONTINUE C EN = IGH T = 0.0D0 ITN = 30*N C .......... search for next eigenvalues .......... 60 IF (EN .LT. LOW) GO TO 1001 ITS = 0 NA = EN - 1 ENM2 = NA - 1 C .......... look for single small sub-diagonal element C for l=en step -1 until low do -- .......... 70 DO 80 LL = LOW, EN L = EN + LOW - LL IF (L .EQ. LOW) GO TO 100 S = DABS(H(L-1,L-1)) + DABS(H(L,L)) IF (S .EQ. 0.0D0) S = NORM TST1 = S TST2 = TST1 + DABS(H(L,L-1)) IF (TST2 .EQ. TST1) GO TO 100 80 CONTINUE C .......... form shift .......... 100 X = H(EN,EN) IF (L .EQ. EN) GO TO 270 Y = H(NA,NA) W = H(EN,NA) * H(NA,EN) IF (L .EQ. NA) GO TO 280 IF (ITN .EQ. 0) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 C .......... form exceptional shift .......... T = T + X C DO 120 I = LOW, EN 120 H(I,I) = H(I,I) - X C S = DABS(H(EN,NA)) + DABS(H(NA,ENM2)) X = 0.75D0 * S Y = X W = -0.4375D0 * S * S 130 ITS = ITS + 1 ITN = ITN - 1 C .......... look for two consecutive small C sub-diagonal elements. C for m=en-2 step -1 until l do -- .......... DO 140 MM = L, ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R * S - W) / H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = DABS(P) + DABS(Q) + DABS(R) P = P / S Q = Q / S R = R / S IF (M .EQ. L) GO TO 150 TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1))) TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R)) IF (TST2 .EQ. TST1) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.0D0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0D0 160 CONTINUE C .......... double qr step involving rows l to en and C columns m to en .......... DO 260 K = M, NA NOTLAS = K .NE. NA IF (K .EQ. M) GO TO 170 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0D0 IF (NOTLAS) R = H(K+2,K-1) X = DABS(P) + DABS(Q) + DABS(R) IF (X .EQ. 0.0D0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P) IF (K .EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 190 P = P + S X = P / S Y = Q / S ZZ = R / S Q = Q / P R = R / P IF (NOTLAS) GO TO 225 C .......... row modification .......... DO 200 J = K, EN P = H(K,J) + Q * H(K+1,J) H(K,J) = H(K,J) - P * X H(K+1,J) = H(K+1,J) - P * Y 200 CONTINUE C J = MIN0(EN,K+3) C .......... column modification .......... DO 210 I = L, J P = X * H(I,K) + Y * H(I,K+1) H(I,K) = H(I,K) - P H(I,K+1) = H(I,K+1) - P * Q 210 CONTINUE GO TO 255 225 CONTINUE C .......... row modification .......... DO 230 J = K, EN P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J) H(K,J) = H(K,J) - P * X H(K+1,J) = H(K+1,J) - P * Y H(K+2,J) = H(K+2,J) - P * ZZ 230 CONTINUE C J = MIN0(EN,K+3) C .......... column modification .......... DO 240 I = L, J P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2) H(I,K) = H(I,K) - P H(I,K+1) = H(I,K+1) - P * Q H(I,K+2) = H(I,K+2) - P * R 240 CONTINUE 255 CONTINUE C 260 CONTINUE C GO TO 70 C .......... one root found .......... 270 WR(EN) = X + T WI(EN) = 0.0D0 EN = NA GO TO 60 C .......... two roots found .......... 280 P = (Y - X) / 2.0D0 Q = P * P + W ZZ = DSQRT(DABS(Q)) X = X + T IF (Q .LT. 0.0D0) GO TO 320 C .......... real pair .......... ZZ = P + DSIGN(ZZ,P) WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ WI(NA) = 0.0D0 WI(EN) = 0.0D0 GO TO 330 C .......... complex pair .......... 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 C .......... set error -- all eigenvalues have not C converged after 30*n iterations .......... 1000 IERR = EN 1001 RETURN END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: DFLRD.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C C DRIVER PROGRAM FOR DSEVAS C C This program prompts the user for a file name that C contains the input, reads the file, allocates the desired C eigenvalues, and finds the eigenvalues of the C matrix (A - b*f). The input values, output vector f, C and the resulting eigenvalues are printed. C C INPUT FORMAT: C ------------- C C 1) Line 1 of the file contains an integer in free format C showing the size of the system. C C 2) Line 2 of the file contains a logical value in free format C showing the form of the system. A value of "T" C indicates a system in controllability form, whereas C a value of "F" indicates a system in general form. C C 3) Line 3 of the file contains a real value in free format C giving the machine unit used for testing for C controllability. If the value 0.0 is given C the subroutine calculates the machine unit. C C 4) The matrix A is read next, one element per line, C in free format. The matrix is read by row. C C 5) The vector b is read next, one element per line, C in free format. C C 6) The desired eigenvalues are read next. Each line C contains an ordered pair of two real, free format C values which represent the real and the C imaginary part of a complex eigenvalue. C Real eigenvalues should obviously be represented C by an ordered pair with a zero imaginary part. C Complex conjugates should appear both, and C contiguously. Real eigenvalues should appear in C groups of even numbers, unless they are at the end C of the list. A good way to list the eigenvalues C would be, to first list all the complex eigenvalues C followed be the real eigenvalues. C C A typical example could be the following: C C 4 C T C 0.0 C 1 For A= 1 2 3 4 C 2 5 6 7 8 C 3 0 9 1 2 C 4 0 0 3 4 C 5 C 6 C 7 C 8 C 0 C 9 C 1 C 2 C 0 C 0 C 3 C 4 C 10 For b= 10 C 0 0 C 0 0 C 0 0 C 1 2 C 1 -2 C 3 0 C 4 0 C C Note that in this test driver program we have set a limit C of 25 for the size of the matrix A. This obviously C can be altered by the user. C C********************************************************************* C********************************************************************* PROGRAM SAMPLE PARAMETER (ISL = 25, JSL = ISL*ISL+2*ISL) DOUBLE PRECISION A(ISL,ISL),AA(ISL,ISL),B(ISL),BB(ISL) DOUBLE PRECISION F(ISL),RST(4,JSL),CST(JSL) INTEGER IST(JSL),N,IEAL COMPLEX*16 EIG(ISL) DOUBLE PRECISION L(2,ISL), RMU LOGICAL FORM CHARACTER*20 FNAME C C.......read matrix A, vector b, and eigenvalue matrix L C WRITE(*,*) WRITE(*,*) '*****************************' WRITE(*,*) ' DSEVAS TEST PROGRAM' WRITE(*,*) '*****************************' WRITE(*,*) WRITE(*,4000) 4000 FORMAT(1X,'DATA FILE: ',/) READ(*,7000) FNAME 7000 FORMAT(A) OPEN(5, FILE=FNAME, STATUS='OLD') READ(5,*) N IF(N .GT. ISL) THEN WRITE(*,*) 'INPUT SYSTEM TOO LARGE.' STOP ENDIF READ(5,*) FORM READ(5,*) RMU DO 100 I=1,N DO 100 J=1,N READ(5,*) A(I,J) AA(I,J) = A(I,J) 100 CONTINUE DO 110 I=1,N READ(5,*) B(I) BB(I) = B(I) 110 CONTINUE DO 120 I=1,N READ(5,*) L(1,I), L(2,I) 120 CONTINUE CLOSE(5) C C.......Echo input C WRITE(*,*) WRITE(*,*) WRITE(*,*) 'SYSTEM MATRIX (b,A):' WRITE(*,*) DO 130 I=1,N WRITE(*,5000) B(I),(A(I,J),J=1,N) 130 CONTINUE WRITE(*,*) WRITE(*,*) 'EIGENVALUES TO BE ALLOCATED:' WRITE(*,*) DO 135 I=1,N WRITE(*,5000) L(1,I),L(2,I) 135 CONTINUE C C.......call DSEVAS routine C CALL DSEVAS(A,N,ISL,B,L,RST,IST,CST,FORM,RMU,IEAL,F) C C.......create (A - b*f) C WRITE(*,*) WRITE(*,*) 'RESULT VECTOR F:' WRITE(*,*) DO 140 I=1,N WRITE(*,*) F(I) AA(1,I) = AA(1,I) - BB(1)*F(I) 140 CONTINUE C C.......find eigenvalues of (A - b*f) C CALL DEIGVL(N,AA,ISL,EIG) WRITE(*,*) WRITE(*,*) 'EIGENVALUES OF (A - B*F): ' WRITE(*,*) DO 150 I=1,N WRITE(*,*) EIG(I) 150 CONTINUE 5000 FORMAT(1X,10E15.6) STOP END CXXXXXXXXXXXXXXXXXXXXX START OF FILE: DSEVAS.F XXXXXXXXXXXXXXXXXXXXX C********************************************************************* C********************************************************************* C********************************************************************* C C Title: DSEVAS C C Authors: G. Miminis and M. Reid C C Purpose: Calculates f such that (A - b*f) has C eigenvalues eigvl. A and b may be C general, or in "controllability" form, C where A is unreduced upper hessenberg C and b has one non-zero element in b(1). C The eigenvalues stored in eigvl may be C real or complex conjugate pairs. C C C Algorithm: C ---------- C C This subroutine implements an algorithm by G. Miminis and C C. Paige, submitted to the SIAM Journal on Matrix Analysis C and Applications. For details please see the paper. C C C Input Arguments: C ---------------- C C [ (+) means argument is altered by the subroutine. ] C C (+) A real*8(lda,n) Matrix A C C n integer Order of matrix A C C lda integer First dimension of matrix A as C defined in the calling program C C (+) binp real*8(n) Input vector b. C C eigvl real*8(2,n) Desired eigenvalues. EIGVL(1,*) C holds the real part; EIGVL(2,*) C holds the imaginary. C C (+) rstor real*8(4,*) Working storage for reflectors. C Must be dimensioned as C RSTOR(4,k), where: C C k >= (n**2 - 2*n + 1)/4 C C C (+) irpos integer(*) Position of reflectors. Must C be dimensioned as IRPOS(k), C where k is defined as for C RSTOR. C C (+) cstor real*8(*) Working storage for C conversion to controllability C form. Must be dimensioned as C CSTOR(w), where: C C w >= (n**2 + 3*n - 4)/2 C C C iform logical Input form. C .TRUE. = controllability C form C .FALSE. = general C C C (+) rmu real*8 The machine unit. C If positive, the subroutine C will use the given value C as the machine unit. C If the given value of rmu is C zero or negative the subroutine C will automatically calculate the C machine unit. (See note 3 C below.) The calculated C machine unit is returned via C this parameter. C C The tolerance used in testing C for controllability is: C C tol = ||(b,A)|| * u , C C where u is the machine unit C and ||.|| is the 1-norm. C C C Output Arguments: C -------------- C C (+) f real*8(n) Solution vector C C (+) ieigal integer Number of eigenvalues actually C allocated. C C Local Variables: C ---------------- C C isize integer Current size of matrix A. C C le integer Position of leading edge of C matrix A in array A. C C irefc integer Index of last elementary C reflector stored in rstor. C C eps real*8 Tolerance used in testing C for controllability. C C C C Subroutines Used: C ----------------- C C APPSH, CRH, PREM, GETEPS, EIGCHK, CONTMK, POST, DSTORE C C C Notes on Usage: C --------------- C C 1) DSEVAS only takes eigenvalues that are real or in C complex conjugate pairs. These are passed to the C subroutine via EIGVL, dimensioned as EIGVL(2,N). C EIGVL(1,x) holds the real part, while EIGVL(2,x) holds C the imaginary part. Real valued eigenvalues can be in C any order. Conjugate pairs of eigenvalues must be C stored contiguously, with the first of the pair on an C odd-numbered index in the array EIGVL. For example, C the eigenvalues 1,2,3, and 1 +/- 4 can be stored: C C [1.0,0.0] C [2.0,0.0] C [1.0,+4.0] C [1.0,-4.0] C [3.0,0.0] C C but the following would generate an error: C C [1.0,0.0] C [1.0,+4.0] C [1.0,-4.0] C [2.0,0.0] C [3.0,0.0] C C since the subroutine working with two eigenvalues C at a time, will work with eigenvalues 1.0 and C 1.0 + i*4.0 , which is not allowed. An easy way C to list the eigenvalues could be to first list C all the complex eigenvalues, followed by the real C eigenvalues. C C 2) Before allocating a pair of eigenvalues, the subroutine C checks to see if the current matrix is controllable. C If this is not the case, no further allocations will C take place, and the subroutine will print a warning C indicating the number of successful allocations. C The resulting matrix (A-b*f) will have those C eigenvalues allocated; the remainder will be undefined. C C 3) The DSEVAS routine uses the machine unit as part of the C tolerance for controllability determination. If C the input parameter RMU is zero or negative, the C internal routine GETEPS automatically calculates C the machine unit. The user should note, however, C that such a calculation can fail and give the C double precision machine unit, even in the single C precision version. (See Gentleman and Marovich C [CACM 17, p. 276, 1974]). For this reason, the C user has the option of supplying the machine unit, C possibly from the PORTS machine constants sub- C routines. (See TOMS 4 (1978), p. 177). The C parameter RMU returns the machine unit actually C used in the calculations. C C 4) This subroutine makes extensive use of the level-1 BLAS C (Basic Linear Algebra Subroutines) package. The C BLAS-1 subroutines used are: C C DASUM, DAXPY, DCOPY, DDOT, DNRM2, DSCAL, IDAMAX C C These routines are available from various sources, C such as the IMSL subroutine library. C C 5) A named common block (MACEPS) is used in this subroutine. C C********************************************************************* C********************************************************************* C********************************************************************* SUBROUTINE DSEVAS(A,N,LDA,BINP,EIGVL,RSTOR,IRPOS, + CSTOR,IFORM,RMU,IEIGAL,F) C C........input arguments C INTEGER N,LDA,IRPOS(*),IEIGAL DOUBLE PRECISION A(LDA,*),BINP(N),F(N),RSTOR(4,*),CSTOR(*) DOUBLE PRECISION EIGVL(2,N),RMU LOGICAL IFORM C C.......local variables C INTEGER ISIZE,LE,IREFC,IOLDRF DOUBLE PRECISION S1,S2,LM,EPS,GETEPS LOGICAL EIGCHK C C.......convert general form to controllability form C IF(.NOT. IFORM) CALL CONTMK(A,N,LDA,BINP,CSTOR) C C.......calculate epsilon for controllability determination C EPS = GETEPS(A,N,LDA,BINP,RMU) C C.......clear zero elements of system C DO 50 J=1,N-2 F(J) = 0.0 DO 50 I=J+2,N A(I,J) = 0.0 50 CONTINUE F(N-1) = 0.0 F(N) = 0.0 C C.......start deflation: continue by double-step until problem is C smaller than size 3 C IEIGAL = 0 IREFC = 0 ISIZE = N 100 IF ( ISIZE .LE. 2) GOTO 200 LE = N - ISIZE + 1 IOLDRF = IREFC C C.......check eigenvalues for complex conjugates C IF(.NOT. EIGCHK(EIGVL,LE,EPS)) THEN WRITE(*,5010) WRITE(*,5005) N-ISIZE GOTO 900 ENDIF C C.......apply reflector containing eigenvalues to be allocated. C CALL APPSH(A,LE,N,LDA,EIGVL,RSTOR,IRPOS,IREFC) C C.......transform A to block upper hessenberg. C CALL CRH(A,LE,N,LDA,RSTOR,IRPOS,IREFC) C C.......apply latest reflector to binp (vector of inputs) C BINP(LE+1) = 0.0 BINP(LE+2) = 0.0 CALL PREM(BINP,LDA,LE,1,1,RSTOR(1,IREFC),3) C C......check system for controllability C IF(DABS(BINP(LE+2)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS WRITE(*,5005) LE-1 IREFC = IOLDRF F(LE) = 0.0 F(LE+1) = 0.0 GOTO 900 C C.......specify le,le+1 elements of f C ELSE F(LE) = A(LE+2,LE) / BINP(LE+2) F(LE+1) = A(LE+2,LE+1) / BINP(LE+2) ENDIF C C.......reduce problem size for next step C ISIZE = ISIZE - 2 IEIGAL = IEIGAL + 2 GOTO 100 200 CONTINUE C C.......allocate eigenvalues for matrices of size one or two. C IF (ISIZE .EQ. 1) THEN C C...........check for controllability C IF(DABS(BINP(N)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS WRITE(*,5005) N-1 F(N) = 0.0 GOTO 900 ELSE F(N) = (A(N,N) - EIGVL(1,N))/BINP(N) IEIGAL = IEIGAL + 1 ENDIF ELSE C C...........check for controllability and conjugate C eigenvalues (size 2) C IF(.NOT. EIGCHK(EIGVL,LE,EPS) .OR. + DABS(BINP(N-1)) .LT. EPS .OR. + DABS(A(N,N-1)) .LT. EPS) THEN IF(DABS(BINP(N-1)) .LT. EPS .OR. + DABS(A(N,N-1)) .LT. EPS) THEN WRITE(*,5000) WRITE(*,5002) EPS ELSE WRITE(*,5010) ENDIF WRITE(*,5005) N-2 F(N-1) = 0.0 F(N) = 0.0 GOTO 900 ELSE LM = EIGVL(1,N-1)*EIGVL(1,N) - EIGVL(2,N-1)*EIGVL(2,N) S1 = EIGVL(1,N-1)+EIGVL(1,N) - A(N,N) S2 = (S1*A(N,N) - LM)/A(N,N-1) F(N-1) = (A(N-1,N-1) - S1)/BINP(N-1) F(N) = (A(N-1,N) - S2)/BINP(N-1) IEIGAL = IEIGAL + 2 ENDIF ENDIF 900 CONTINUE C C.......apply accumulated reflectors to f vector C DO 300 K=IREFC,1,-1 CALL PREM(F,LDA,IRPOS(K),1,1,RSTOR(1,K),3) 300 CONTINUE C C.......reverse transformation to controllability form for C vector f C IF(.NOT. IFORM) THEN II=N*(N+1)/2 + N - 4 DO 350 K=2,N CALL POST(F,1,N-K+1,N-K+1,1,CSTOR(II),K) II = II - K - 2 350 CONTINUE ENDIF 5000 FORMAT(1X,'SYSTEM IS CLOSE TO UNCONTROLLABLE') 5002 FORMAT(1X,'USING TOLERANCE ',E20.12) 5005 FORMAT(1X,I4,' EIGENVALUES ALLOCATED.') 5010 FORMAT(1X,'COMPLEX EIGENVALUES ARE NOT COMPLEX CONJUGATES.') RETURN END C********************************************************************* C C Title: CONTMK C C Purpose: Converts the system (b,A) into C controllability form. C C C Input Arguments: C ---------------- C C A real*8(lda,n) Control matrix A C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C binp real*8(n) Input vector C C cstor real*8(*) Working storage for C conversion to controllability C form. C C C C Local Variables: C ---------------- C C ipnt integer Pointer to current reflector C in CSTOR() C C C Subroutines Used: C ----------------- C C DCOPY, DNRM2, VHOUSE, PREM, POST C C********************************************************************* SUBROUTINE CONTMK(A,N,LDA,BINP,CSTOR) INTEGER N,LDA,IPNT DOUBLE PRECISION A(LDA,*),BINP(N),CSTOR(*),DNRM2 IPNT = 1 DO 100 I=N,2,-1 IF(I .EQ. N) THEN CALL DCOPY(I, BINP(1), 1, CSTOR(IPNT+1), 1) BINP(1) = -DSIGN(DNRM2(I, BINP(1), 1),BINP(1)) DO 150 J=2,N BINP(J) = 0.0 150 CONTINUE ELSE CALL DCOPY(I, A(N-I+1,N-I), 1, CSTOR(IPNT+1), 1) A(N-I+1,N-I) = -DSIGN(DNRM2(I, A(N-I+1,N-I), 1), + A(N-I+1,N-I)) DO 175 J=N-I+2,N A(J,N-I) = 0.0 175 CONTINUE ENDIF CALL VHOUSE(CSTOR(IPNT), I) CALL PREM(A,LDA,N-I+1,N-I+1,N,CSTOR(IPNT),I) CALL POST(A,LDA,1,N,N-I+1,CSTOR(IPNT),I) IPNT = IPNT + I + 1 100 CONTINUE RETURN END C********************************************************************* C C Title: VHOUSE C C Purpose: Given a n-element vector stored in v(2)- C v(n+1), HOUSE calculates pi and u such that C I - u * uT / pi will introduce zeroes into C first two elements of v. Pi is returned in C v(1), while u overwrites the input vector C v(2)-v(n+1). C C Input Arguments: C ---------------- C C v real*8(*) The last n elements of v contain C the input vector to be zeroed. C C C Output Values: C -------------- C C v v(1) = pi C v(2)-v(n+1) = u C C C Subroutines used: C ----------------- C C DSCAL, DNRM2, IDAMAX C C********************************************************************* C C.......overwrites v with u of elem. reflector; returns p C SUBROUTINE VHOUSE(V,N) INTEGER N,I DOUBLE PRECISION V(*),SIGMA,ALFA,DNRM2 I = IDAMAX( N, V(2), 1) ALFA = 1.0/DABS(V(I+1)) CALL DSCAL( N, ALFA, V(2), 1) SIGMA = DSIGN(DNRM2(N, V(2), 1),V(2)) V(2) = V(2) + SIGMA V(1) = SIGMA*V(2) RETURN END C********************************************************************* C C Title: EIGCHK C C Purpose: Checks that complex eigenvalues are C complex conjugates. The function C returns .TRUE. if the eigenvalues are C real or complex conjugates, .FALSE. C otherwise. C C C C Input Arguments: C ---------------- C C eigvl real*8(2,n) Array of eigenvalues to C be allocated. C C le integer Position of leading edge of C matrix A in array A. C C eps real*8 Epsilon value for C comparison to zero. C C eps = ||(b,A)||u C C where ||.|| is the C 1-norm of the system, C and u is the machine unit. C C Local Variables: C ---------------- C C rd1 real*8 Relative difference in real C parts. C C rd2 real*8 Relative difference in complex C parts. C C norm,tnorm,eps C C C Subroutines Used: C ----------------- C C DASUM C C********************************************************************* LOGICAL FUNCTION EIGCHK(EIGVL,LE,EPS) DOUBLE PRECISION EIGVL(2,*),EPS,RD1,RD2 INTEGER LE EIGCHK = .TRUE. C C.......check if complex eigenvalues are complex conjugates C IF(DABS(EIGVL(2,LE)) .GT. EPS .OR. + DABS(EIGVL(2,LE+1)) .GT. EPS) THEN RD1 = DABS((EIGVL(1,LE) - EIGVL(1,LE+1))/EIGVL(1,LE+1)) RD2 = DABS((EIGVL(2,LE) + EIGVL(2,LE+1))/EIGVL(2,LE+1)) IF(RD1 .GT. EPS .OR. RD2 .GT. EPS) THEN EIGCHK = .FALSE. ENDIF ENDIF RETURN END C********************************************************************* C C Title: GETEPS C C Purpose: Calculates the tolerance used in C determining the uncontrollability of C a system. To compute the tolerance, the C routine computes the 1-norm of the system C (b,A) and multiplies this value by the C machine unit: C C eps = ||(b,A)||*u C C where ||.|| is the 1-norm, and u is the C machine unit. C C If the input parameter RMU is positive, C that value is used as the machine unit in C calculating EPS. If RMU is zero or negative, C the routine calculates the machine unit. C C C Input Arguments: C ---------------- C C A real*8(lda,n) Control matrix A C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C b real*8(n) Input vector b C C rmu real*8 User value for machine unit C C Local Variables: C ---------------- C C norm,tnorm,eps C C C Subroutines Used: C ----------------- C C DASUM, DSTORE C C********************************************************************* DOUBLE PRECISION FUNCTION GETEPS(A,N,LDA,B,RMU) INTEGER N,LDA DOUBLE PRECISION A(LDA,*),B(N),RMU DOUBLE PRECISION NORM,TNORM,EEPS,EPS,DASUM COMMON/MACEPS/EEPS C C.......use user value for machine unit if non-zero C IF(RMU .GT. 0.0) THEN EPS = RMU ELSE C.......calculate machine unit EPS = 1.0 30 EPS = EPS / 2.0 EEPS=EPS + 1.0 CALL DSTORE(EEPS) IF(EEPS .NE. 1.0) GOTO 30 EPS = EPS * 2.0 RMU = EPS ENDIF C C.......calculate 1 norm of system (b,A) C NORM = DASUM(N, B(1), 1) DO 50 I=1,N TNORM = DASUM(N, A(1,I), 1) IF(TNORM .GT. NORM) NORM = TNORM 50 CONTINUE GETEPS = EPS * NORM RETURN END C C.......This subroutine ensures that the value E is actually C stored in memory, with appropriate rounding errors. C (From Gentleman and Marovich, Communications of the ACM, C vol. 17, no. 5, May 1974. p. 277) C SUBROUTINE DSTORE(E) DOUBLE PRECISION E,V COMMON/MACEPS/V V = E RETURN END C********************************************************************* C C Title: CRH C C Purpose: Uses a series of elementary reflectors C to restore the matrix A to a "near" upper C Hessenberg form. C C C Input Arguments: C ---------------- C C A real*8(lda,n) Control matrix A C C le integer Position of leading edge of C matrix A in array A. C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C rstor real*8(4,*) Working storage for C reflectors C C irpos integer(*) Position of reflectors C C irefc integer Index of last elementary C reflector stored in rstor. C C C Output Values: C -------------- C C A A is returned, transformed into near C upper Hessenberg. C C rstor,irpos,irefc Contain reflectors used to transform C A to near upper Hessenberg. C C Local Variables: C ---------------- C C C C Subroutines Used: C ----------------- C C DCOPY, HOUSE, DNRM2, POST, PREM C C********************************************************************* SUBROUTINE CRH(A,LE,N,LDA,RSTOR,IRPOS,IREFC) INTEGER LE,N,LDA,I,IREFC,IRPOS(*) DOUBLE PRECISION A(LDA,*),P,RSTOR(4,*),DNRM2 DO 100 I=N,LE+3,-1 IREFC = IREFC + 1 IRPOS(IREFC) = I-3 CALL DCOPY(3, A(I,I-3), LDA, RSTOR(2,IREFC), 1) CALL HOUSE(RSTOR(1,IREFC)) P = -DSIGN(DNRM2(3, A(I,I-3), LDA),A(I,I-1)) CALL POST(A,LDA,LE,I-1,I-3,RSTOR(1,IREFC),3) A(I,I-1) = P A(I,I-2) = 0.0 A(I,I-3) = 0.0 CALL PREM(A,LDA,I-3,LE,N,RSTOR(1,IREFC),3) 100 CONTINUE RETURN END C********************************************************************* C C Title: APPSH C C Purpose: Calculates an elementary reflector C containing the implicit shift and pre- C and post multiplies A. C C Input Arguments: C ---------------- C C A real*8(lda,n) Control matrix A. C C le integer Position of leading edge of C matrix A in array A. C C n integer Order of matrix A C C lda integer First dimension of array A as C defined in the calling program C C eigvl real*8(2,n) Vector of eigenvalues to be C allocated. C C rstor real*8(4,*) Working storage for reflectors. C C irpos integer(*) Position of reflectors. C C irefc integer Index of last reflector in C irpos and rstor. C C C Output Values: C -------------- C C A Returns the original matrix A pre- and C post-multiplied by the elementary reflector. C C rstor,irpos Stores the value and position of the C elementary reflector. C C la,lm C C Subroutines used: C ----------------- C C HOUSE, PREM, POST C C********************************************************************* SUBROUTINE APPSH(A,LE,N,LDA,EIGVL,RSTOR,IRPOS,IREFC) INTEGER LE,N,LDA,IRPOS(*),IREFC DOUBLE PRECISION A(LDA,*),RSTOR(4,*) DOUBLE PRECISION EIGVL(2,*) DOUBLE PRECISION LA,LM IREFC = IREFC + 1 IRPOS(IREFC) = N-2 C C.......construct bottom row of A~ C LA = EIGVL(1,LE) + EIGVL(1,LE+1) LM = EIGVL(1,LE)*EIGVL(1,LE+1) - EIGVL(2,LE)*EIGVL(2,LE+1) RSTOR(2,IREFC) = A(N,N-1)*A(N-1,N-2) RSTOR(3,IREFC) = A(N,N-1)*A(N-1,N-1) + A(N,N)*A(N,N-1) - + LA*A(N,N-1) RSTOR(4,IREFC) = A(N,N-1)*A(N-1,N) + A(N,N)*A(N,N) - + LA*A(N,N) + LM C C.......convert A to PtAP C CALL HOUSE(RSTOR(1,IREFC)) CALL POST(A,LDA,LE,N,N-2,RSTOR(1,IREFC),3) CALL PREM(A,LDA,N-2,LE,N,RSTOR(1,IREFC),3) RETURN END C********************************************************************* C C Title: PREM C C Purpose: Premultiplies the matrix contained in array A C by the 3-element elementary reflector stored C in vector v. C C Input Arguments: C ---------------- C C A real*8(lda,nv) Control matrix A. C C lda integer First dimension of array A as C defined in the calling program C C irow integer First row of A affected by C the reflector. (Reflector will C be applied to rows *irow* to C *irow+nv-1*.) C C j,k integer Range of columns of A on C which to apply reflector. C C v real*8(4) Contains the elementary C reflector. v(1) = pi and C v(2)->v(nv+1) = u, where C I - u*uT/pi is the reflector. C C nv integer Dimension of vector v. C C C Output Values: C -------------- C C A Returns the original matrix A pre- C multiplied by the elementary reflector. C C C C Subroutines used: C ----------------- C C DDOT, DAXPY C C********************************************************************* SUBROUTINE PREM(A,LDA,IROW,J,K,V,NV) INTEGER IROW,J,K,LDA,I,NV DOUBLE PRECISION A(LDA,*),V(*),ALPHA,DDOT DO 100 I=J,K ALPHA = DDOT(NV, A(IROW,I), 1, V(2), 1)/V(1) CALL DAXPY(NV, -ALPHA, V(2), 1, A(IROW,I), 1) 100 CONTINUE RETURN END C********************************************************************* C C Title: POST C C Purpose: post-multiplies the matrix contained in C array A by the 3-element elementary reflector C stored in vector v. C C Input Arguments: C ---------------- C C A real*8(lda,n) Control matrix A. C C lda integer First dimension of array A as C defined in the calling program C C j,k integer Range of rows of A on C which to apply reflector. C C icol integer First column of A affected by C the reflector. (Reflector will C be applied to columns *icol* C to *icol+nv-1*.) C C v real*8(4) Contains the elementary C reflector. v(1) = pi and C v(2)->v(nv+1) = u, where C I - u*uT/pi is the reflector. C C nv integer Dimension of vector v. C C C Output Values: C -------------- C C A Returns the original matrix A C post-multiplied by the elementary reflector. C C C C Subroutines used: C ----------------- C C DDOT, DAXPY C C********************************************************************* SUBROUTINE POST(A,LDA,J,K,ICOL,V,NV) INTEGER J,K,ICOL,LDA,I,NV DOUBLE PRECISION A(LDA,*),V(*),ALPHA,DDOT DO 100 I=J,K ALPHA = DDOT(NV, A(I,ICOL), LDA, V(2), 1)/V(1) CALL DAXPY(NV, -ALPHA, V(2), 1, A(I,ICOL), LDA) 100 CONTINUE RETURN END C********************************************************************* C C Title: HOUSE C C Purpose: Given a 3-element vector stored in v(2)- C v(4), HOUSE calculates pi and u such that C I - u * uT / pi will introduce zeroes into C first two elements of v. Pi is returned in C v(1), while u overwrites the input vector C v(2)-v(4). C C Input Arguments: C ---------------- C C v real*8(4) The last three elements of v contain C the input vector to be zeroed. C C C Output Values: C -------------- C C v v(1) = pi C v(2)-v(4) = u C C C Subroutines used: C ----------------- C C DSCAL, DNRM2, IDAMAX C C********************************************************************* C C.......overwrites v with u of elem. reflector; returns p C SUBROUTINE HOUSE(V) INTEGER I DOUBLE PRECISION V(4),SIGMA,ALFA,DNRM2 I = IDAMAX( 3, V(2), 1) ALFA = 1.0/DABS(V(I+1)) CALL DSCAL( 3, ALFA, V(2), 1) SIGMA = DSIGN(DNRM2(3, V(2), 1),V(4)) V(4) = V(4) + SIGMA V(1) = SIGMA*V(4) RETURN END