]>
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 LSQQR (A,X,B,M,N,IP,IM,IN,ERR,WORK) | |
11 | DIMENSION A(IM,N),X(IN,IP),B(IM,IP),WORK(1) | |
12 | LOGICAL ERR | |
13 | #if defined(CERNLIB_CDC) | |
14 | DATA ETAE/1.0E-14/ | |
15 | #endif | |
16 | #if !defined(CERNLIB_CDC) | |
17 | DATA ETAE/1.E-6/ | |
18 | #endif | |
19 | C | |
20 | C THE ARRAY A(M,N) CONTAINS THE GIVEN MATRIX OF AN | |
21 | C OVERDETERMINED SYSTEM OF M LINEAR EQUATIONS IN N | |
22 | C UNKNOWNS (M \ N). FOR THE IP RIGHT HAND SIDES GIVEN | |
23 | C AS THE COLUMNS OF THE ARRAY B(M,IP), THE LEAST SQUARES | |
24 | C SOLUTIONS ARE COMPUTED AND STORED AS THE COLUMNS OF THE | |
25 | C ARRAY X(N,IP). IF RANK(A) < N, THEN THE PROBLEM IS LEFT | |
26 | C UNSOLVED AND THE LOGICAL VARIABLE ERR IS SET TO .FALSE. | |
27 | C AND CONTROL IS RETURNED TO THE MAIN PROGRAM. | |
28 | C IN EITHER CASE, A AND B ARE LEFT INTACT. ETA IS THE | |
29 | C RELATIVE MACHINE PRECISION. | |
30 | C | |
31 | C THE ARRAY WORK (LENGTH 7*N+M*N+M) IS USED TO HOLD THE | |
32 | C FOLLOWING ARRAYS, WHICH ARE USED IN THIS AND THE CALLED | |
33 | C SUBROUTINES :- | |
34 | C | |
35 | C REAL ALPHA(N) WHERE WORK(1) = ALPHA(1). | |
36 | C REAL E(N) WHERE WORK(N+1) = E(1). | |
37 | C REAL Y(N) WHERE WORK(2*N+1) = Y(1). | |
38 | C REAL SUM(N) WHERE WORK(3*N+1) = SUM(1). | |
39 | C REAL Z(N) WHERE WORK(4*N+1) = Z(1). | |
40 | C REAL YA(N) WHERE WORK(5*N+1) = YA(1). | |
41 | C REAL R(M) WHERE WORK(6*N+1) = R(1). | |
42 | C REAL QR(M,N) WHERE WORK(6*N+M+1) = QR(1,1). | |
43 | C INTEGER IPIVOT(N) WHERE WORK(6*N+M*N+M+1) = IPIVOT(1). | |
44 | C | |
45 | C THUS, IN THIS SUBROUTINE, THESE ARRAYS ARE REFERENCED | |
46 | C BY THE ONE ARRAY NAME, WORK. | |
47 | C | |
48 | IMP1=M*N+6*N+M+1 | |
49 | IMP2=6*N+M+1 | |
50 | IMP3=6*N+1 | |
51 | IMP4=4*N+1 | |
52 | JO=IMP3-1 | |
53 | JO1=2*N | |
54 | C | |
55 | CALL PROC1(A,WORK(IMP2),M,N,IM) | |
56 | CALL DECOMP(M,N,WORK(IMP2),WORK(1),WORK(IMP1),ERR,WORK(5*N+1), | |
57 | 1 WORK(3*N+1)) | |
58 | IF(.NOT. ERR) RETURN | |
59 | C | |
60 | ETA=ETAE**2 | |
61 | DO 55 K=1,IP | |
62 | DO 10 I=1,M | |
63 | JU=JO+I | |
64 | 10 WORK(JU)=B(I,K) | |
65 | CALL SLV (M,N,WORK(IMP2),WORK(1),WORK(IMP1),WORK(IMP3), | |
66 | 1 WORK(2*N+1),WORK(IMP4)) | |
67 | DO 15 I=1,M | |
68 | JU=JO+I | |
69 | 15 WORK(JU)=-PROD2(A,WORK(2*N+1),IM,1,N,I,-B(I,K)) | |
70 | CALL SLV (M,N,WORK(IMP2),WORK(1),WORK(IMP1),WORK(IMP3), | |
71 | 1 WORK(N+1),WORK(IMP4)) | |
72 | YNORM0=0 | |
73 | ENORM1=0 | |
74 | DO 20 I=1,N | |
75 | JU=JO1+I | |
76 | YNORM0=YNORM0+WORK(JU)**2 | |
77 | JU=N+I | |
78 | 20 ENORM1=ENORM1+WORK(JU)**2 | |
79 | C | |
80 | C NO ATTEMPT AT OBTAINING THE SOLUTION IS MADE, UNLESS | |
81 | C THE NORM OF THE FIRST CORRECTION IS SIGNIFICANTLY SMALLER | |
82 | C THAN THE NORM OF THE INITIAL SOLUTION. | |
83 | IF(ENORM1.LE. 0.0625*YNORM0) GO TO 25 | |
84 | ERR=.FALSE. | |
85 | RETURN | |
86 | 25 DO 30 I=1,N | |
87 | JU=JO1+I | |
88 | JU1=N+I | |
89 | 30 WORK(JU)=WORK(JU)+WORK(JU1) | |
90 | C | |
91 | C TERMINATE THE ITERATION IF THE CORRECTION WAS OF LITTLE | |
92 | C SIGNIFICANCE. | |
93 | IF(ENORM1 .LT. ETA*YNORM0) GO TO 45 | |
94 | DO 35 I=1,M | |
95 | JU=JO+I | |
96 | 35 WORK(JU)=-PROD2(A,WORK(2*N+1),IM,1,N,I,-B(I,K)) | |
97 | CALL SLV (M,N,WORK(IMP2),WORK(1),WORK(IMP1),WORK(IMP3), | |
98 | 1 WORK(N+1),WORK(IMP4)) | |
99 | ENORM0=ENORM1 | |
100 | ENORM1=0 | |
101 | DO 40 I=1,N | |
102 | JU=N+I | |
103 | 40 ENORM1=ENORM1+WORK(JU)**2 | |
104 | C | |
105 | C TERMINATE THE ITERATION ALSO IF THE NORM OF THE | |
106 | C CORRECTION FAILED TO DECREASE SUFFICIENTLY, AS COMPARED | |
107 | C WITH THE NORM OF THE PREVIOUS CORRECTION. | |
108 | IF(ENORM1 .LE. 0.0625*ENORM0) GO TO 25 | |
109 | 45 DO 50 I=1,N | |
110 | JU=JO1+I | |
111 | 50 X(I,K)=WORK(JU) | |
112 | 55 CONTINUE | |
113 | RETURN | |
114 | END |