[CONTACT]

[ABOUT]

[POLICY]

[ADVERTISE]

C ALGORITHM COLLECTED ALGORITHMS FROM

Found at: ftp.icm.edu.pl:70/packages/netlib/toms-2014-06-10/696

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
*
*
*

		
.


AD:

NEW PAGES:

[ODDNUGGET]

[GOPHER]