C C ________________________________________________________ C | | C | REDUCE A GENERAL MATRIX TO BIDIAGONAL FORM | C | | C | INPUT: | C | | C | D --ARRAY WITH AT LEAST N ELEMENTS | C | | C | B --ARRAY WITH AT LEAST M ELEMENTS | C | | C | A --ARRAY CONTAINING COEFFICIENT MATRIX | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | M --ROW DIMENSION OF MATRIX STORED IN A | C | | C | N --COLUMN DIMENSION OF MATRIX STORED IN A | C | | C | OUTPUT: | C | | C | D --DIAGONAL OF BIDIAGONAL FORM | C | | C | B --SUPERDIAGONAL (IF M .GE. N) OR SUBDIAG-| C | ONAL (IF N .GT. M) OF BIDIAGONAL FORM | C | | C | A --THE HOUSEHOLDER VECTORS USED IN THE | C | REDUCTION PROCESS | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C |________________________________________________________| C SUBROUTINE BIDAG(D,B,A,LA,M,N) INTEGER H,I,J,K,L,LA,M,N REAL A(LA,1),B(1),D(1),R,S,T,U K = 1 H = 2 IF ( M .LT. N ) GOTO 220 IF ( M .GT. 1 ) GOTO 10 D(1) = A(1,1) RETURN 10 J = K K = H DO 20 I = K,M 20 IF ( A(I,J) .NE. 0. ) GOTO 40 D(J) = A(J,J) A(J,J) = 0. DO 30 I = K,N 30 D(I) = A(J,I) GOTO 110 40 T = ABS(A(J,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 60 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 50 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 60 50 R = R + (S*U)**2 60 CONTINUE S = T*SQRT(R) R = A(J,J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(J) = -S A(J,J) = U*(R+S) DO 70 I = K,M 70 A(I,J) = A(I,J)*U IF ( K .GT. N ) RETURN DO 100 L = K,N T = 0. DO 80 I = J,M 80 T = T + A(I,J)*A(I,L) A(J,L) = A(J,L) - T*A(J,J) D(L) = A(J,L) DO 90 I = K,M 90 A(I,L) = A(I,L) - T*A(I,J) 100 CONTINUE 110 H = K + 1 IF ( K .LT. N ) GOTO 120 IF ( K .GT. N ) RETURN IF ( M .EQ. N ) GOTO 210 B(J) = A(J,N) GOTO 10 120 DO 130 I = H,N 130 IF ( D(I) .NE. 0. ) GOTO 140 B(J) = D(K) A(J,K) = 0. GOTO 10 140 T = ABS(D(K)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 160 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 150 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 160 150 R = R + (S*U)**2 160 CONTINUE S = T*SQRT(R) R = D(K) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(K) = U*(R+S) DO 170 I = H,N 170 D(I) = D(I)*U B(J) = -S DO 180 I = K,M 180 B(I) = 0. DO 190 L = K,N T = D(L) A(J,L) = T DO 190 I = K,M 190 B(I) = B(I) + T*A(I,L) DO 200 L = K,N T = D(L) DO 200 I = K,M 200 A(I,L) = A(I,L) - T*B(I) GOTO 10 210 D(N) = A(N,N) B(N-1) = A(N-1,N) RETURN 220 DO 230 I = K,N 230 D(I) = A(K,I) 240 J = K K = H DO 250 I = K,N 250 IF ( D(I) .NE. 0. ) GOTO 260 D(J) = A(J,J) A(J,J) = 0. GOTO 330 260 T = ABS(D(J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 280 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 270 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 280 270 R = R + (S*U)**2 280 CONTINUE S = T*SQRT(R) R = D(J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(J) = U*(R+S) DO 290 I = K,N 290 D(I) = D(I)*U U = -S IF ( K .GT. M ) GOTO 340 DO 300 I = K,M 300 B(I) = 0. DO 310 L = J,N T = D(L) A(J,L) = T DO 310 I = K,M 310 B(I) = B(I) + T*A(I,L) DO 320 L = J,N T = D(L) DO 320 I = K,M 320 A(I,L) = A(I,L) - T*B(I) 330 H = K + 1 D(J) = U IF ( K .LT. M ) GOTO 360 IF ( K .GT. M ) RETURN B(J) = A(M,J) GOTO 220 340 DO 350 I = J,N 350 A(J,I) = D(I) GOTO 330 360 DO 370 I = H,M 370 IF ( A(I,J) .NE. 0. ) GOTO 380 B(J) = A(K,J) A(K,J) = 0. GOTO 240 380 T = ABS(A(K,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 400 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 390 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 400 390 R = R + (S*U)**2 400 CONTINUE S = T*SQRT(R) R = A(K,J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S B(J) = -S A(K,J) = U*(R+S) DO 410 I = H,M 410 A(I,J) = A(I,J)*U DO 440 L = K,N T = 0. DO 420 I = K,M 420 T = T + A(I,J)*A(I,L) A(K,L) = A(K,L) - T*A(K,J) D(L) = A(K,L) DO 430 I = H,M 430 A(I,L) = A(I,L) - T*A(I,J) 440 CONTINUE GOTO 240 END