5 * Revision 1.1.1.1 1996/04/01 15:02:27 mclareni
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)
13 #if defined(CERNLIB_CDC)
16 #if !defined(CERNLIB_CDC)
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.
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
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).
45 C THUS, IN THIS SUBROUTINE, THESE ARRAYS ARE REFERENCED
46 C BY THE ONE ARRAY NAME, WORK.
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),
65 CALL SLV (M,N,WORK(IMP2),WORK(1),WORK(IMP1),WORK(IMP3),
66 1 WORK(2*N+1),WORK(IMP4))
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))
76 YNORM0=YNORM0+WORK(JU)**2
78 20 ENORM1=ENORM1+WORK(JU)**2
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
89 30 WORK(JU)=WORK(JU)+WORK(JU1)
91 C TERMINATE THE ITERATION IF THE CORRECTION WAS OF LITTLE
93 IF(ENORM1 .LT. ETA*YNORM0) GO TO 45
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))
103 40 ENORM1=ENORM1+WORK(JU)**2
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