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