CTITLESACLEVN --MODIFIED LEVINSON RECURSION ROUTINE FOR COMPLEX MATRIX 00010006 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR C MOSER 00020006 CA DESIGNER C MOSER 00030006 CA LANGUAGE FORTRAN 77 00040006 CA SYSTEM IBM / CRAY 00050006 CA WRITTEN AUGUST 1988 00060006 C REVISED 01-15-89 JJC FOR SPARC PRODUCTION. 00070006 C REVISED 08-15-90 ESN IF R(1) = 0, SET H TO 0 AND RETURN 00080006 CA 00090006 CA 00100006 CA CALL SACLEVN (M,R,C,H,A,B) 00110006 CA 00120006 CA SUBROUTINES CALLED : NONE 00130006 CA 00140006 CA 00150006 C SOLVE THE MATRIX EQUATION R * H = C, WHERE R IS A COMPLEX 00160006 C TOEPPLITZ MATRIX (HERMITIAN, CONSTANT ALONG ALL DIAGONALS). 00170006 C USES A MODIFIED LEVINSON RECURSION 00180006 C 00190006 C ARGUMENTS: 00200006 C M INPUT ORDER OF EOUATION TO BE SOLVED (NUMBER OF 00210006 C DATA POINTS). 00220006 C R INPUT COMPLEX VECTOR OF LENGTH N CONTAINNG THE 00230006 C FIRST ROW OF THE MATRIX R. 00240006 C C INPUT COMPLEX VECTOR OF LENGTH N CONTAINING THE 00250006 C RIGHT HAND SIDE OF THE EQUATION.(CROSS- 00260006 C CORRELATION OF THE INPUT AND THE DESIRED 00270006 C OUTPUT FROM FILTER H). 00280006 C A INPUT COMPLEX WORK VECTOR OF LENGTH M. 00290006 C B INPUT COMPLEX WORK VECTOR OF LENGTH M 00300006 C H OUTPUT COMPLEX VECTOR OF LENGTH M CONTAINING THE 00310006 C SOLUTION. 00320006 C---------------------------------------------------------------------- 00330006 C 00340006 SUBROUTINE SACLEVN(M,R,C,H,A,B) 00350006 C 00360006 IMPLICIT INTEGER (A-Z) 00370006 COMPLEX R(M),C(M),H(M),A(M),B(M) 00380006 COMPLEX ALPHA,BETA,GAMMA,XK,Q 00390006 C CHECK FOR ZERO INPUT 00400006 IF (R(1) .EQ. 0.0) THEN 00410006 DO 5 I = 1, M 00420006 H(I) = 0.0 00430006 5 CONTINUE 00440006 RETURN 00450006 ENDIF 00460006 C ORDER 1 SOLUTION 00470006 A(1)=1. 00480006 BETA=R(2) 00490006 ALPHA=R(1) 00500006 H(1)=C(1)/R(1) 00510006 IF (M.EQ.1) RETURN 00520006 C ORDER 2 SOLUTION 00530006 GAMMA=H(1)*R(2) 00540006 XK=-BETA/ALPHA 00550006 A(2)=XK 00560006 ALPHA=ALPHA+XK*CONJG(BETA) 00570006 Q=(C(2)-GAMMA)/CONJG(ALPHA) 00580006 H(1)=H(1)+Q*CONJG(A(2)) 00590006 H(2)=Q 00600006 IF (M.EQ.2) RETURN 00610006 C RECURSION TAKES PLACE OVER LOOP 50 FOR ORDERS .GE. 3. 00620006 BETA=R(3)+A(2)*R(2) 00630006 GAMMA=H(1)*R(3)+H(2)*R(2) 00640006 M1=M-1 00650006 DO 50 N=2,M1 00660006 XK=-BETA/CONJG(ALPHA) 00670006 DO 10 I=2,N 00680006 B(I)=A(I)+XK*CONJG(A(N+2-I)) 00690006 10 CONTINUE 00700006 DO 20 I=2,N 00710006 A(I)=B(I) 00720006 20 CONTINUE 00730006 A(N+1)=XK 00740006 ALPHA=ALPHA+XK*CONJG(BETA) 00750006 Q=(C(N+1)-GAMMA)/CONJG(ALPHA) 00760006 H(1)=H(1)+Q*CONJG(A(N+1)) 00770006 DO 30 I=2,N 00780006 H(I)=H(I)+Q*CONJG(A(N+2-I)) 00790006 30 CONTINUE 00800006 H(N+1)=Q 00810006 IF (N.EQ.M1) RETURN 00820006 GAMMA=0. 00830006 BETA=0. 00840006 DO 40 I=1,N+1 00850006 BETA=BETA+A(I)*R(N-I+3) 00860006 GAMMA=GAMMA+H(I)*R(N-I+3) 00870006 40 CONTINUE 00880006 50 CONTINUE 00890006 RETURN 00900006 END 00910006