This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / e / lsqqr.F
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