C ALGORITHM 640 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.12, NO. 1, C MAR., 1986, P. 26. C PROGRAM SFREQT C C THIS IS A SAMPLE MAIN CALLING PROGRAM TO ILLUSTRATE AND TEST THE C USE OF SFRMG IN OBTAINING THE DATA TO PRODUCE A BODE PLOT FROM C A STATE SPACE REPRESENTATION OF A SINGLE-INPUT-SINGLE-OUTPUT C (M=1, L=1) CONTINUOUS-TIME LINEAR SYSTEM WITH AT MOST N=25 STATES C (NUMBER OF STATES IS SET BY THE PARAMETER NS BELOW). C THE CODE IS FORTRAN 77 COMPATIBLE BUT THE I/O MAY BE BETTER C MODIFIED FOR YOUR PARTICULAR SYSTEM C INTEGER NS,NIN,NOUT PARAMETER (NS=25) C C SET UNIT OR DEVICE NUMBERS C NIN UNIT NUMBER FOR READING THE INPUT DATA FROM FILE INTEST.DAT C NOUT UNIT NUMBER FOR WRITING THE RESULTS TO FILE OUTEST.DAT C PARAMETER (NIN=2) PARAMETER (NOUT=3) C INTEGER I,IERR,INFO,J,K,KOUNT,L,M,N,NA,NB,NC,NTESTS INTEGER IPVT(NS) REAL ALOGFR,AMAG,PHASE,RCOND,X,Y,Z REAL A(NS,NS),B(NS,1),C(1,NS),EVIM(NS),EVRE(NS),ZR(NS) COMPLEX FREQ COMPLEX AINVB(NS,1),G(1,1),H(NS,NS),ZC(NS) LOGICAL FIRST CHARACTER *50 FIN,FOUT C C FORTRAN FUNCTIONS CALLED C C REAL ABS,AIMAG,ATAN,REAL,SQRT C COMPLEX CMPLX C C READ LOCAL FILE OR PATH NAMES FROM UNIT=1. C THIS FILE HAS TWO NAMES: C INDATA C OUTDATA C THE FILE 'INDATA' IS READ BY THIS DRIVER PROGRAM. C THE FILE 'OUTDATA' IS WRITTEN BY THIS DRIVER PROGRAM. READ (1,'(A)') FIN,FOUT OPEN (UNIT=NIN,FILE=FIN,STATUS='OLD') OPEN (UNIT=NOUT,FILE=FOUT,STATUS='UNKNOWN') C C READ AN INTEGER NTESTS FROM FILE 'INTEST.DAT' WHICH C TELLS HOW MANY DIFFERENT TEST CASES WILL BE RUN C READ (NIN,*) NTESTS DO 130 KOUNT = 1,NTESTS WRITE (NOUT,9001) WRITE (NOUT,9011) KOUNT WRITE (NOUT,9001) FIRST = .TRUE. NA = NS NB = NS NC = 1 L = 1 M = 1 C C READ INPUT DATA FROM FILE 'INTEST.DAT' C N .LE. 25; SEE ABOVE C IERR,INFO SEE SFRMG DOCUMENTATION C A N X N STATE MATRIX; INPUT BY ROWS C B N-VECTOR; INPUT ITS TRANSPOSE AS A ROW VECTOR C C N-VECTOR; INPUT AS A ROW VECTOR C READ (NIN,*) N READ (NIN,*) IERR,INFO DO 10 I = 1,N READ (NIN,*) (A(I,J),J=1,N) 10 CONTINUE READ (NIN,*) (B(I,1),I=1,N) READ (NIN,*) (C(1,J),J=1,N) C FREQ = CMPLX(0.0E0,0.0E0) C C INITIAL CALL TO SFRMG PERFORMS THE HESSENBERG REDUCTION PLUS C WHATEVER OPTIONS ARE REQUESTED BY CHOICE OF THE PARAMETERS C IERR AND INFO AS INPUT C CALL SFRMG(FIRST,NA,NB,NC,L,M,N,A,B,C,FREQ,AINVB,G,H,IPVT,ZC, . ZR,EVRE,EVIM,IERR,INFO,RCOND) C C CHECK ERROR RETURN CONDITIONS, ETC. FOR INITIAL DATA C (I.E., ZERO FREQUENCY); SEE SFRMG DOCUMENTATION FOR DETAILS C THE REAL COMPUTING STARTS BELOW AFTER STATEMENT NO. 200 C IF (IERR.LT.0) GO TO 40 IF (IERR.EQ.0) GO TO 20 WRITE (NOUT,9021) IERR GO TO 40 20 CONTINUE WRITE (NOUT,9031) DO 30 I = 1,N WRITE (NOUT,*) EVRE(I),EVIM(I) 30 CONTINUE 40 CONTINUE IF (INFO.LE.0) GO TO 50 WRITE (NOUT,9041) INFO GO TO 80 50 IF (INFO.EQ.0) GO TO 70 IF ((1.0E0+RCOND).GT.1.0E0) GO TO 60 WRITE (NOUT,9051) GO TO 80 60 CONTINUE RCOND = 1.0E0/RCOND WRITE (NOUT,9061) RCOND 70 CONTINUE WRITE (NOUT,9071) WRITE (NOUT,*) G(1,1) WRITE (NOUT,9131) 80 CONTINUE C C READY TO COMPUTE FREQUENCY RESPONSE PLOT DATA; THIS EXAMPLE C COMPUTES MAGNITUDE AND PHASE OF G(FREQ) FOR VALUES OF C LOG10(FREQ) RANGING OVER 31 EQUALLY SPACED VALUES FROM -3 TO +3 C DO 120 K = 1,31 ALOGFR = 0.2E0*K - 3.2E0 WRITE (NOUT,9081) ALOGFR FREQ = CMPLX(0.0E0, (10.0E0**ALOGFR)) CALL SFRMG(FIRST,NA,NB,NC,L,M,N,A,B,C,FREQ,AINVB,G,H,IPVT, . ZC,ZR,EVRE,EVIM,IERR,INFO,RCOND) IF (INFO.LE.0) GO TO 90 WRITE (NOUT,9091) INFO GO TO 120 90 CONTINUE IF (INFO.EQ.0) GO TO 110 IF ((1.0E0+RCOND).GT.1.0E0) GO TO 100 WRITE (NOUT,9101) GO TO 120 100 CONTINUE RCOND = 1.0E0/RCOND WRITE (NOUT,9111) RCOND 110 CONTINUE X = REAL(G(1,1)) Y = AIMAG(G(1,1)) Z = ABS(X) + ABS(Y) X = X/Z Y = Y/Z AMAG = Z*SQRT(X*X+Y*Y) PHASE = ATAN(Y/X)* (45.0E0/ATAN(1.0E0)) WRITE (NOUT,9121) WRITE (NOUT,*) G(1,1),AMAG,PHASE WRITE (NOUT,9131) 120 CONTINUE 130 CONTINUE STOP 9001 FORMAT (/' **********') 9011 FORMAT (////' TEST CASE NUMBER',I4////) 9021 FORMAT (/' IERR =',I5/) 9031 FORMAT (/' THE EIGENVALUES OF A ARE:'//) 9041 FORMAT (/' INFO =',I5/) 9051 FORMAT (/' THE A MATRIX IS SINGULAR TO WORKING PRECISION'/) 9061 FORMAT (/' CONDITION ESTIMATE OF A (POSSIBLY BALANCED) IS:', . F10.6/) 9071 FORMAT (/' THE VALUE OF -C*(A-INVERSE)*B IS:'/) 9081 FORMAT (/' LOG10(FREQ) =',F7.3) 9091 FORMAT (/' INFO =',I5/) 9101 FORMAT (/' ((FREQ*I)-A) IS SINGULAR TO WORKING PRECISION'/) 9111 FORMAT (/' CONDITION ESTIMATE OF ((FREQ*I)-A) IS:',F10.6/) 9121 FORMAT (/' G(FREQ), ITS MAGNITUDE, AND ITS PHASE:'/) 9131 FORMAT (///) END C C ATTACHED ARE TWO EISPACK SUBBROUTINES (BALANC AND HQR) C WHICH MAY BE NECESSARY IF IT IS DESIRED TO COMPUTE THE C EIGENVALUES OF A AT THE INITIAL STEP (SEE SFRMG DOCUM.) C ***** NOTE: THESE ROUTINES CONTAIN MACHINE DEPENDENT C PARAMETERS (RADIX IN BALANC, MACHEP IN ORTHES) C WHICH MAY NEED MODIFICATION FOR YOUR SYSTEM ***** C SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE) C INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC REAL A(NM,N),SCALE(N) REAL C,F,G,R,S,B2,RADIX REAL ABS LOGICAL NOCONV C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). C C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES C EIGENVALUES WHENEVER POSSIBLE. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C A CONTAINS THE INPUT MATRIX TO BE BALANCED. C C ON OUTPUT- C C A CONTAINS THE BALANCED MATRIX, C C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) C IS EQUAL TO ZERO IF C (1) I IS GREATER THAN J AND C (2) J=1,...,LOW-1 OR I=IGH+1,...,N, C C SCALE CONTAINS INFORMATION DETERMINING THE C PERMUTATIONS AND SCALING FACTORS USED. C C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN C SCALE(J) = P(J), FOR J = 1,...,LOW-1 C = D(J,J), J = LOW,...,IGH C = P(J) J = IGH+1,...,N. C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, C THEN 1 TO LOW-1. C C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. C C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS C K,L HAVE BEEN REVERSED.) C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C ********** RADIX IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE BASE OF THE MACHINE FLOATING POINT REPRESENTATION. C C ********** RADIX = 2. C B2 = RADIX*RADIX K = 1 L = N GO TO 60 C ********** IN-LINE PROCEDURE FOR ROW AND C COLUMN EXCHANGE ********** 10 SCALE(M) = J IF (J.EQ.M) GO TO 40 C DO 20 I = 1,L F = A(I,J) A(I,J) = A(I,M) A(I,M) = F 20 CONTINUE C DO 30 I = K,N F = A(J,I) A(J,I) = A(M,I) A(M,I) = F 30 CONTINUE C 40 GO TO (50,90),IEXC C ********** SEARCH FOR ROWS ISOLATING AN EIGENVALUE C AND PUSH THEM DOWN ********** 50 IF (L.EQ.1) GO TO 230 L = L - 1 C ********** FOR J=L STEP -1 UNTIL 1 DO -- ********** 60 DO 80 JJ = 1,L J = L + 1 - JJ C DO 70 I = 1,L IF (I.EQ.J) GO TO 70 IF (A(J,I).NE.0.0) GO TO 80 70 CONTINUE C M = L IEXC = 1 GO TO 10 80 CONTINUE C GO TO 100 C ********** SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE C AND PUSH THEM LEFT ********** 90 K = K + 1 C 100 DO 120 J = K,L C DO 110 I = K,L IF (I.EQ.J) GO TO 110 IF (A(I,J).NE.0.0) GO TO 120 110 CONTINUE C M = K IEXC = 2 GO TO 10 120 CONTINUE C ********** NOW BALANCE THE SUBMATRIX IN ROWS K TO L ********** DO 130 I = K,L SCALE(I) = 1.0 130 CONTINUE C ********** ITERATIVE LOOP FOR NORM REDUCTION ********** 140 NOCONV = .FALSE. C DO 220 I = K,L C = 0.0 R = 0.0 C DO 150 J = K,L IF (J.EQ.I) GO TO 150 C = C + ABS(A(J,I)) R = R + ABS(A(I,J)) 150 CONTINUE C ********** GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ********** IF (C.EQ.0.0 .OR. R.EQ.0.0) GO TO 220 G = R/RADIX F = 1.0 S = C + R 160 IF (C.GE.G) GO TO 170 F = F*RADIX C = C*B2 GO TO 160 170 G = R*RADIX 180 IF (C.LT.G) GO TO 190 F = F/RADIX C = C/B2 GO TO 180 C ********** NOW BALANCE ********** 190 IF ((C+R)/F.GE.0.95*S) GO TO 220 G = 1.0/F SCALE(I) = SCALE(I)*F NOCONV = .TRUE. C DO 200 J = K,N A(I,J) = A(I,J)*G 200 CONTINUE C DO 210 J = 1,L A(J,I) = A(J,I)*F 210 CONTINUE C 220 CONTINUE C IF (NOCONV) GO TO 140 C 230 LOW = K IGH = L RETURN C ********** LAST CARD OF BALANC ********** END C C ------------------------------------------------------------------ C SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR) C INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITS,LOW,MP2,ENM2,IERR REAL H(NM,N),WR(N),WI(N) REAL P,Q,R,S,T,W,X,Y,ZZ,NORM,MACHEP REAL SQRT,ABS,SIGN INTEGER MIN0 LOGICAL NOTLAS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL C UPPER HESSENBERG MATRIX BY THE QR METHOD. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, C SET LOW=1, IGH=N, C C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. C C ON OUTPUT- C C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED, C C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT C FOR INDICES IERR+1,...,N, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** MACHEP = 2.** (-24) C IERR = 0 NORM = 0.0 K = 1 C ********** STORE ROOTS ISOLATED BY BALANC C AND COMPUTE MATRIX NORM ********** DO 20 I = 1,N C DO 10 J = K,N NORM = NORM + ABS(H(I,J)) 10 CONTINUE C K = I IF (I.GE.LOW .AND. I.LE.IGH) GO TO 20 WR(I) = H(I,I) WI(I) = 0.0 20 CONTINUE C EN = IGH T = 0.0 C ********** SEARCH FOR NEXT EIGENVALUES ********** 30 IF (EN.LT.LOW) GO TO 250 ITS = 0 NA = EN - 1 ENM2 = NA - 1 C ********** LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT C FOR L=EN STEP -1 UNTIL LOW DO -- ********** 40 DO 50 LL = LOW,EN L = EN + LOW - LL IF (L.EQ.LOW) GO TO 60 S = ABS(H(L-1,L-1)) + ABS(H(L,L)) IF (S.EQ.0.0) S = NORM IF (ABS(H(L,L-1)).LE.MACHEP*S) GO TO 60 50 CONTINUE C ********** FORM SHIFT ********** 60 X = H(EN,EN) IF (L.EQ.EN) GO TO 200 Y = H(NA,NA) W = H(EN,NA)*H(NA,EN) IF (L.EQ.NA) GO TO 210 IF (ITS.EQ.30) GO TO 240 IF (ITS.NE.10 .AND. ITS.NE.20) GO TO 80 C ********** FORM EXCEPTIONAL SHIFT ********** T = T + X C DO 70 I = LOW,EN H(I,I) = H(I,I) - X 70 CONTINUE C S = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75*S Y = X W = -0.4375*S*S 80 ITS = ITS + 1 C ********** LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS. C FOR M=EN-2 STEP -1 UNTIL L DO -- ********** DO 90 MM = L,ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R*S-W)/H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = ABS(P) + ABS(Q) + ABS(R) P = P/S Q = Q/S R = R/S IF (M.EQ.L) GO TO 100 IF (ABS(H(M,M-1))* (ABS(Q)+ABS(R)).LE.MACHEP*ABS(P)* . (ABS(H(M-1,M-1))+ABS(ZZ)+ABS(H(M+1,M+1)))) GO TO 100 90 CONTINUE C 100 MP2 = M + 2 C DO 110 I = MP2,EN H(I,I-2) = 0.0 IF (I.EQ.MP2) GO TO 110 H(I,I-3) = 0.0 110 CONTINUE C ********** DOUBLE QR STEP INVOLVING ROWS L TO EN AND C COLUMNS M TO EN ********** DO 190 K = M,NA NOTLAS = K .NE. NA IF (K.EQ.M) GO TO 120 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X.EQ.0.0) GO TO 190 P = P/X Q = Q/X R = R/X 120 S = SIGN(SQRT(P*P+Q*Q+R*R),P) IF (K.EQ.M) GO TO 130 H(K,K-1) = -S*X GO TO 140 130 IF (L.NE.M) H(K,K-1) = -H(K,K-1) 140 P = P + S X = P/S Y = Q/S ZZ = R/S Q = Q/P R = R/P C ********** ROW MODIFICATION ********** DO 160 J = K,EN P = H(K,J) + Q*H(K+1,J) IF ( .NOT. NOTLAS) GO TO 150 P = P + R*H(K+2,J) H(K+2,J) = H(K+2,J) - P*ZZ 150 H(K+1,J) = H(K+1,J) - P*Y H(K,J) = H(K,J) - P*X 160 CONTINUE C J = MIN0(EN,K+3) C ********** COLUMN MODIFICATION ********** DO 180 I = L,J P = X*H(I,K) + Y*H(I,K+1) IF ( .NOT. NOTLAS) GO TO 170 P = P + ZZ*H(I,K+2) H(I,K+2) = H(I,K+2) - P*R 170 H(I,K+1) = H(I,K+1) - P*Q H(I,K) = H(I,K) - P 180 CONTINUE C 190 CONTINUE C GO TO 40 C ********** ONE ROOT FOUND ********** 200 WR(EN) = X + T WI(EN) = 0.0 EN = NA GO TO 30 C ********** TWO ROOTS FOUND ********** 210 P = (Y-X)/2.0 Q = P*P + W ZZ = SQRT(ABS(Q)) X = X + T IF (Q.LT.0.0) GO TO 220 C ********** REAL PAIR ********** ZZ = P + SIGN(ZZ,P) WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ.NE.0.0) WR(EN) = X - W/ZZ WI(NA) = 0.0 WI(EN) = 0.0 GO TO 230 C ********** COMPLEX PAIR ********** 220 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 230 EN = ENM2 GO TO 30 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 240 IERR = EN 250 RETURN C ********** LAST CARD OF HQR ********** END SUBROUTINE CHECO(H,LDH,N,IPVT,RCOND,Z) INTEGER LDH,N,IPVT(N) REAL RCOND COMPLEX H(LDH,N),Z(N) C C *** PURPOSE: C C CHECO FACTORS A COMPLEX UPPER HESSENBERG MATRIX BY GAUSSIAN C ELIMINATION AND ESTIMATES THE CONDITION OF THE MATRIX. THE C LINPACK ANALOGUE OF CHECO IS CGECO. C C IF RCOND IS NOT NEEDED, CHEFA MAY BE SUBSTANTIALLY FASTER. TO C SOLVE THE LINEAR SYSTEM H*X = B (I.E., TO COMPUTE C (INVERSE(H)*B)), FOLLOW CHECO BY CHESL. C C ON ENTRY: C C H COMPLEX(LDH,N) C THE MATRIX TO BE FACTORED. C C LDH INTEGER C THE LEADING OR ROW DIMENSION OF THE ARRAY H AS DECLARED C IN THE MAIN CALLING PROGRAM. C C N INTEGER C THE ORDER OF THE MATRIX H. C C ON RETURN: C C H AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE C USED TO OBTAIN IT. THE FACTORIZATION CAN BE WRITTEN C H = L*U WHERE L IS A PRODUCT OF PERMUTATION AND UNIT C LOWER BIDIAGONAL MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND REAL C AN ESTIMATE OF THE RECIPROCAL OF THE CONDITION OF H. C FOR THE SYSTEM H*X = B , RELATIVE PERTURBATIONS C IN H AND B OF SIZE EPSILON MAY CAUSE RELATIVE PERTURBA- C TIONS IN X OF SIZE EPSILON/RCOND. IF RCOND IS SO SMALL C THAT THE LOGICAL EXPRESSION (1.0+RCOND) .EQ. 1.0 C IS TRUE, THEN H MAY BE SINGULAR TO WORKING PRECISION. C IN PARTICULAR, RCOND IS ZERO IF EXACT SINGULARITY IS C DETECTED OR THE ESTIMATE UNDERFLOWS. C C Z COMPLEX(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF H IS CLOSE TO A SINGULAR MATRIX, THEN Z IS AN C APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(H*Z) = RCOND*NORM(H)*NORM(Z) C C THIS VERSION DATED MAY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C SCL1,CHEFA,CL1NRM C REAL SCL1 C C INTERNAL VARIABLES: C INTEGER I,INFO,J,K,KB,KM1,KP1,L REAL HNORM,S,SM,YNORM,ZNORM COMPLEX EK,T,WK,WKM C C FORTRAN FUNCTIONS CALLED: C C COMPLEX CMPLX,CONJG C C INTERNAL FUNCTION: C COMPLEX Z1,Z2,CSIGN CSIGN(Z1,Z2) = SCL1(Z1)* (Z2/SCL1(Z2)) C IF (N.GT.1) GO TO 10 C C SPECIAL CASE OF SCALAR (N = 1) SYSTEM C HNORM = SCL1(H(1,1)) YNORM = HNORM GO TO 250 10 CONTINUE C C COMPUTE COMPLEX 1-NORM OF H C HNORM = 0.0E0 DO 30 J = 1,N S = 0.0E0 DO 20 I = 1,N S = S + SCL1(H(I,J)) 20 CONTINUE IF (S.GT.HNORM) HNORM = S 30 CONTINUE C C FACTOR C CALL CHEFA(H,LDH,N,IPVT,INFO) C C RCOND = 1/(NORM(H)*(ESTIMATE OF NORM(INVERSE(H))). C ESTIMATE = NORM(Z)/NORM(Y) WHERE H*Z =Y C AND TRANSPOSE(H)*Y = E. THE COMPONENTS OF E ARE C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS C OF W WHERE TRANSPOSE(U)*W = E. THE VECTORS ARE C FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE TRANSPOSE(U)*W = E C EK = (1.0E0,0.0E0) DO 40 J = 1,N Z(J) = (0.0E0,0.0E0) 40 CONTINUE DO 130 K = 1,N IF (SCL1(Z(K)).NE.0.0E0) EK = CSIGN(EK,-Z(K)) IF (SCL1(EK-Z(K)).LE.SCL1(H(K,K))) GO TO 60 S = SCL1(H(K,K))/SCL1(EK-Z(K)) DO 50 I = 1,N Z(I) = CMPLX(S,0.0E0)*Z(I) 50 CONTINUE EK = CMPLX(S,0.0E0)*EK 60 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = SCL1(WK) SM = SCL1(WKM) IF (SCL1(H(K,K)).EQ.0.0E0) GO TO 70 WK = WK/CONJG(H(K,K)) WKM = WKM/CONJG(H(K,K)) GO TO 80 70 CONTINUE WK = (1.0E0,0.0E0) WKM = (1.0E0,0.0E0) 80 CONTINUE KP1 = K + 1 IF (KP1.GT.N) GO TO 120 DO 90 J = KP1,N SM = SM + SCL1(Z(J)+WKM*CONJG(H(K,J))) Z(J) = Z(J) + WK*CONJG(H(K,J)) S = S + SCL1(Z(J)) 90 CONTINUE IF (S.GE.SM) GO TO 110 T = WKM - WK WK = WKM DO 100 J = KP1,N Z(J) = Z(J) + T*CONJG(H(K,J)) 100 CONTINUE 110 CONTINUE 120 CONTINUE Z(K) = WK 130 CONTINUE CALL CL1NRM(N,Z,ZNORM) C C SOLVE TRANSPOSE(L)*Y = W C DO 160 KB = 1,N K = N + 1 - KB IF (K.LT.N) Z(K) = Z(K) - CONJG(H(K+1,K))*Z(K+1) IF (SCL1(Z(K)).LE.1.0E0) GO TO 150 S = 1.0E0/SCL1(Z(K)) DO 140 I = 1,N Z(I) = CMPLX(S,0.0E0)*Z(I) 140 CONTINUE 150 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 160 CONTINUE CALL CL1NRM(N,Z,ZNORM) C YNORM = 1.0E0 C C SOLVE L*V = Y C DO 190 K = 1,N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K.LT.N) Z(K+1) = Z(K+1) - T*H(K+1,K) IF (SCL1(Z(K)).LE.1.0E0) GO TO 180 S = 1.0E0/SCL1(Z(K)) DO 170 I = 1,N Z(I) = CMPLX(S,0.0E0)*Z(I) 170 CONTINUE YNORM = S*YNORM 180 CONTINUE 190 CONTINUE CALL CL1NRM(N,Z,ZNORM) YNORM = YNORM/ZNORM C C SOLVE U*Z = V C DO 240 KB = 1,N K = N + 1 - KB IF (SCL1(Z(K)).LE.SCL1(H(K,K))) GO TO 210 S = SCL1(H(K,K))/SCL1(Z(K)) DO 200 I = 1,N Z(I) = CMPLX(S,0.0E0)*Z(I) 200 CONTINUE YNORM = S*YNORM 210 CONTINUE IF (SCL1(H(K,K)).NE.0.0E0) Z(K) = Z(K)/H(K,K) IF (SCL1(H(K,K)).EQ.0.0E0) Z(K) = (1.0E0,0.0E0) T = -Z(K) IF (K.LE.1) GO TO 230 KM1 = K - 1 DO 220 I = 1,KM1 Z(I) = Z(I) + T*H(I,K) 220 CONTINUE 230 CONTINUE 240 CONTINUE CALL CL1NRM(N,Z,ZNORM) YNORM = YNORM/ZNORM 250 CONTINUE IF (HNORM.NE.0.0E0) RCOND = YNORM/HNORM IF (HNORM.EQ.0.0E0) RCOND = 0.0E0 RETURN END SUBROUTINE CHEFA(H,LDH,N,IPVT,INFO) INTEGER LDH,N,IPVT(N),INFO COMPLEX H(LDH,N) C C *** PURPOSE: C C CHEFA FACTORS A COMPLEX UPPER HESSENBERG MATRIX BY GAUSSIAN C ELIMINATION. THE LINPACK ANALOGUE OF CHEFA IS CGEFA. C C CHEFA IS USUALLY CALLED BY CHECO, BUT IT CAN BE CALLED DIRECTLY C WITH A SUBSTANTIAL SAVING IN TIME IF RCOND IS NOT NEEDED. C C ON ENTRY: C C H COMPLEX(LDH,N) C THE MATRIX TO BE FACTORED. C C LDH INTEGER C THE LEADING OR ROW DIMENSION OF THE ARRAY H AS DECLARED C IN THE MAIN CALLING PROGRAM. C C N INTEGER C THE ORDER OF THE MATRIX H. C C ON RETURN: C C H AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH C WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN H = L*U WHERE L IS C A PRODUCT OF PERMUTATION AND UNIT LOWER BIDIAGONAL C MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE C = K IF U(K,K) .EQ. (0.0,0.0) THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES INDICATE C THAT CHESL WILL ATTEMPT TO DIVIDE BY ZERO IF CALLED. C USE RCOND IN CHECO FOR A RELIABLE INDICATION OF C SINGULARITY. C C THIS VERSION DATED MAY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C SCL1 C REAL SCL1 C C INTERNAL VARIABLES: C INTEGER J,K,KP1,L,NM1 COMPLEX T C C FORTRAN FUNCTIONS CALLED: C C NONE C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1.LT.1) GO TO 70 DO 60 K = 1,NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = K IF (SCL1(H(K+1,K)).GT.SCL1(H(K,K))) L = KP1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED C IF (SCL1(H(L,K)).EQ.0.0E0) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L.EQ.K) GO TO 10 T = H(L,K) H(L,K) = H(K,K) H(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS (NOTE: LINPACK COMPUTES NEGATIVE OF C WHAT IS COMPUTED HERE) C T = (1.0E0,0.0E0)/H(K,K) H(KP1,K) = T*H(KP1,K) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1,N T = H(L,J) IF (L.EQ.K) GO TO 20 H(L,J) = H(K,J) H(K,J) = T 20 CONTINUE H(KP1,J) = H(KP1,J) - T*H(KP1,K) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (SCL1(H(N,N)).EQ.0.0E0) INFO = N RETURN END SUBROUTINE CHESL(H,LDH,N,IPVT,B) INTEGER LDH,N,IPVT(N) COMPLEX H(LDH,N),B(N) C C *** PURPOSE: C C CHESL SOLVES THE COMPLEX UPPER HESSENBERG SYSTEM C H*X = B C USING THE FACTORS COMPUTED BY CHECO OR CHEFA. C THE LINPACK ANALOGUE OF CHESL IS CGESL. C C ON ENTRY: C C H COMPLEX(LDH,N) C THE OUTPUT FROM CHECO OR CHEFA. C C LDH INTEGER C THE LEADING OR ROW DIMENSION OF THE ARRAY H AS C DECLARED IN THE MAIN CALLING PROGRAM. C C N INTEGER C THE ORDER OF THE MATRIX H. C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM CHECO OR CHEFA. C C B COMPLEX(N) C THE RIGHT HAND SIDE VECTOR. C C ON RETURN: C C B THE SOLUTION VECTOR X (= (INVERSE(H)*B). C C ERROR CONDITION: C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY, THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDH. IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF CHECO HAS SET RCOND .GT. 0.0 C OR CHEFA HAS SET INFO .EQ. 0. C C TO COMPUTE INVERSE(H)*G WHERE G IS A MATRIX WITH M COLUMNS C CALL CHECO (H,LDH,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO ... C DO 10 J = 1,M C CALL CHESL (H,LDH,N,IPVT,G(1,J)) C 10 CONTINUE C C THIS VERSION DATED MAY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C NONE C C INTERNAL VARIABLES: C INTEGER I,K,KB,L,NM1 COMPLEX T C C FORTRAN FUNCTIONS CALLED: C C NONE C IF (N.EQ.1) GO TO 50 NM1 = N - 1 C C FIRST SOLVE L*Y = B C DO 20 K = 1,NM1 L = IPVT(K) T = B(L) IF (L.EQ.K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE B(K+1) = B(K+1) - T*H(K+1,K) 20 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1,NM1 K = N - KB B(K+1) = B(K+1)/H(K+1,K+1) T = -B(K+1) DO 30 I = 1,K B(I) = B(I) + T*H(I,K+1) 30 CONTINUE 40 CONTINUE 50 CONTINUE B(1) = B(1)/H(1,1) RETURN END SUBROUTINE CL1NRM(N,Z,ZNORM) INTEGER N REAL ZNORM COMPLEX Z(N) C C *** PURPOSE: C C CL1NRM COMPUTES AN L1-TYPE NORM OF A COMPLEX N-VECTOR Z AND C SCALES THE VECTOR Z BY THIS NORM. SPECIFICALLY, THE NORM OF Z C IS DETERMINED AS THE SUMMATION OF ABS(RE(Z(I)))+ABS(IM(Z(I))) C FOR I FROM 1 TO N WHERE RE(Z(I)) AND IM(Z(I)) ARE THE REAL AND C IMAGINARY PARTS, RESPECTIVELY, OF THE I-TH COMPONENT OF Z. C C ON ENTRY: C C N INTEGER C THE LENGTH OF THE VECTOR Z. C C Z COMPLEX(N) C THE VECTOR WHICH IS TO BE SCALED BY ITS NORM. C C ON RETURN: C C ZNORM REAL C THE NORM OF THE VECTOR Z. C C Z THE VECTOR Z SCALED BY ZNORM. C C THIS VERSION DATED JULY 1981. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C INTERNAL VARIABLES: C INTEGER I COMPLEX T C C FORTRAN FUNCTIONS CALLED: C REAL ABS,AIMAG,REAL C COMPLEX CMPLX C ZNORM = 0.0E0 DO 10 I = 1,N ZNORM = ZNORM + ABS(REAL(Z(I))) + ABS(AIMAG(Z(I))) 10 CONTINUE IF (ZNORM.EQ.0.0E0) RETURN T = CMPLX((1.0E0/ZNORM),0.0E0) DO 20 I = 1,N Z(I) = T*Z(I) 20 CONTINUE RETURN END REAL FUNCTION SCL1(Z) COMPLEX Z C C *** PURPOSE: C C SCL1 IS A REAL FUNCTION SUBPROGRAM WHICH RETURNS THE L1-TYPE C MAGNITUDE OF A COMPLEX NUMBER Z, I.E., THE NUMBER C ABS(RE(Z))+ABS(IM(Z)) WHERE RE(Z) AND IM(Z) ARE THE REAL AND C IMAGINARY PARTS, RESPECTIVELY, OF Z. C C THIS VERSION DATED MAY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C FORTRAN FUNCTIONS CALLED: C REAL ABS,AIMAG,REAL C SCL1 = ABS(REAL(Z)) + ABS(AIMAG(Z)) RETURN END SUBROUTINE SFRMG(FIRST,NA,NB,NC,L,M,N,A,B,C,FREQ,AINVB,G,H,IPVT, . ZC,ZR,EVRE,EVIM,IERR,INFO,RCOND) INTEGER NA,NB,NC,L,M,N,IPVT(N),IERR,INFO REAL A(NA,N),B(NB,M),C(NC,N),ZR(N),EVRE(N),EVIM(N),RCOND COMPLEX FREQ,AINVB(NB,M),G(NC,M),H(NA,N),ZC(N) LOGICAL FIRST C C *** PURPOSE: C C SFRMG TAKES REAL MATRICES A (N X N), B (N X M), AND C (L X N) C AND FORMS THE COMPLEX FREQUENCY RESPONSE MATRIX C G(FREQ) := C * (((FREQ * I) - A)-INVERSE) * B C WHERE I = (N X N) IDENTITY MATRIX AND FREQ IS A COMPLEX C SCALAR PARAMETER TAKING VALUES ALONG THE IMAGINARY AXIS FOR C CONTINUOUS-TIME SYSTEMS AND ON THE UNIT CIRCLE FOR DISCRETE- C TIME SYSTEMS. C C ON ENTRY: C C FIRST LOGICAL C SET = .TRUE. FOR THE FIRST CALL OF SFRMG WHEREUPON C IT IS SET TO .FALSE. FOR ALL SUBSEQUENT CALLS; MAY BE C SET = .FALSE. INITIALLY IF A IS ALREADY IN UPPER C HESSENBERG FORM AND NEITHER BALANCING NOR EIGENVALUES C OF A ARE DESIRED. C C NA INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY A C (AND THE COMPLEX ARRAY H) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C NB INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY B C (AND THE COMPLEX ARRAY AINVB) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C NC INTEGER C THE LEADING OR ROW DIMENSION OF THE REAL ARRAY C C (AND THE COMPLEX ARRAY G) AS DECLARED IN THE MAIN C CALLING PROGRAM. C C L INTEGER C THE NUMBER OF ROWS OF C (THE NUMBER OF OUTPUTS). C C M INTEGER C THE NUMBER OF COLUMNS OF B (THE NUMBER OF INPUTS). C C N INTEGER C THE ORDER OF THE MATRIX A (THE NUMBER OF STATES); C ALSO = NUMBER OF COLUMNS OF C = NUMBER OF ROWS OF B. C C A REAL(NA,N) C A REAL N X N MATRIX (THE SYSTEM MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C B REAL(NB,M) C A REAL N X M MATRIX (THE INPUT MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C C REAL(NC,N) C A REAL L X N MATRIX (THE OUTPUT MATRIX); NOT NEEDED AS C INPUT IF FIRST .EQ. .FALSE. C C FREQ COMPLEX C A COMPLEX SCALAR (THE FREQUENCY PARAMETER). C C IERR INTEGER C --FIRST = .FALSE. : -IERR IS IGNORED. C --FIRST = .TRUE. : -IF IERR .GE. 0, A WILL BE BALANCED C AND THE EIGENVALUES OF A WILL BE C COMPUTED (RECOMMENDED OPTION). C -IF IERR .LT. 0, A WILL NOT BE BAL- C ANCED AND THE EIGENVALUES OF A WILL C NOT BE COMPUTED. C C INFO INTEGER C --IF INFO .LT. 0, THE CONDITION OF (FREQ*I - A) WITH C RESPECT TO INVERSION WILL BE ESTIMATED (RECOMMENDED C OPTION). C --IF INFO .GE. 0, THE CONDITION OF (FREQ*I - A) WITH C RESPECT TO INVERSION WILL NOT BE ESTIMATED. C C ON RETURN: C C G COMPLEX(NC,M) C THE FREQUENCY RESPONSE MATRIX G(FREQ). C C A,B,C A IS IN UPPER HESSENBERG FORM WHILE B AND C HAVE BEEN C SUITABLY TRANSFORMED; IF FIRST .EQ. .FALSE. THESE C ARRAYS ARE NOT FURTHER MODIFIED. C C EVRE, REAL(N),REAL(N) C EVIM THE REAL AND IMAGINARY PARTS, RESPECTIVELY, OF THE C EIGENVALUES OF THE MATRIX A IF SFRMG WAS INITIALLY C CALLED WITH FIRST = .TRUE. .AND. IERR .GE. 0. C C IERR IF (FIRST = .TRUE. .AND. IERR .GE. 0 ON INPUT) THEN C IERR IS OUTPUT AS THE IERR PARAMETER OF THE EISPACK C SUBROUTINE HQR (HQR DOCUMENTATION MAY BE CONSULTED FOR C DETAILS); NORMAL RETURN IS 0. C C INFO --IGNORE IF INFO .LT. 0 ON INPUT. C --IF INFO .GE. 0 ON INPUT, THIS IS OUTPUT AS THE INFO C PARAMETER OF SUBROUTINE CHEFA (CHEFA MAY BE CONSULTED C FOR DETAILS); NORMAL RETURN IS 0. C C RCOND REAL C --IGNORE IF INFO .GE. 0 ON INPUT. C --IF INFO .LT. 0 ON INPUT, THIS IS OUTPUT AS THE RCOND C PARAMETER OF SUBROUTINE CHECO (CHECO MAY BE CONSULTED C FOR DETAILS); NORMAL RETURN IS WHEN C (1.0 + RCOND) .GT. 1.0. C C AINVB COMPLEX(NB,M) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C H COMPLEX(NA,N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C IPVT INTEGER(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C ZC COMPLEX(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C ZR REAL(N) C SCRATCH ARRAY USED FOR INTERNAL COMPUTATIONS. C C THIS VERSION DATED JUNE 1982. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C BALANC(EISPACK),CHECO,CHEFA,CHESL,HQR(EISPACK),SHETR C C INTERNAL VARIABLES: C INTEGER I,IGH,J,K,KK,KP,LOW REAL T C C FORTRAN FUNCTIONS CALLED: C REAL REAL C COMPLEX CMPLX C IF ( .NOT. FIRST) GO TO 150 IF (IERR.GE.0) GO TO 10 LOW = 1 IGH = N GO TO 90 10 CONTINUE CALL BALANC(NA,N,A,LOW,IGH,ZR) C C ADJUST B AND C MATRICES BASED ON INFORMATION IN THE VECTOR C ZR WHICH DESCRIBES THE BALANCING OF A AND IS DEFINED IN THE C SUBROUTINE BALANC C DO 40 K = 1,N KK = N - K + 1 IF (KK.GE.LOW .AND. KK.LE.IGH) GO TO 40 IF (KK.LT.LOW) KK = LOW - KK KP = ZR(KK) IF (KP.EQ.KK) GO TO 40 C C PERMUTE ROWS OF B C DO 20 J = 1,M T = B(KK,J) B(KK,J) = B(KP,J) B(KP,J) = T 20 CONTINUE C C PERMUTE COLUMNS OF C C DO 30 I = 1,L T = C(I,KK) C(I,KK) = C(I,KP) C(I,KP) = T 30 CONTINUE C 40 CONTINUE IF (IGH.EQ.LOW) GO TO 80 DO 70 K = LOW,IGH T = ZR(K) C C SCALE COLUMNS OF PERMUTED C C DO 50 I = 1,L C(I,K) = C(I,K)*T 50 CONTINUE C C SCALE ROWS OF PERMUTED B C DO 60 J = 1,M B(K,J) = B(K,J)/T 60 CONTINUE C 70 CONTINUE 80 CONTINUE 90 CONTINUE C C REDUCE A TO HESSENBERG FORM BY ORTHOGONAL SIMILARITIES AND C ACCUMULATE THE ORTHOGONAL TRANSFORMATIONS INTO B AND C C CALL SHETR(NA,NB,NC,L,M,N,LOW,IGH,A,B,C,ZR) IF (IERR.LT.0) GO TO 140 C C TEMPORARILY STORE HESSENBERG FORM OF A IN THE ARRAY H C DO 110 J = 1,N DO 100 I = 1,N H(I,J) = CMPLX(A(I,J),0.0E0) 100 CONTINUE 110 CONTINUE C C COMPUTE THE EIGENVALUES OF A IF THAT OPTION IS REQUESTED C CALL HQR(NA,N,LOW,IGH,A,EVRE,EVIM,IERR) C C RESTORE UPPER HESSENBERG FORM OF A TO THE ARRAY A C DO 130 J = 1,N DO 120 I = 1,N A(I,J) = REAL(H(I,J)) 120 CONTINUE 130 CONTINUE 140 CONTINUE FIRST = .FALSE. IF (IERR.GT.0) RETURN C C UPDATE H := (FREQ * I) - A WITH APPROPRIATE VALUE OF FREQ C 150 CONTINUE DO 170 J = 1,N DO 160 I = 1,N H(I,J) = -CMPLX(A(I,J),0.0E0) 160 CONTINUE H(J,J) = FREQ + H(J,J) 170 CONTINUE IF (INFO.LT.0) GO TO 180 C C FACTOR THE COMPLEX HESSENBERG MATRIX WITH NO CONDITION C ESTIMATE C CALL CHEFA(H,NA,N,IPVT,INFO) IF (INFO.EQ.0) GO TO 190 RETURN 180 CONTINUE C C FACTOR THE COMPLEX HESSENBERG MATRIX AND ESTIMATE ITS C CONDITION C CALL CHECO(H,NA,N,IPVT,RCOND,ZC) T = 1.0E0 + RCOND IF (T.EQ.1.0E0) RETURN 190 CONTINUE C C COMPUTE (H-INVERSE)*B C DO 220 J = 1,M DO 200 I = 1,N ZC(I) = CMPLX(B(I,J),0.0E0) 200 CONTINUE CALL CHESL(H,NA,N,IPVT,ZC) DO 210 I = 1,N AINVB(I,J) = ZC(I) 210 CONTINUE 220 CONTINUE C C COMPUTE C*(H-INVERSE)*B C DO 260 J = 1,M DO 230 I = 1,L G(I,J) = (0.0E0,0.0E0) 230 CONTINUE DO 250 K = 1,N DO 240 I = 1,L G(I,J) = G(I,J) + CMPLX(C(I,K),0.0E0)*AINVB(K,J) 240 CONTINUE C C G NOW CONTAINS THE DESIRED FREQUENCY RESPONSE MATRIX C 250 CONTINUE 260 CONTINUE RETURN END SUBROUTINE SHETR(NA,NB,NC,L,M,N,LOW,IGH,A,B,C,ORT) INTEGER NA,NB,NC,L,M,N,LOW,IGH REAL A(NA,N),B(NB,M),C(NC,N),ORT(N) C C *** PURPOSE: C C GIVEN A REAL GENERAL MATRIX A, SHETR REDUCES A SUBMATRIX C OF A IN ROWS AND COLUMNS LOW THROUGH IGH TO UPPER HESSENBERG C FORM BY ORTHOGONAL SIMILARITY TRANSFORMATIONS. THESE C ORTHOGONAL TRANSFORMATIONS ARE FURTHER ACCUMULATED INTO ROWS C LOW THROUGH IGH OF AN N X M MATRIX B AND COLUMNS LOW C THROUGH IGH OF AN L X N MATRIX C BY PREMULTIPLICATION AND C POSTMULTIPLICATION, RESPECTIVELY. C C THE CODE IS ADAPTED FROM THE EISPACK SUBROUTINE ORTHES C AND THE ORTHES DOCUMENTATION MAY BE CONSULTED FOR FURTHER C NUMERICAL DETAILS CONCERNING HOUSEHOLDER TRANSFORMATIONS. C C ON ENTRY: C C NA,NB,NC INTEGERS C THE LEADING OR ROW DIMENSIONS OF THE ARRAYS C A, B, AND C, RESPECTIVELY, AS DECLARED IN THE C MAIN CALLING PROGRAM. C C L INTEGER C NUMBER OF ROWS OF C. C C M INTEGER C NUMBER OF COLUMNS OF B. C C N INTEGER C ORDER OF A AND NUMBER OF ROWS OF B AND COLUMNS OF C. C C LOW,IGH INTEGERS C INTEGERS DETERMINED BY THE BALANCING SUBROUTINE C BALANC. IF BALANC HAS NOT BEEN USED, SET LOW = 1 C AND IGH = N. C C A REAL(NA,N) C AN N X N REAL GENERAL MATRIX TO BE REDUCED TO C UPPER HESSENBERG FORM. C C B REAL(NB,M) C AN N X M REAL MATRIX. C C C REAL(NC,N) C AN L X N REAL MATRIX. C C ON RETURN: C C A AN UPPER HESSENBERG MATRIX SIMILAR TO (VIA AN C ORTHOGONAL MATRIX CONSISTING OF A SEQUENCE OF C HOUSEHOLDER TRANSFORMATIONS) THE ORIGINAL MATRIX A; C FURTHER INFORMATION CONCERNING THE ORTHOGONAL C TRANSFORMATIONS USED IN THE REDUCTION IS CONTAINED C IN THE ELEMENTS BELOW THE FIRST SUBDIAGONAL; SEE C ORTHES DOCUMENTATION FOR DETAILS. C C B THE ORIGINAL B MATRIX PREMULTIPLIED BY THE TRANSPOSE C OF THE ORTHOGONAL TRANSFORMATION USED TO REDUCE A. C C C THE ORIGINAL C MATRIX POSTMULTIPLIED BY THE ORTHOGONAL C TRANSFORMATION USED TO REDUCE A. C C ORT REAL(N) C A WORK VECTOR CONTAINING INFORMATION ABOUT THE C ORTHOGONAL TRANSFORMATIONS; SEE ORTHES DOCUMENTATION C FOR DETAILS. C C THIS VERSION DATED JULY 1980. C ALAN J. LAUB, UNIVERSITY OF SOUTHERN CALIFORNIA. C C SUBROUTINES AND FUNCTIONS CALLED: C C NONE C C INTERNAL VARIABLES: C INTEGER I,II,J,JJ,K,KP1,KPN,LA REAL F,G,H,SCALE C C FORTRAN FUNCTIONS CALLED: C REAL ABS,SIGN,SQRT C LA = IGH - 1 KP1 = LOW + 1 IF (LA.LT.KP1) GO TO 170 DO 160 K = KP1,LA H = 0.0E0 ORT(K) = 0.0E0 SCALE = 0.0E0 C C SCALE COLUMN C DO 10 I = K,IGH SCALE = SCALE + ABS(A(I,K-1)) 10 CONTINUE IF (SCALE.EQ.0.0E0) GO TO 150 KPN = K + IGH DO 20 II = K,IGH I = KPN - II ORT(I) = A(I,K-1)/SCALE H = H + ORT(I)*ORT(I) 20 CONTINUE G = -SIGN(SQRT(H),ORT(K)) H = H - ORT(K)*G ORT(K) = ORT(K) - G C C FORM (I-(U*TRANSPOSE(U))/H)*A C DO 50 J = K,N F = 0.0E0 DO 30 II = K,IGH I = KPN - II F = F + ORT(I)*A(I,J) 30 CONTINUE F = F/H DO 40 I = K,IGH A(I,J) = A(I,J) - F*ORT(I) 40 CONTINUE 50 CONTINUE C C FORM (I-(U*TRANSPOSE(U))/H)*B C DO 80 J = 1,M F = 0.0E0 DO 60 II = K,IGH I = KPN - II F = F + ORT(I)*B(I,J) 60 CONTINUE F = F/H DO 70 I = K,IGH B(I,J) = B(I,J) - F*ORT(I) 70 CONTINUE 80 CONTINUE C C FORM (I-(U*TRANSPOSE(U))/H)*A*(I-(U*TRANSPOSE(U))/H) C DO 110 I = 1,IGH F = 0.0E0 DO 90 JJ = K,IGH J = KPN - JJ F = F + ORT(J)*A(I,J) 90 CONTINUE F = F/H DO 100 J = K,IGH A(I,J) = A(I,J) - F*ORT(J) 100 CONTINUE 110 CONTINUE C C FORM C*(I-(U*TRANSPOSE(U))/H) C DO 140 I = 1,L F = 0.0E0 DO 120 JJ = K,IGH J = KPN - JJ F = F + ORT(J)*C(I,J) 120 CONTINUE F = F/H DO 130 J = K,IGH C(I,J) = C(I,J) - F*ORT(J) 130 CONTINUE 140 CONTINUE ORT(K) = SCALE*ORT(K) A(K,K-1) = SCALE*G 150 CONTINUE 160 CONTINUE 170 CONTINUE RETURN END 5 3 1 -1 -1.0 -1.0 0.0 0.0 -2.0 -0.5 -2.0 -1.0 0.0 4.0 6.0 5.0 4.0 1.0 -3.0 3 1 -1 -1.25 -1.0 -1.25 -1.5 0.0 2.5 1.25 0.0 -0.75 0.75 -0.5 0.25 2.0 2.0 2.0 1 1 -1 -2.0 1.0 1.0 1 -1 1 -2.0 1.0 1.0 12 1 1 -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. -1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.