C C ________________________________________________________ C | | C | REDUCE A GENERAL MATRIX A TO BIDIAGONAL FORM | C | A = Q TIMES BIDIAGONAL MATRIX TIMES P TRANSPOSE | C | | C | INPUT: | C | | C | D --ARRAY WITH AT LEAST N ELEMENTS | C | | C | B --ARRAY WITH AT LEAST M ELEMENTS | C | | C | LQ --LEADING (ROW) DIMENSION OF ARRAY Q | C | | C | IQ --AN INTEGER WHICH INDICATES WHICH COL- | C | UMNS OF Q TO COMPUTE (= 0 MEANS NONE, | C | = 1 MEANS FIRST L, = 2 MEANS LAST M-L, | C | = 3 MEANS ALL M WHERE L = MIN(M,N)) | C | | C | LP --LEADING (ROW) DIMENSION OF ARRAY P | C | | C | IP --AN INTEGER (LIKE IQ) WHICH INDICATES | C | WHICH COLUMNS OF P TO COMPUTE | 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 M .LT. N) OF BIDIAGONAL FORM | C | | C | A --THE HOUSEHOLDER VECTORS USED IN THE | C | REDUCTION PROCESS | C | | C | Q --THE Q FACTOR IN THE BIDIAGONALIZATION | C | | C | P --THE P FACTOR IN THE BIDIAGONALIZATION | C | | C | NOTE: | C | | C | EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT | C | BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,| C | THEN THE HOUSEHOLDER VECTORS IN A ARE DESTROYED| C | | C | BUILTIN FUNCTIONS: ABS,MIN0,SQRT | C | PACKAGE SUBROUTINES: HSR1-HSR5 | C |________________________________________________________| C SUBROUTINE BIDAG2(D,B,Q,LQ,IQ,P,LP,IP,A,LA,M,N) INTEGER H,I,IP,IQ,J,JP,JQ,K,L,LA,LP,LQ,M,N REAL A(LA,1),B(1),D(1),Q(LQ,1),P(LP,1),R,S,T,U L = MIN0(M,N) IF ( IQ .GE. 0 ) GOTO 20 10 WRITE(6,*) 'ERROR: INPUT PARAMETER IQ FOR SUBROUTINE BIDAG2' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 20 IF ( IQ .GT. 3 ) GOTO 10 JQ = IQ IF ( IQ .LE. 1 ) GOTO 30 IF ( IQ .EQ. 3 ) GOTO 30 IF ( M .EQ. L ) JQ = 0 30 IF ( IP .GE. 0 ) GOTO 50 40 WRITE(6,*) 'ERROR: INPUT PARAMETER IP FOR SUBROUTINE BIDAG2' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 50 IF ( IP .GT. 3 ) GOTO 40 JP = IP IF ( IP .LE. 1 ) GOTO 60 IF ( IP .EQ. 3 ) GOTO 60 IF ( N .EQ. L ) JP = 0 60 K = 1 H = 2 IF ( M .LT. N ) GOTO 330 IF ( M .GT. 1 ) GOTO 70 D(1) = A(1,1) IF ( IQ .GT. 0 ) Q(1,1) = 1. IF ( IP .GT. 0 ) P(1,1) = 1. RETURN 70 J = K K = H DO 80 I = K,M 80 IF ( A(I,J) .NE. 0. ) GOTO 110 D(J) = A(J,J) A(J,J) = 0. DO 90 I = K,N 90 D(I) = A(J,I) IF ( JQ .EQ. 0 ) GOTO 200 DO 100 I = J,M 100 Q(I,J) = 0. GOTO 200 110 T = ABS(A(J,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 130 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 120 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 130 120 R = R + (S*U)**2 130 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 140 I = K,M 140 A(I,J) = A(I,J)*U IF ( JQ .EQ. 0 ) GOTO 160 DO 150 I = J,M 150 Q(I,J) = A(I,J) 160 IF ( K .GT. N ) GOTO 620 DO 190 L = K,N T = 0. DO 170 I = J,M 170 T = T + A(I,J)*A(I,L) A(J,L) = A(J,L) - T*A(J,J) D(L) = A(J,L) DO 180 I = K,M 180 A(I,L) = A(I,L) - T*A(I,J) 190 CONTINUE 200 H = K + 1 IF ( K .LT. N ) GOTO 210 IF ( K .GT. N ) GOTO 620 IF ( M .EQ. N ) GOTO 610 B(J) = A(J,N) GOTO 70 210 DO 220 I = H,N 220 IF ( D(I) .NE. 0. ) GOTO 240 B(J) = D(K) A(J,K) = 0. IF ( IP .EQ. 0 ) GOTO 70 DO 230 I = K,N 230 P(I,J) = 0. GOTO 70 240 T = ABS(D(K)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 260 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 250 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 260 250 R = R + (S*U)**2 260 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 270 I = H,N 270 D(I) = D(I)*U IF ( IP .EQ. 0 ) GOTO 290 DO 280 I = K,N 280 P(I,J) = D(I) 290 B(J) = -S DO 300 I = K,M 300 B(I) = 0. DO 310 L = K,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 = K,N T = D(L) DO 320 I = K,M 320 A(I,L) = A(I,L) - T*B(I) GOTO 70 330 DO 340 I = K,N 340 D(I) = A(K,I) 350 J = K K = H DO 360 I = K,N 360 IF ( D(I) .NE. 0. ) GOTO 370 U = D(J) D(J) = 0. GOTO 440 370 T = ABS(D(J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 390 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 380 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 390 380 R = R + (S*U)**2 390 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 400 I = K,N 400 D(I) = D(I)*U U = -S IF ( K .GT. M ) GOTO 470 DO 410 I = K,M 410 B(I) = 0. DO 420 L = J,N T = D(L) A(J,L) = T DO 420 I = K,M 420 B(I) = B(I) + T*A(I,L) DO 430 L = J,N T = D(L) DO 430 I = K,M 430 A(I,L) = A(I,L) - T*B(I) 440 H = K + 1 IF ( IP .EQ. 0 ) GOTO 460 DO 450 I = J,N 450 P(I,J) = D(I) 460 D(J) = U IF ( K .LT. M ) GOTO 490 IF ( K .GT. M ) GOTO 620 B(J) = A(M,J) GOTO 330 470 DO 480 I = J,N 480 A(J,I) = D(I) GOTO 440 490 DO 500 I = H,M 500 IF ( A(I,J) .NE. 0. ) GOTO 520 B(J) = A(K,J) A(K,J) = 0. IF ( IQ .EQ. 0 ) GOTO 330 DO 510 I = K,M 510 Q(I,J) = 0. GOTO 330 520 T = ABS(A(K,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 540 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 530 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 540 530 R = R + (S*U)**2 540 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 550 I = H,M 550 A(I,J) = A(I,J)*U IF ( IQ .EQ. 0 ) GOTO 570 DO 560 I = K,M 560 Q(I,J) = A(I,J) 570 DO 600 L = K,N T = 0. DO 580 I = K,M 580 T = T + A(I,J)*A(I,L) A(K,L) = A(K,L) - T*A(K,J) D(L) = A(K,L) DO 590 I = H,M 590 A(I,L) = A(I,L) - T*A(I,J) 600 CONTINUE GOTO 350 610 D(N) = A(N,N) B(N-1) = A(N-1,N) 620 IF ( JQ .EQ. 0 ) GOTO 650 IF ( N .GT. M ) GOTO 640 IF ( N .EQ. M ) GOTO 630 IF ( JQ .EQ. 1 ) CALL HSR3(Q,LQ,M,N) IF ( JQ .EQ. 2 ) CALL HSR4(Q,LQ,M,N) IF ( JQ .EQ. 3 ) CALL HSR5(Q,LQ,M,N) GOTO 650 630 CALL HSR2(Q,LQ,M) GOTO 650 640 CALL HSR1(Q,LQ,M) 650 IF ( JP .EQ. 0 ) RETURN IF ( N .LE. M ) GOTO 660 IF ( JP .EQ. 1 ) CALL HSR3(P,LP,N,M) IF ( JP .EQ. 2 ) CALL HSR4(P,LP,N,M) IF ( JP .EQ. 3 ) CALL HSR5(P,LP,N,M) RETURN 660 CALL HSR1(P,LP,N) RETURN END