]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:27 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE DECOMP(M,N,QR,ALPHA,IPIVOT,ERR,YA,SUM) | |
11 | DIMENSION QR(M,N),ALPHA(N),IPIVOT(N),YA(N),SUM(N) | |
12 | LOGICAL ERR | |
13 | C | |
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. | |
29 | C | |
30 | DO 5 J=1,N | |
31 | C *** J TH. COLUMN SUM *** | |
32 | SUM(J)=PROD1(QR(1,J),QR(1,J),1,M) | |
33 | 5 IPIVOT(J)=J | |
34 | C | |
35 | DO 40 K=1,N | |
36 | C *** K TH. HOUSEHOLDER TRANSFORMATION *** | |
37 | SIGMA=SUM(K) | |
38 | JBAR=K | |
39 | L=K+1 | |
40 | IF(L.GT.N) GO TO 11 | |
41 | DO 10 J=L,N | |
42 | IF(SIGMA .GE. SUM(J)) GO TO 10 | |
43 | SIGMA=SUM(J) | |
44 | JBAR=J | |
45 | 10 CONTINUE | |
46 | 11 CONTINUE | |
47 | IF(JBAR .EQ. K) GO TO 20 | |
48 | C *** COLUMN INTERCHANGE *** | |
49 | I=IPIVOT(K) | |
50 | IPIVOT(K)=IPIVOT(JBAR) | |
51 | IPIVOT(JBAR)=I | |
52 | SUM(JBAR)=SUM(K) | |
53 | SUM(K)=SIGMA | |
54 | DO 15 I=1,M | |
55 | SIGMA=QR(I,K) | |
56 | QR(I,K)=QR(I,JBAR) | |
57 | 15 QR(I,JBAR)=SIGMA | |
58 | 20 CONTINUE | |
59 | SIGMA=PROD1(QR(1,K),QR(1,K),K,M) | |
60 | IF(SIGMA .NE. 0.0) GO TO 22 | |
61 | ERR=.FALSE. | |
62 | RETURN | |
63 | 22 QRKK=QR(K,K) | |
64 | ALPHA(K)=SQRT(SIGMA) | |
65 | IF(QRKK .GE. 0.0) ALPHA(K)=-ALPHA(K) | |
66 | ALPHAK=ALPHA(K) | |
67 | BETA=1.0/(SIGMA-QRKK*ALPHAK) | |
68 | QR(K,K)=QRKK-ALPHAK | |
69 | IF (L.GT.N) GO TO 40 | |
70 | DO 25 J=L,N | |
71 | 25 YA(J)=BETA*PROD1(QR(1,K),QR(1,J),K,M) | |
72 | DO 35 J=L,N | |
73 | DO 30 I=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 | |
76 | 40 CONTINUE | |
77 | RETURN | |
78 | END |