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