5 * Revision 1.1.1.1 1996/04/01 15:02:27 mclareni
10 SUBROUTINE DECOMP(M,N,QR,ALPHA,IPIVOT,ERR,YA,SUM)
11 DIMENSION QR(M,N),ALPHA(N),IPIVOT(N),YA(N),SUM(N)
14 C THIS ROUTINE REDUCES THE MATRIX GIVEN IN THE ARRAY
15 C QR(M,N), WHERE M \ N, TO UPPER RIGHT TRIANGULAR FORM
16 C BY MEANS OF N ELEMENTARY ORTHOGONAL TRANSFORMATIONS
17 C ( I - BETA UT , WHERE T IS THE TRANSPOSE OF U).
18 C THE DIAGONAL ELEMENTS OF THE REDUCED MATRIX ARE STORED
19 C IN THE ARRAY ALPHA(N), THE OFF-DIAGONAL ELEMENTS IN THE
20 C UPPER RIGHT TRIANGULAR PART OF QR. THE NON-ZERO
21 C COMPONENTS OF THE VECTORS U ARE STORED ON AND BELOW
22 C THE LEADING DIAGONAL OF QR. PIVOTING IS DONE BY
23 C CHOOSING AT EACH STEP, THE COLUMN WITH THE LARGEST SUM
24 C OF SQUARES TO BE REDUCED NEXT. THESE INTERCHANGES ARE
25 C RECORDED IN THE ARRAY IPIVOT(N). IF AT ANY STAGE, THE
26 C SUM OF SQUARES OF THE COLUMN TO BE REDUCED IS EXACTLY
27 C EQUAL TO ZERO, THEN THE LOGICAL VARIABLE ERR IS SET TO
28 C .FALSE. AND CONTROL IS RETURNED TO THE MAIN PROGRAM.
31 C *** J TH. COLUMN SUM ***
32 SUM(J)=PROD1(QR(1,J),QR(1,J),1,M)
36 C *** K TH. HOUSEHOLDER TRANSFORMATION ***
42 IF(SIGMA .GE. SUM(J)) GO TO 10
47 IF(JBAR .EQ. K) GO TO 20
48 C *** COLUMN INTERCHANGE ***
50 IPIVOT(K)=IPIVOT(JBAR)
59 SIGMA=PROD1(QR(1,K),QR(1,K),K,M)
60 IF(SIGMA .NE. 0.0) GO TO 22
65 IF(QRKK .GE. 0.0) ALPHA(K)=-ALPHA(K)
67 BETA=1.0/(SIGMA-QRKK*ALPHAK)
71 25 YA(J)=BETA*PROD1(QR(1,K),QR(1,J),K,M)
74 30 QR(I,J)=QR(I,J)-QR(I,K)*YA(J)
75 35 SUM(J)=SUM(J)-QR(K,J)**2