C ALGORITHM 696 , COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 17, NO. 3, SEPTEMBER, 1991, PP. 335-340.
C***********************************************************************
C
C TESTDRIVER FOR THE SUBROUTINE GIRI JULY 16, 1990
C
C THIS PROGRAM IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN IMPLE-
C MENTATION OF A TESTDRIVER FOR THE SUBROUTINE GIRI.
C
C THERE ARE FOUR TESTS BASED ON AN ASYMMETRIC, COMPLEX BAND MATRIX.
C THE REAL PART AR OF THE MATRIX IS:
C
C
C 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. 0. 0. ...
C
C -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. 0. ...
C
C 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. 0. ...
C
C -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. 0. ...
C
C 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. 9.-10. ...
C
C -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. -8. 9. ...
C
C 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. 7. -8. ...
C
C -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. -6. 7. ...
C
C 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. 5. -6. ...
C
C -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. -4. 5. ...
C
C 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. 3. -4. ...
C
C 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. -2. 3. ...
C
C 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. 1. -2. ...
C
C 0. 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 0. 1. ...
C
C 0. 0. 0. 0. 10. -9. 8. -7. 6. -5. 4. -3. 2. -1. 20. ...
C . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . .
C
C
C THE IMAGINARY PART AI OF THE MATRIX IS:
C
C
C 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. 10. 0. 0. 0. 0. ...
C
C 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. 10. 0. 0. 0. ...
C
C -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. 10. 0. 0. ...
C
C 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. 10. 0. ...
C
C -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. 10. ...
C
C 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. -9. ...
C
C -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. 8. ...
C
C 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. -7. ...
C
C -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. 6. ...
C
C 9. -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. -5. ...
C
C -10. 9. -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. 4. ...
C
C 0.-10. 9. -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. -3. ...
C
C 0. 0.-10. 9. -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. 2. ...
C
C 0. 0. 0.-10. -9. -8. 7. -6. 5. -4. 3. -2. 1. 0. -1. ...
C
C 0. 0. 0. 0.-10. 9. -8. 7. -6. 5. -4. 3. -2. 1. 0. ...
C . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . .
C . . . . . . . . . . . . . . .
C
C THE REAL AND THE IMAGINARY PARTS OF THE MATRIX ARE CALCULATED IN
C THE REQUIRED FORM BY THE SUBROUTINE BNDMAT WHICH IS SUPPLIED
C TOGETHER WITH THIS TESTDRIVER.
C
C THE SPECTRUM OF THE MATRIX IS:
C
C . . .
C . . .
C . . . XX
C . . . XXX
C 90 ..................................................XXX....
C . . XXX
C . . X.
C . . XXX .
C . . XX .
C . . X* .
C . . XX .
C . . XX .
C . . XXX .
C . . XXXXXXX .
C 0 ............................XXXXXX.......................
C . XXXXXXX .
C . XXX . .
C . XX . .
C . XX . .
C . XX . .
C . XX . .
C . XXX . .
C . X . .
C . XXX . .
C -90 .........XXX.............................................
C XXX . .
C XX. . .
C . . .
C . . .
C -90 0 90
C
C NOTE THAT AN "X" IN THE ABOVE PLOT SIMULATION CAN REPRESENT MORE
C THAN ONE EIGENVALUE.
C
C THE TEST PROGRAM CALCULATES THE COMPLEX EIGENVALUE
C
C (OMEGA_R,OMEGA_I) = (60.55191,48.46971)
C
C MARKED WITH AN ASTERISK IN THE ABOVE PLOT SIMULATION IN THREE
C DIFFERENT WAYS: IN THE FIRST TEST, THE EIGENVALUE IS COMPUTED
C USING THE INTERNAL INITIAL VECTORS AND AN ITERATION SCHEME WITH
C TWO SIMPLIFIED ITERATION STEPS AFTER THE FIRST FULL ONE. THIS
C SHOULD BE THE STANDARD PROCEDURE, BECAUSE THE TWO SIMPLIFIED STEPS
C IMPROVE THE INITIAL VECTORS CONSIDERABLY. IN THE SECOND TEST, THE
C RIGHT AND LEFT EIGENVECTORS COMPUTED IN THE FIRST TEST ARE TAKEN
C AS INITIAL VECTORS. ONE FULL ITERATION CYCLE IS SAVED. THIS DEMON-
C STRATES THAT THE USER SHOULD APPLY BETTER INITIAL VECTORS IF THEY
C ARE AVAILABLE. IN THE THIRD TEST, THE SIMPLIFIED STEPS ARE OMITTED
C SIX FULL ITERATION STEPS ARE NEEDED INSTEAD OF THREE IN THE FIRST
C TEST. THIS SHOWS THE EFFECT OF THE TWO SIMPLIFIED STEPS AFTER THE
C FIRST FULL STEP. THE FOURTH TEST CHECKS THE ERROR CONDITION
C "IERROR=1."
C
C THE RAYLEIGH ITERATION CONVERGES SAFELEY IF THE INITIAL VALUE IS
C CLOSE ENOUGH TO AN ISOLATED EIGENVALUE. HOWEVER, IN THE NEIGHBOR-
C HOOD OF A PAIR OF EIGENVALUES, A SMALL CHANGE IN THE INITIAL VALUE
C CAN CAUSE THE ALGORITHM TO CONVERGE TO DIFFERENT EIGENVALUES. THIS
C CAN BE SEEN, FOR EXAMPLE, WITH THE EIGENVALUES
C
C (OMERAR,OMEGAI) = (106.6606,95.5099),
C (OMERAR,OMEGAI) = (106.3353,95.1810),
C
C WHICH CAN BE OBTAINED BY USING EITHER (106.5,95.4) OR (106.4,95.4)
C AS INITIAL VALUES.
C
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
PARAMETER ( NP = 500 )
C
DOUBLE PRECISION AR(21,NP),AI(21,NP),AWR(21,NP),AWI(21,NP),
. XLR(10,NP),XLI(10,NP),
. UR(NP),UI(NP),VR(NP),VI(NP),WR(NP),WI(NP),
. EPS,OMEGAR,OMEGAI
INTEGER JP(NP),IERROR,INIT,IUNIT,M,N
C
IUNIT = 6
C
M = 21
N = NP
C
EPS = 1.D-8
C
WRITE (6,9000) N,M
C
CALL BNDMAT(AR,AI,N)
C
C-----------------------------------------------------------------------
C F I R S T T E S T: EIGENVALUE COMPUTATION USING
C INTERNAL INITIAL VECTORS
C-----------------------------------------------------------------------
C
WRITE (6,9001)
C
OMEGAR = 60.3
OMEGAI = 48.3
C
WRITE (6,9010) OMEGAR,OMEGAI
C
INIT = 0
C
CALL GIRI(AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
. INIT,WR,WI,OMEGAR,OMEGAI,EPS,2,1,
. IERROR,1,IUNIT)
C
WRITE (6,9020) OMEGAR,OMEGAI
C
WRITE (6,9030) IERROR
C
C-----------------------------------------------------------------------
C S E C O N D T E S T: EIGENVALUE COMPUTATION USING THE CALCU-
C LATED EIGENVECTORS AS INITIAL VECTORS
C-----------------------------------------------------------------------
C
WRITE (6,9002)
C
OMEGAR = 60.3
OMEGAI = 48.3
C
WRITE (6,9010) OMEGAR,OMEGAI
C
INIT = 1
C
CALL GIRI(AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
. INIT,WR,WI,OMEGAR,OMEGAI,EPS,2,1,
. IERROR,1,IUNIT)
C
WRITE (6,9020) OMEGAR,OMEGAI
C
WRITE (6,9030) IERROR
C
C-----------------------------------------------------------------------
C T H I R D T E S T: EIGENVALUE COMPUTATION WITHOUT
C THE SIMPLIFIED ITERATION STEPS
C-----------------------------------------------------------------------
C
WRITE (6,9003)
C
OMEGAR = 60.3
OMEGAI = 48.3
C
WRITE (6,9010) OMEGAR,OMEGAI
C
INIT = 0
C
CALL GIRI(AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
. INIT,WR,WI,OMEGAR,OMEGAI,EPS,0,1,
. IERROR,1,IUNIT)
C
WRITE (6,9020) OMEGAR,OMEGAI
C
WRITE (6,9030) IERROR
C
C-----------------------------------------------------------------------
C F O U R T H T E S T: M IS NOT ODD. ERROR CODE "IERROR=1"
C-----------------------------------------------------------------------
C
WRITE (6,9004)
C
M = 20
C
CALL GIRI(AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
. INIT,WR,WI,OMEGAR,OMEGAI,EPS,0,1,
. IERROR,1,IUNIT)
C
WRITE (6,9030) IERROR
C
STOP
C
9000 FORMAT(//' ',77('-')/
. ' TEST OF THE GENERALIZED INVERSE RAYLEIGH ITERATION',
. ' USING A COMPLEX BAND'/
. 19(' '),'MATRIX OF DIMENSION',I4,' AND BAND WIDTH',I3/
. ' ',77('-')////)
9001 FORMAT(' ',77('-')/
. ' F I R S T T E S T: EIGENVALUE COMPUTATION USING',
. ' INTERNAL INITIAL VECTORS'/
. ' ',77('-')///)
9002 FORMAT(' ',77('-')/
. ' S E C O N D T E S T: EIGENVALUE COMPUTATION',
. ' USING THE CALCULATED'/
. 30(' '),'EIGENVECTORS AS INITIAL VECTORS'/
. ' ',77('-')///)
9003 FORMAT(' ',77('-')/
. ' T H I R D T E S T: EIGENVALUE COMPUTATION',
. ' WITHOUT THE SIMPLIFIED STEPS'/
. ' ',77('-')///)
9004 FORMAT(' ',77('-')/
. ' F O U R T H T E S T: ERROR CODE "IERROR=1"'/
. ' ',77('-')///)
9010 FORMAT(' INITIAL GUESS FOR EIGENVALUE: (OMEGA_R,OMEGA_I) = (',
. E10.3,',',E10.3,')'/)
9020 FORMAT(///' COMPUTED EIGENVALUE: (OMEGA_R,OMEGA_I) = (',
. E14.7,',',E14.7,')'///)
9030 FORMAT(' ERROR CODE: IERROR =',I2////)
END
C
C
C
C
C
C***********************************************************************
C
C NAME: BNDMAT ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C BNDMAT COMPUTES THE ELEMENTS OF THE BAND MATRIX FOR THE TEST OF
C THE SUBROUTINE GIRI.
C
C CALLING SEQUENCE: CALL BNDMAT(AR,AI,N)
C
C PARAMETERS:
C
C AR,AI: DOUBLE PRECISION, OUTPUT. AR(21,N), AI(21,N) CONTAIN
C REAL AND IMAGINARY PARTS OF THE BAND MATRIX.
C N: INTEGER, INPUT, DIMENSION OF THE MATRIX.
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE BNDMAT(AR,AI,N)
INTEGER K,L,N
DOUBLE PRECISION AR(21,N),AI(21,N), AROWR(21),AROWI(21)
C
DO 20 L=1,N
CALL ROW(AROWR,AROWI,L,N)
DO 10 K=1,21
AR(K,L) = AROWR(K)
AI(K,L) = AROWI(K)
10 CONTINUE
20 CONTINUE
C
RETURN
END
C
C
C
C
C
SUBROUTINE ROW(AROWR,AROWI,L,N)
INTEGER I,ISIGN,ISIGNI,ITEMP,L,MD,MDMI,MDPI,N
DOUBLE PRECISION AROWR(*),AROWI(*),ATEMP
C
ISIGN = -1
MD = 11
C
ITEMP = L/2
ATEMP = L - ITEMP*2
C
AROWR(MD) = 20.D0 * ATEMP
AROWI(MD) = 0.D0
C
DO 10 I = 1,10
ISIGN = ISIGN* (-1)
MDMI = MD - I
MDPI = MD + I
C
ISIGNI = ISIGN*I
C
AROWR(MDMI) = -ISIGNI
AROWR(MDPI) = ISIGNI
C
AROWI(MDMI) = ISIGNI
AROWI(MDPI) = -ISIGNI
10 CONTINUE
C
IF (L.GT.10) GO TO 30
C
IF (L.LE.10) THEN
AROWR(1) = 0.D0
AROWI(1) = 0.D0
ENDIF
IF (L.LE.9) THEN
AROWR(2) = 0.D0
AROWI(2) = 0.D0
ENDIF
IF (L.LE.8) THEN
AROWR(3) = 0.D0
AROWI(3) = 0.D0
ENDIF
IF (L.LE.7) THEN
AROWR(4) = 0.D0
AROWI(4) = 0.D0
ENDIF
IF (L.LE.6) THEN
AROWR(5) = 0.D0
AROWI(5) = 0.D0
ENDIF
IF (L.LE.5) THEN
AROWR(6) = 0.D0
AROWI(6) = 0.D0
ENDIF
IF (L.LE.4) THEN
AROWR(7) = 0.D0
AROWI(7) = 0.D0
ENDIF
IF (L.LE.3) THEN
AROWR(8) = 0.D0
AROWI(8) = 0.D0
ENDIF
IF (L.LE.2) THEN
AROWR(9) = 0.D0
AROWI(9) = 0.D0
ENDIF
IF (L.GT.1) GO TO 30
AROWR(10) = 0.D0
AROWI(10) = 0.D0
C
30 IF (L.LT.N-9) GO TO 40
C
IF (L.GE.N) THEN
AROWR(12) = 0.D0
AROWI(12) = 0.D0
ENDIF
IF (L.GE.N-1) THEN
AROWR(13) = 0.D0
AROWI(13) = 0.D0
ENDIF
IF (L.GE.N-2) THEN
AROWR(14) = 0.D0
AROWI(14) = 0.D0
ENDIF
IF (L.GE.N-3) THEN
AROWR(15) = 0.D0
AROWI(15) = 0.D0
ENDIF
IF (L.GE.N-4) THEN
AROWR(16) = 0.D0
AROWI(16) = 0.D0
ENDIF
IF (L.GE.N-5) THEN
AROWR(17) = 0.D0
AROWI(17) = 0.D0
ENDIF
IF (L.GE.N-6) THEN
AROWR(18) = 0.D0
AROWI(18) = 0.D0
ENDIF
IF (L.GE.N-7) THEN
AROWR(19) = 0.D0
AROWI(19) = 0.D0
ENDIF
IF (L.GE.N-8) THEN
AROWR(20) = 0.D0
AROWI(20) = 0.D0
ENDIF
IF (L.GE.N-9) THEN
AROWR(21) = 0.D0
AROWI(21) = 0.D0
ENDIF
40 CONTINUE
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: GIRI FORTRAN 77 JULY 17, 1990
C
C PURPOSE:
C THIS ALGORITHM IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN
C IMPLEMENTATION OF A GENERALIZED INVERSE RAYLEIGH ITERATION TO
C CALCULATE AN EIGENVALUE AND ITS ASSOCIATED RIGHT AND LEFT
C EIGENVECTORS OF A COMPLEX BAND MATRIX.
C
C THE ALGORITHM CONSISTS OF THE FOLLOWING STEPS:
C
C (0) CHOOSE INITIAL VECTORS U AND V FOR THE RIGHT AND LEFT
C EIGENVECTORS ASSOCIATED WITH THE APPROXIMATELY KNOWN
C EIGENVALUE OMEGA.
C
C (1) SOLVE THE TWO LINEAR SYSTEMS:
C
C ( A - OMEGA *I ) * X = U
C Y * ( A - OMEGA *I ) = V
C
C (2) CALCULATE THE IMPROVED APPROXIMATION OF THE EIGENVALUE:
C
C ( Y, A*X )
C OMEGA = OMEGA + ----------
C ( Y, X )
C
C (3) NORMALIZE X AND Y:
C
C X Y
C X = ------, Y = ------
C // X // // Y //
C
C (4) GO BACK TO STEP (1) IF THE CONVERGENCE CRITERION IS NOT
C SATISFIED, I.E. IF THE INCREMENT OF STEP (2) IS TOO LARGE.
C
C HEREIN, ( , ) IS THE COMPLEX EUCLEDIAN SCALAR PRODUCT AND //.//
C IS A NORM FOR COMPLEX VECTORS. THE TWO LINEAR SYSTEMS IN STEP (1)
C ARE SOLVED BY CALCULATING ONLY ONE LU-DECOMPOSITION OF THE MATRIX
C A - OMEGA*I.
C
C THE ALGORITHM USES REAL, DOUBLE-PRECISION ARITHMETIC. ALL COMPLEX
C VARIABLES ARE EXPRESSED BY THEIR REAL AND IMAGINARY PARTS.
C
C
C THE ELEMENTS OF THE REAL AND THE IMAGINARY PART OF THE BAND MATRIX
C ARE STORED IN ROWS AS ILLUSTRATED IN THE FOLLOWING EXAMPLE:
C
C
C
C * * 31 41 51 0 0 0 0 0 0
C * 22 32 42 52 0 0 0 0 0
C 13 23 33 43 53 0 0 0 0
C 0 14 24 34 44 54 0 0 0
C 0 0 15 25 35 45 55 0 0
C 0 0 0 16 26 36 46 56 0
C 0 0 0 0 17 27 37 47 57
C 0 0 0 0 0 18 28 38 48 *
C 0 0 0 0 0 0 19 29 39 * *
C
C
C THE FIRST INDEX SHOWS THE POSITION OF THE ELEMENT IN THE BAND,
C AND THE SECOND ONE INDICATES THE ROW. ONLY THE BAND IS STORED,
C AND THE USER HAS TO SUPPLY ZEROS FOR THE POSITIONS INDICATED WITH
C AN ASTERISK IN THE ABOVE EXAMPLE. THE MATRICES AR AND AI DECLARED
C IN GIRI CONTAIN THE FOLLOWING ELEMENTS:
C
C 0 0 31 41 51
C 0 22 32 42 52
C 13 23 33 43 53
C 14 24 34 44 54
C 15 25 35 45 55
C 16 26 36 46 56
C 17 27 37 47 57
C 18 28 38 48 0
C 19 29 39 0 0
C
C
C THE ELEMENTS OF THE ROWS OF AR AND AI ARE STORED IN CONSECUTIVE
C ADDRESSES SO THAT THE PROGRAM VECTORIZES ON MOST MACHINES.
C
C CALLING SEQUENCE:
C
C CALL GIRI( AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
C INIT,WR,WI,OMEGAR,OMEGAI,EPS,NSIMPL,LSTINV,
C IERROR,INFO,IUNIT )
C
C PARAMETERS:
C
C AR,AI: DOUBLE PRECISION, INPUT. AR(M,N), AI(M,N) CONTAIN THE
C REAL AND IMAGINARY PARTS OF THE MATRIX. THEIR VALUES
C REMAIN UNCHANGED.
C
C AWR,AWI: DOUBLE PRECISION. AWR(M,N) ,AWI(M,N) ARE TEMPORARY
C ARRAYS USED TO STORE REAL AND IMAGINARY PARTS OF THE
C UPPER-TRIANGULAR MATRIX OF THE LU-DECOMPOSITION OF
C THE MATRIX A - OMEGA*I.
C
C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C BAND. N MUST BE AN ODD NUMBER.
C
C N: INTEGER, INPUT, NUMBER OF EQUATIONS.
C
C XLR,XLI: DOUBLE PRECISION. XLR((M-1)/2,N) ,XLI((M-1)/2,N) ARE
C TEMPORARY ARRAYS USED TO STORE REAL AND IMAGINARY PARTS
C OF THE PSEUDO-LOWER-TRIANGULAR MATRIX OF THE LU-DECOM-
C POSITION OF THE MATRIX A - OMEGA*I.
C
C JP: INTEGER. JP(N) IS A TEMPORARY ARRAY USED TO STORE THE
C PIVOTING INFORMATION.
C
C UR,UI: DOUBLE PRECISION. UR(N),UI(N) CONTAIN REAL AND IMAGI-
C NARY PARTS OF THE INITIAL APPROXIMATION FOR THE RIGHT
C EIGENVECTOR AS INPUT (USED IF INIT=1) AND THE REAL
C AND IMAGINARY PARTS OF THE COMPUTED RIGHT EIGENVECTOR
C AS OUTPUT.
C
C VR,VI: DOUBLE PRECISION. VR(N),VI(N) CONTAIN REAL AND IMAGI-
C NARY PARTS OF THE INITIAL APPROXIMATION FOR THE LEFT
C EIGENVECTOR AS INPUT (USED IF INIT=1) AND THE REAL
C AND IMAGINARY PARTS OF THE COMPUTED LEFT EIGENVECTOR
C AS OUTPUT.
C
C INIT: INTEGER, INPUT:
C INIT=0: NO INITIAL VECTORS ARE PROVIDED. THE VECTORS
C UR = UI = VR = VI = (1, ... ,1) ARE TAKEN AS
C INITIAL VECTORS FOR THE ITERATIVE COMPUTATION
C OF THE RIGHT AND LEFT EIGENVECTORS.
C INIT=1: INITIAL APPROXIMATIONS FOR THE RIGHT AND LEFT
C EIGENVECTORS ARE PROVIDED. THEIR REAL AND
C IMAGINARY PARTS ARE STORED IN THE ARRAYS
C UR,UI,VR,VI.
C
C WR,WI: DOUBLE PRECISION. WR(N),WI(N) ARE WORK ARRAYS.
C
C OMEGAR, DOUBLE PRECISION. OMEGAR AND OMEGAI CONTAIN AN INITIAL
C OMEGAI: APPROXIMATION OF THE REAL AND IMAGINARY PARTS OF THE
C EIGENVALUE ON INPUT AND THE REAL AND IMAGINARY PARTS
C OF THE COMPUTED EIGENVALUE ON OUTPUT.
C
C EPS: DOUBLE PRECISION, INPUT. ABSOLUTE ERROR TOLERANCE FOR
C THE COMPUTED EIGENVALUE: THE ITERATION IS TERMINATED IF
C THE ABSOLUTE VALUE OF THE COMPUTED INCREMENT FOR THE
C EIGENVALUE IS LOWER THAN EPS. AN APPROPRIATE VALUE FOR
C EPS IS, FOR EXAMPLE, THE SQUARE ROOT OF MACHEP WHICH IS
C THE SMALLEST DOUBLE-PRECISION CONSTANT FOR WHICH
C 1.D0 + MACHEP .GT. 1.D0 HOLDS.
C
C NSIMPL: INTEGER, INPUT. NUMBER OF SIMPLIFIED ITERATION STEPS
C WITHOUT CALCULATING A NEW LU-DECOMPOSITION FOLLOWING
C A FULL ITERATION STEP. THE RECOMMENDED VALUE FOR NSIMPL
C IS 2.
C
C LAST: INTEGER, INPUT. NUMBER OF THE LAST ITERATION STEP THAT
C IS FOLLOWED BY NSIMPL SIMPLIFIED ITERATION STEPS. THE
C RECOMMENDED VALUE FOR LAST IS 1.
C
C IERROR: INTEGER, OUTPUT, ERROR CODE:
C IERROR=0: NO ERRORS OCCURRED.
C IERROR=1: M IS NOT ODD.
C IERROR=2: THE ITERATION DID NOT CONVERGE.
C
C INFO: INTEGER, INPUT.
C INFO = 0: CONVERGENCE INFORMATION IS NOT DISPLAYED.
C INFO = 1: CONVERGENCE INFORMATION IS DISPLAYED.
C
C IUNIT: INTEGER, INPUT. UNIT NUMBER ON THAT THE CONVERGENCE
C INFORMATION IS WRITTEN.
C
C
C GIRI CONTAINS CALLS TO THE FOLLOWING SUBROUTINES THAT ARE SUPPLIED
C WITH THE PROGRAM:
C
C GCBLU, TO PERFORM A LU-DECOMPOSITION OF A COMPLEX BAND MATRIX,
C GCSOL, TO SOLVE THE LINEAR SYSTEM USING THE LU-DECOMPOSTION,
C GCSOLA, TO SOLVE THE LINEAR SYSTEM WITH THE ADJOINT MATRIX
C WITHOUT CALCULATING A NEW LU-DECOMPOSITION,
C NORMAL, TO NORMALIZE A COMPLEX VECTOR,
C DOTPR, TO CALCULATE THE COMPLEX EUCLIDIAN SCALAR PRODUCT, AND
C CDIV, TO DIVIDE TWO COMPLEX NUMBERS.
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE GIRI(AR,AI,AWR,AWI,M,N,XLR,XLI,JP,UR,UI,VR,VI,
. INIT,WR,WI,OMEGAR,OMEGAI,EPS,NSIMPL,LAST,
. IERROR,INFO,IUNIT)
INTEGER I,IERROR,INFO,INIT,IPOS,IUNIT,J,LAST,
. M,MD,MNGIRI,N,NCOUNT,NGIRI,NSIMPL
INTEGER JP(N)
DOUBLE PRECISION AR(M,N),AI(M,N),AWR(M,N),AWI(M,N),
. XLR((M-1)/2,N),XLI((M-1)/2,N),
. UR(N),UI(N),VR(N),VI(N),WR(N),WI(N),
. DELTA,DELTAR,DELTAI,EPS,OMEGAR,OMEGAI,
. PIVMAX,PIVMIN,VUR,VUI,VWR,VWI
C
IERROR = 0
C MAXIMAL NUMBER OF GENERALIZED INVERSE RAYLEIGH ITERATION STEPS.
MNGIRI = 10
C CHECK ODDNESS OF M.
IF ( ((M+1)/2)*2 .NE. M+1 ) THEN
IERROR = 1
RETURN
ENDIF
C
MD = (M+1)/2
C
IF ( INIT .EQ. 1 ) GO TO 20
DO 10 J=1,N
UR(J) = 1.D0
UI(J) = 1.D0
VR(J) = 1.D0
VI(J) = 1.D0
10 CONTINUE
20 CONTINUE
C
DO 90 NGIRI =1,MNGIRI
C
DO 30 J=1,N
AR(MD,J) = AR(MD,J) - OMEGAR
AI(MD,J) = AI(MD,J) - OMEGAI
30 CONTINUE
C
DO 50 J=1,N
DO 40 I=1,M
AWR(I,J) = AR(I,J)
AWI(I,J) = AI(I,J)
40 CONTINUE
50 CONTINUE
C
CALL GCBLU(AWR,AWI,M,N,XLR,XLI,JP,PIVMAX,PIVMIN,IPOS)
C
DO 60 NCOUNT=0,NSIMPL
C
CALL GCSOL (AWR,AWI,M,N,XLR,XLI,UR,UI,JP)
CALL GCSOLA(AWR,AWI,M,N,XLR,XLI,VR,VI,JP)
C
CALL NORMAL(UR,UI,N)
CALL NORMAL(VR,VI,N)
C
CALL MATVEC(AR,AI,UR,UI,WR,WI,M,N)
C
CALL DOTPR(VR,VI,WR,WI,N,VWR,VWI)
CALL DOTPR(VR,VI,UR,UI,N,VUR,VUI)
C
CALL CDIV(VWR,VWI,VUR,VUI,DELTAR,DELTAI)
C
DELTA = SQRT( DELTAR**2 + DELTAI**2 )
C
C***** TO SUPRESS THE CONVERGENCE INFORMATION, THE
C***** FOLLOWING 20 LINES CAN BE COMMENTED OUT:
IF ( INFO .EQ. 1 ) THEN
IF ( NGIRI .EQ. 1 .AND. NCOUNT .EQ. 0 ) THEN
WRITE(IUNIT,900) OMEGAR,OMEGAI,
. NGIRI,NCOUNT,PIVMIN,PIVMAX,
. DELTAR,DELTAI,
. OMEGAR+DELTAR,OMEGAI+DELTAI
ELSE
WRITE(IUNIT,910) NGIRI,NCOUNT,PIVMIN,PIVMAX,
. DELTAR,DELTAI,
. OMEGAR+DELTAR,OMEGAI+DELTAI
ENDIF
ENDIF
C
900 FORMAT(//' N M PIVMIN PIVMAX ',
. ' INCREMENT EIGENVALUE'/
. 53(' '),' (',E10.3,',',E10.3,')'/
. 2I3,2E11.3,' (',E10.3,',',E10.3,')',
. ' (',E10.3,',',E10.3,')')
910 FORMAT(2I3,2E11.3,' (',E10.3,',',E10.3,')',
. ' (',E10.3,',',E10.3,')')
C*****
IF ( DELTA .LT. EPS ) GO TO 100
IF ( NGIRI .GT. LAST ) GO TO 70
60 CONTINUE
70 CONTINUE
C
DO 80 J=1,N
AR(MD,J) = AR(MD,J) + OMEGAR
AI(MD,J) = AI(MD,J) + OMEGAI
80 CONTINUE
C
OMEGAR = OMEGAR + DELTAR
OMEGAI = OMEGAI + DELTAI
C
90 CONTINUE
C ITERATION DID NOT CONVERGE WITHIN MNGIRI STEPS.
IERROR = 2
RETURN
C ITERATION CONVERGED. RESTORE AR AND AI.
100 DO 110 J=1,N
AR(MD,J) = AR(MD,J) + OMEGAR
AI(MD,J) = AI(MD,J) + OMEGAI
110 CONTINUE
C
OMEGAR = OMEGAR + DELTAR
OMEGAI = OMEGAI + DELTAI
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: GCBLU ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C GCBLU IS A 1977 AMERICAN NATIONAL STANDARD FORTRAN IMPLEMENTAION
C OF A LU-DECOMPOSITION OF A COMPLEX BAND MATRIX WITH PARTIAL
C PIVOTING.
C
C THE ELEMENTS OF THE REAL AND THE IMAGINARY PARTS OF THE BAND
C MATRIX ARE STORED IN THE ARRAYS AR AND AI AS DESCRIBED FOR GIRI.
C
C A SYSTEM WITH THIS MATRIX CAN BE SOLVED BY A SUBSEQUENT CALL OF
C THE SUBROUTINE GCSOL; A SYSTEM WITH THE ADJOINT MATRIX BY CALLING
C GCSOLA.
C
C CALLING SEQUENCE:
C
C CALL GCBLU(AR,AI,M,N,XLR,XLI,JP,PIVMAX,PIVMIN,IPOS)
C
C PARAMETERS:
C
C AR,AI: DOUBLE PRECISION. AR(M,N), AI(M,N) CONTAIN REAL AND
C IMAGINARY PARTS OF THE MATRIX ON INPUT AND ON OUTPUT
C REAL AND IMAGINARY PARTS OF THE UPPER TRIANGULAR
C MATRIX.
C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C BAND.
C N: INTEGER, INPUT, NUMBER OF EQUATIONS.
C XLR,XLI: DOUBLE PRECISION, OUTPUT. XLR((M-1)/2,N),XLI((M-1)/2,N)
C CONTAIN THE REAL AND IMAGINARY PARTS OF THE PSEUDO-
C LOWER-TRIANGULAR MATRIX OF THE LU-DECOMPOSITION.
C JP: INTEGER, OUTPUT. JP(N) CONTAINS PIVOTING IFORMATION.
C PIVMAX: DOUBLE PRECISION, OUTPUT, ( /RE(.)/+/IM(.)/ )-NORM
C OF THE LARGEST PIVOT ELEMENT.
C VALUES OF THE PIVOT ELEMENTS.
C PIVMIN: DOUBLE PRECISION, OUTPUT, ( /RE(.)/+/IM(.)/ )-NORM
C OF THE SMALLEST PIVOT.
C IPOS: INTEGER, OUTPUT, POSITION OF THE SMALLEST PIVOT.
C
C AUTHOR: GEZA SCHRAUF, CAL TECH/DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE GCBLU(AR,AI,M,N,XLR,XLI,JP,PIVMAX,PIVMIN,IPOS)
INTEGER I,IH,IPOS,ISTART,J,JPABS,JPREL,L,L1,L2,L3,LP1,
. M,M1,MD,MD1,N,N1
INTEGER JP(N)
DOUBLE PRECISION AR(M,N),AI(M,N),XLR((M-1)/2,N),XLI((M-1)/2,N),
. ABSPIV,EPS,PIVMAX,PIVMIN,RPIVR,RPIVI,SR,SI,TEMP
C EPS REPLACES ZERO PIVOT.
EPS = 1.D-20
C
MD = (M + 1)/2
MD1 = MD - 1
M1 = M - 1
N1 = N - 1
C SHIFT THE ELEMENTS OF THE FIRST MD1 ROWS TO THE LEFT.
IH = MD
ISTART = IH
C
DO 30 J=1,MD1
DO 10 I=1,IH
AR(I,J) = AR(ISTART+I-1,J)
AI(I,J) = AI(ISTART+I-1,J)
10 CONTINUE
IH = IH + 1
DO 20 I=IH,M
AR(I,J)=0.D0
AI(I,J)=0.D0
20 CONTINUE
ISTART = ISTART - 1
30 CONTINUE
C INITIALIZE PIVMAX AND PIVMIN.
PIVMAX = 0.D0
PIVMIN = SQRT( AR(1,1)**2 + AI(1,1)**2 )
C
DO 90 L=1,N1
L1 = L - 1
L2 = MIN(L+MD1,N)
L3 = L2 - L1
C PIVOT SEARCH.
ABSPIV = 0.D0
DO 40 J=1,L3
TEMP = ABS( AR(1,J+L1) ) + ABS( AI(1,J+L1) )
IF ( TEMP.LE.ABSPIV ) GOTO 40
JPREL = J
ABSPIV = TEMP
40 CONTINUE
IF ( ABSPIV .GT. PIVMAX ) PIVMAX = ABSPIV
IF ( ABSPIV .GE. PIVMIN ) GO TO 50
PIVMIN = ABSPIV
IPOS = L
50 CONTINUE
C
JP(L) = JPREL
JPABS = JPREL + L1
C INTERCHANGE THE JPABS-TH AND THE L-TH ROW.
DO 60 I=1,M
SR = AR(I,JPABS)
SI = AI(I,JPABS)
AR(I,JPABS) = AR(I,L)
AI(I,JPABS) = AI(I,L)
AR(I,L) = SR
AI(I,L) = SI
60 CONTINUE
C
C ADD A SUITABLE MULTIPLE OF OF THE L-TH ROW TO THE FOLLOWING
C ROWS, AND SHIFT THOSE ROWS TO THE LEFT, SO THAT THE NEXT PI-
C VOT IS TO BE SEARCHED AGAIN IN THE FIRST COLUMN OF A.
C
TEMP = AR(1,L)**2 + AI(1,L)**2
C ZERO PIVOT IS REPLACED BY EPS.
IF ( TEMP .EQ. 0.D0 ) TEMP = EPS
C
TEMP = 1.D0 / TEMP
RPIVR = AR(1,L) * TEMP
RPIVI =-AI(1,L) * TEMP
C
IH = 0
LP1 = L + 1
DO 80 J=LP1,L2
SR = -AR(1,J)*RPIVR + AI(1,J)*RPIVI
SI = -AI(1,J)*RPIVR - AR(1,J)*RPIVI
IH = IH + 1
XLR(IH,L) = SR
XLI(IH,L) = SI
DO 70 I=1,M1
AR(I,J) = AR(I+1,J) + AR(I+1,L)*SR - AI(I+1,L)*SI
AI(I,J) = AI(I+1,J) + AI(I+1,L)*SR + AR(I+1,L)*SI
70 CONTINUE
AR(M,J) = 0.D0
AI(M,J) = 0.D0
80 CONTINUE
90 CONTINUE
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: GCSOL ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C GCSOL SOLVES A SYSTEM WITH A COMPLEX BAND MATRIX, THE LU-DECOMPO-
C SITION OF WHICH HAS BEEN COMPUTED USING GCBLU.
C
C CALLING SEQUENCE:
C
C CALL GCSOL(AR,AI,M,N,XLR,XLI,BR,BI,JP)
C
C PARAMETERS:
C
C AR,AI: DOUBLE PRECISION, INPUT. AR(M,N), AI(M,N) CONTAIN REAL
C AND IMAGINARY PARTS OF THE UPPER TRIANGULAR MATRIX AS
C COMPUTED BY GCBLU. THEIR VALUES REMAIN UNCHANGED.
C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C BAND.
C N: INTEGER, INPUT, NUMBER OF EQUATIONS.
C XLR,XLI: DOUBLE PRECISION, INPUT. XLR((M-1)/2,N) ,XLI((M-1)/2,N)
C CONTAIN THE REAL AND IMAGINARY PARTS OF THE PSEUDO-
C LOWER-TRIANGULAR MATRIX OF THE LU-DECOMPOSITION.
C THEIR VALUES REMAIN UNCHANGED.
C BR,BI: DOUBLE PRECISION. BR(N), BI(N) ARE THE REAL AND IMAGI
C NARY PARTS OF THE RIGHT HAND SIDE ON INPUT AND OF THE
C SOLUTION VECTOR ON OUTPUT.
C JP: INTEGER, INPUT. JP(N) CONTAINS THE PIVOTING INFORMATION
C
C AUTHOR: GEZA SCHRAUF, CAL TECH/DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE GCSOL(AR,AI,M,N,XLR,XLI,BR,BI,JP)
INTEGER I,JPABS,L,L2,LEN,LH,LI,M,M1,N,N1
INTEGER JP(N)
DOUBLE PRECISION AR(M,N),AI(M,N),XLR((M-1)/2,N),XLI((M-1)/2,N),
. BR(N),BI(N),
. RA1LR,RA1LI,SR,SI,TEMP,TEMPR,TEMPI
C
M1 = M - 1
N1 = N - 1
C
C INTERCHANGE THE ELEMENTS OF B CORRESPONDING TO THE PERMUTATIONS
C OF THE ROWS OF A.
DO 20 L=1,N1
JPABS = L-1 + JP(L)
SR = BR(JPABS)
SI = BI(JPABS)
BR(JPABS) = BR(L)
BI(JPABS) = BI(L)
BR(L) = SR
BI(L) = SI
C AND ADD A MULTIPLE OF B(L) TO B(L+1),B(L+2), ... ,B(L+L2).
L2 = MIN((M-1)/2,N-L)
DO 10 I=1,L2
LI = L + I
BR(LI) = BR(LI) + XLR(I,L)*SR - XLI(I,L)*SI
BI(LI) = BI(LI) + XLI(I,L)*SR + XLR(I,L)*SI
10 CONTINUE
20 CONTINUE
C SOLVE THE SYSTEM WITH THE UPPER TRIANGULAR MATRIX.
C
TEMP = 1.D0 / ( AR(1,N)**2 + AI(1,N)**2 )
SR = ( BR(N)*AR(1,N) + BI(N)*AI(1,N) ) * TEMP
SI = ( BI(N)*AR(1,N) - BR(N)*AI(1,N) ) * TEMP
BR(N) = SR
BI(N) = SI
C
DO 40 LH=1,N1
LEN = MIN(LH,M1)
L = N - LH
SR = 0.D0
SI = 0.D0
DO 30 I=1,LEN
TEMPR = AR(I+1,L)*BR(I+L) - AI(I+1,L)*BI(I+L)
TEMPI = AI(I+1,L)*BR(I+L) + AR(I+1,L)*BI(I+L)
SR = SR + TEMPR
SI = SI + TEMPI
30 CONTINUE
C
TEMP = 1.D0 / ( AR(1,L)**2 + AI(1,L)**2 )
RA1LR = AR(1,L) * TEMP
RA1LI =-AI(1,L) * TEMP
C
TEMPR = ( BR(L) - SR )*RA1LR - ( BI(L) - SI )*RA1LI
TEMPI = ( BR(L) - SR )*RA1LI + ( BI(L) - SI )*RA1LR
BR(L) = TEMPR
BI(L) = TEMPI
40 CONTINUE
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: GCSOLA ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C GCSOL SOLVES A SYSTEM WITH A COMPLEX BAND MATRIX, THE LU-DECOMPO-
C SITION OF WHICH HAS BEEN COMPUTED USING GCBLU.
C
C CALLING SEQUENCE:
C
C CALL GCSOLA(AR,AI,M,N,XLR,XLI,BR,BI,JP)
C
C PARAMETERS:
C
C AR,AI: DOUBLE PRECISION, INPUT. AR(M,N), AI(M,N) CONTAIN THE
C REAL AND IMAGINARY PARTS OF THE UPPER TRIANGULAR MATRIX
C AS COMPUTED BY GCBLU. THEIR VALUES REMAIN UNCHANGED.
C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C BAND.
C N: INTEGER, INPUT, NUMBER OF EQUATIONS.
C XLR,XLI: DOUBLE PRECISION, INPUT. XLR((M-1)/2,N) ,XLI((M-1)/2,N)
C CONTAIN REAL AND IMAGINARY PARTS OF THE PSEUDO-
C LOWER-TRIANGULAR MATRIX OF THE LU-DECOMPOSITION.
C THEIR VALUES STAY UNCHANGED.
C BR,BI: DOUBLE PRECISION. BR(N), BI(N) ARE THE REAL AND IMAGI-
C NARY PARTS OF THE RIGHT HAND SIDE ON INPUT AND OF THE
C SOLUTION VECTOR ON OUTPUT.
C JP: INTEGER, INPUT. JP(N) CONTAINS THE PIVOTING INFORMATION
C
C AUTHOR: GEZA SCHRAUF, CAL TECH/DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE GCSOLA(AR,AI,M,N,XLR,XLI,BR,BI,JP)
INTEGER I,I1,IH,J,JH,JJ,JPABS,L,L2,LL,M,N,NM1,NP1
INTEGER JP(N)
DOUBLE PRECISION AR(M,N),AI(M,N),XLR((M-1)/2,N),XLI((M-1)/2,N),
. BR(N),BI(N),
. RA1LR,RA1LI,SR,SI,TEMP,TEMPR,TEMPI
C
NP1 = N + 1
C TRANSPOSE BY SHIFT.
DO 30 I=2,M
I1 = I - 1
JH = NP1 - I
IF (JH.EQ.0) GOTO 40
DO 10 JJ=1,JH
J = NP1 - JJ
AR(I,J) = AR(I,J-I1)
AI(I,J) = AI(I,J-I1)
10 CONTINUE
JH = JH + 1
DO 20 JJ = JH,N
J = NP1 - JJ
AR(I,J) = 0.D0
AI(I,J) = 0.D0
20 CONTINUE
30 CONTINUE
40 CONTINUE
C SOLVE TRANS(U)*Y = B.
C
TEMP = 1.D0 / ( AR(1,1)**2 + AI(1,1)**2 )
TEMPR = ( BR(1)*AR(1,1) - BI(1)*AI(1,1) )*TEMP
TEMPI = ( BI(1)*AR(1,1) + BR(1)*AI(1,1) )*TEMP
BR(1) = TEMPR
BI(1) = TEMPI
C
DO 60 L=2,N
L2 = MIN0(L,M)
SR = 0.D0
SI = 0.D0
DO 50 I=2,L2
IH = L+1 - I
TEMPR = AR(I,L)*BR(IH) + AI(I,L)*BI(IH)
TEMPI = - AI(I,L)*BR(IH) + AR(I,L)*BI(IH)
SR = SR + TEMPR
SI = SI + TEMPI
50 CONTINUE
C
TEMP = 1.D0 / ( AR(1,L)**2 + AI(1,L)**2 )
RA1LR = AR(1,L) * TEMP
RA1LI = AI(1,L) * TEMP
C
TEMPR = ( BR(L) - SR )*RA1LR - ( BI(L) - SI )*RA1LI
TEMPI = ( BI(L) - SI )*RA1LR + ( BR(L) - SR )*RA1LI
BR(L) = TEMPR
BI(L) = TEMPI
C
60 CONTINUE
C SHIFT BACK IN ORDER TO RESTORE THE MATRIX.
DO 90 I=2,M
I1 = I - 1
JH = NP1 - I
IF (JH.EQ.0) GOTO 100
DO 70 J=1,JH
AR(I,J) = AR(I,J+I1)
AI(I,J) = AI(I,J+I1)
70 CONTINUE
JH = JH + 1
DO 80 J = JH,N
AR(I,J) = 0.D0
AI(I,J) = 0.D0
80 CONTINUE
90 CONTINUE
100 CONTINUE
C SOLVE TRANS(XL)*X = Y.
NM1 = N - 1
DO 120 LL=1,NM1
L = N - LL
L2 = MIN0(LL,(M-1)/2)
SR = 0.D0
SI = 0.D0
DO 110 I=1,L2
IH = L + I
TEMPR = XLR(I,L)*BR(IH) + XLI(I,L)*BI(IH)
TEMPI = - XLI(I,L)*BR(IH) + XLR(I,L)*BI(IH)
SR = SR + TEMPR
SI = SI + TEMPI
110 CONTINUE
BR(L) = BR(L) + SR
BI(L) = BI(L) + SI
C
JPABS = L-1 + JP(L)
TEMPR = BR(L)
TEMPI = BI(L)
BR(L) = BR(JPABS)
BI(L) = BI(JPABS)
BR(JPABS) = TEMPR
BI(JPABS) = TEMPI
120 CONTINUE
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: MATVEC ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C MATVEC COMPUTES Y = A*X FOR A COMPLEX BAND MATRIX A AND A COMPLEX
C VECTOR X.
C
C CALLING SEQUENCE:
C
C CALL MATVEC(AR,AI,XR,XI,YR,YI,M,N)
C
C PARAMETERS:
C AR,AI: DOUBLE PRECISION, INPUT. AR(M,N), AI(M,N) CONTAIN REAL
C AND THE IMAGINARY PARTS OF THE MATRIX.
C XR,XI: DOUBLE PRECISION, INPUT. XR(N), XI(N) REAL AND IMAGI-
C NARY PARTS OF X.
C YR,YI: DOUBLE PRECISION, OUTPUT. YR(N), YI(N) REAL AND IMAGI-
C NARY PARTS OF Y.
C M: INTEGER, INPUT, NUMBER OF ELEMENTS IN EACH ROW OF THE
C BAND.
C N: INTEGER, INPUT, DIMENSION OF THE MATRIX.
C
C AUTHOR: GEZA SCHRAUF, CAL TECH/DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE MATVEC(AR,AI,XR,XI,YR,YI,M,N)
INTEGER K,KEND,KSTART,L,L0,LL,LMMD,M,MD,N
DOUBLE PRECISION AR(M,N),AI(M,N),XR(N),XI(N),YR(N),YI(N),
. SR,SI,TEMPR,TEMPI
C
MD = (M+1)/2
C
DO 20 L=1,N
LMMD = L - MD
KSTART = MAX(1,1-LMMD)
KEND = MIN(M,N-LMMD)
L0 = MAX(0,LMMD) + 1
SR = 0.D0
SI = 0.D0
DO 10 K=KSTART,KEND
LL = K - KSTART + L0
TEMPR = AR(K,L)*XR(LL) - AI(K,L)*XI(LL)
TEMPI = AI(K,L)*XR(LL) + AR(K,L)*XI(LL)
SR = SR + TEMPR
SI = SI + TEMPI
10 CONTINUE
YR(L) = SR
YI(L) = SI
20 CONTINUE
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: DOTPR ANSI FORTRAN 77 MAY 15, 1900
C
C PURPOSE:
C DOTPR COMPUTES THE COMPLEX EUCLEDIAN SCALAR PRODUCT OF TWO COMPLEX
C VECTORS USING LOOP-UNROLLING.
C
C CALLING SEQUENCE:
C
C CALL DOTPR(XR,XI,YR,YI,N,SPR,SPI)
C
C PARAMETERS:
C XR,XI: DOUBLE PRECISION, INPUT. XR(N), XI(N) REAL AND IMAGI-
C NARY PARTS OF THE FIRST VECTOR X.
C YR,YI: DOUBLE PRECISION, INPUT. YR(N), YI(N) REAL AND IMAGI-
C NARY PARTS OF THE SECOND VECTOR Y.
C N: INTEGER, INPUT, DIMENSION OF THE VECTORS.
C SPR,SPI: DOUBLE PRECISION, OUTPUT, REAL AND IMAGINARY PARTS OF
C THE SCALAR PRODUCT.
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE DOTPR(XR,XI,YR,YI,N,SPR,SPI)
INTEGER I,M,MP1,N
DOUBLE PRECISION XR(N),XI(N),YR(N),YI(N), SPR,SPI,TEMPR,TEMPI
C
SPR = 0.D0
SPI = 0.D0
C
TEMPR = 0.D0
TEMPI = 0.D0
C
M = MOD(N,5)
IF ( M .EQ. 0 ) GO TO 20
C
DO 10 I=1,M
TEMPR = TEMPR + XR(I)*YR(I) + XI(I)*YI(I)
TEMPI = TEMPI - XI(I)*YR(I) + XR(I)*YI(I)
10 CONTINUE
C
IF ( N .LT. 5 ) GO TO 40
C
20 MP1 = M + 1
C
DO 30 I=MP1,N,5
TEMPR = TEMPR + XR(I)*YR(I) + XI(I)*YI(I)
. + XR(I+1)*YR(I+1) + XI(I+1)*YI(I+1)
. + XR(I+2)*YR(I+2) + XI(I+2)*YI(I+2)
. + XR(I+3)*YR(I+3) + XI(I+3)*YI(I+3)
. + XR(I+4)*YR(I+4) + XI(I+4)*YI(I+4)
TEMPI = TEMPI - XI(I)*YR(I) + XR(I)*YI(I)
. - XI(I+1)*YR(I+1) + XR(I+1)*YI(I+1)
. - XI(I+2)*YR(I+2) + XR(I+2)*YI(I+2)
. - XI(I+3)*YR(I+3) + XR(I+3)*YI(I+3)
. - XI(I+4)*YR(I+4) + XR(I+4)*YI(I+4)
30 CONTINUE
C
40 SPR = TEMPR
SPI = TEMPI
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: NORMAL ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C NORMAL NORMALIZES A COMPLEX VECTOR (USING LOOP-UNROLLING).
C
C CALLING SEQUENCE:
C
C CALL NORMAL(XR,XI,N)
C
C PARAMETERS:
C XR,XI: DOUBLE PRECISION, XR(N), XI(N). ON INPUT REAL AND
C IMAGINARY PARTS OF THE VECTOR TO BE NORMALIZED. ON
C OUTPUT REAL AND IMAGINARY PARTS OF THE NORMALIZED
C VECTOR.
C N: INTEGER, INPUT, DIMENSION OF THE VECTOR.
C
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE NORMAL(XR,XI,N)
INTEGER I,M,MP1,N
DOUBLE PRECISION XR(N),XI(N), RFAC,SUP,TEMP, TIR,TII,TIP1R,TIP1I,
. TIP2R,TIP2I,TIP3R,TIP3I,TIP4R,TIP4I
C
SUP = 0.D0
C
DO 10 I=1,N
TEMP = MAX( ABS(XR(I)),ABS(XI(I)) )
IF ( TEMP .GT. SUP ) SUP = TEMP
10 CONTINUE
C
RFAC = 1.D0 / SUP
C
M = MOD(N,5)
IF ( M .EQ. 0 ) GO TO 30
C
DO 20 I=1,M
TIR = RFAC * ( XR(I) - XI(I) )
TII = RFAC * ( XR(I) + XI(I) )
XR(I) = TIR
XI(I) = TII
20 CONTINUE
C
IF ( N .LT. 5 ) RETURN
C
30 MP1 = M + 1
C
DO 40 I=MP1,N,5
TIR = RFAC * ( XR(I) - XI(I) )
TII = RFAC * ( XR(I) + XI(I) )
XR(I) = TIR
XI(I) = TII
C
TIP1R = RFAC * ( XR(I+1) - XI(I+1) )
TIP1I = RFAC * ( XR(I+1) + XI(I+1) )
XR(I+1) = TIP1R
XI(I+1) = TIP1I
C
TIP2R = RFAC * ( XR(I+2) - XI(I+2) )
TIP2I = RFAC * ( XR(I+2) + XI(I+2) )
XR(I+2) = TIP2R
XI(I+2) = TIP2I
C
TIP3R = RFAC * ( XR(I+3) - XI(I+3) )
TIP3I = RFAC * ( XR(I+3) + XI(I+3) )
XR(I+3) = TIP3R
XI(I+3) = TIP3I
C
TIP4R = RFAC * ( XR(I+4) - XI(I+4) )
TIP4I = RFAC * ( XR(I+4) + XI(I+4) )
XR(I+4) = TIP4R
XI(I+4) = TIP4I
40 CONTINUE
C
RETURN
END
C
C
C
C
C
C***********************************************************************
C
C NAME: CDIV ANSI FORTRAN 77 JULY 16, 1990
C
C PURPOSE:
C DIVIDES TWO COMPLEX NUMBER: Z = X / Y.
C
C CALLING SEQUENCE:
C
C CALL CDIV(XR,XI,YR,YI,ZR,ZI)
C
C PARAMETERS:
C XR,XI: DOUBLE PRECISION, INPUT. REAL AND IMAGINARY PARTS OF
C THE DIVIDEND.
C YR,YI: DOUBLE PRECISION, INPUT. REAL AND IMAGINARY PARTS OF
C THE DIVSOR.
C ZR,ZI: DOUBLE PRECISION, OUTPUT. REAL AND IMAGINARY PARTS OF
C THE QUOTIENT.
C
C AUTHOR: GEZA SCHRAUF, DEUTSCHE AIRBUS GMBH
C
C***********************************************************************
SUBROUTINE CDIV(XR,XI,YR,YI,ZR,ZI)
DOUBLE PRECISION A,B,D,RFAC,XR,XI,YR,YI,ZR,ZI
C
IF ( ABS(YR) .GE. ABS(YI) ) THEN
RFAC = 1.D0 / YR
A = RFAC * XR
B = RFAC * XI
D = RFAC * YI
ELSE
RFAC = 1.D0 / YI
A = RFAC * XI
B =-RFAC * XR
D =-RFAC * YR
ENDIF
C
RFAC = 1.D0 / ( 1.D0 + D**2 )
C
ZR = RFAC * ( A + B*D )
ZI = RFAC * ( B - A*D )
C
RETURN
END
*
*
-----------------------------------------------------------------------
TEST OF THE GENERALIZED INVERSE RAYLEIGH ITERATION USING A COMPLEX B
MATRIX OF DIMENSION 500 AND BAND WIDTH 21
-----------------------------------------------------------------------
*
*
*
*
-----------------------------------------------------------------------
F I R S T T E S T: EIGENVALUE COMPUTATION USING INTERNAL INITIAL VE
-----------------------------------------------------------------------
*
*
*
INITIAL GUESS FOR EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.603E+02, 0.48
*
*
*
N M PIVMIN PIVMAX INCREMENT EIGENVALUE
( 0.603E+02, 0.48
1 0 0.418E+02 0.202E+03 (-0.499E+01,-0.294E+02) ( 0.553E+02, 0.18
1 1 0.418E+02 0.202E+03 ( 0.252E+00, 0.169E+00) ( 0.606E+02, 0.48
1 2 0.418E+02 0.202E+03 ( 0.252E+00, 0.170E+00) ( 0.606E+02, 0.48
2 0 0.435E+02 0.197E+03 ( 0.299E-07, 0.308E-07) ( 0.606E+02, 0.48
3 0 0.435E+02 0.197E+03 ( 0.579E-15,-0.178E-15) ( 0.606E+02, 0.48
*
*
*
COMPUTED EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.6055191E+02, 0.484697
*
*
*
ERROR CODE: IERROR = 0
*
*
*
*
-----------------------------------------------------------------------
S E C O N D T E S T: EIGENVALUE COMPUTATION USING THE CALCULATE
EIGENVECTORS AS INITIAL VECTORS
-----------------------------------------------------------------------
*
*
*
INITIAL GUESS FOR EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.603E+02, 0.48
*
*
*
N M PIVMIN PIVMAX INCREMENT EIGENVALUE
( 0.603E+02, 0.48
1 0 0.418E+02 0.202E+03 ( 0.252E+00, 0.170E+00) ( 0.606E+02, 0.48
1 1 0.418E+02 0.202E+03 ( 0.252E+00, 0.170E+00) ( 0.606E+02, 0.48
1 2 0.418E+02 0.202E+03 ( 0.252E+00, 0.170E+00) ( 0.606E+02, 0.48
2 0 0.435E+02 0.197E+03 ( 0.531E-15,-0.852E-15) ( 0.606E+02, 0.48
*
*
*
COMPUTED EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.6055191E+02, 0.484697
*
*
*
ERROR CODE: IERROR = 0
*
*
*
*
-----------------------------------------------------------------------
T H I R D T E S T: EIGENVALUE COMPUTATION WITHOUT THE SIMPLIFIED S
-----------------------------------------------------------------------
*
*
*
INITIAL GUESS FOR EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.603E+02, 0.48
*
*
*
N M PIVMIN PIVMAX INCREMENT EIGENVALUE
( 0.603E+02, 0.48
1 0 0.418E+02 0.202E+03 (-0.499E+01,-0.294E+02) ( 0.553E+02, 0.18
2 0 0.401E+02 0.772E+02 ( 0.286E+02, 0.382E+02) ( 0.839E+02, 0.57
3 0 0.857E+02 0.146E+03 (-0.220E+02,-0.899E+01) ( 0.619E+02, 0.48
4 0 0.478E+02 0.179E+03 (-0.139E+01, 0.349E+00) ( 0.606E+02, 0.48
5 0 0.435E+02 0.197E+03 ( 0.263E-03, 0.477E-03) ( 0.606E+02, 0.48
6 0 0.435E+02 0.197E+03 (-0.557E-13, 0.760E-12) ( 0.606E+02, 0.48
*
*
*
COMPUTED EIGENVALUE: (OMEGA_R,OMEGA_I) = ( 0.6055191E+02, 0.484697
*
*
*
ERROR CODE: IERROR = 0
*
*
*
*
-----------------------------------------------------------------------
F O U R T H T E S T: ERROR CODE "IERROR=1"
-----------------------------------------------------------------------
*
*
*
ERROR CODE: IERROR = 1
*
*
*
.