relibc_openlibm/slatec/du11ls.f

*DECK DU11LS
      SUBROUTINE DU11LS (A, MDA, M, N, UB, DB, MODE, NP, KRANK, KSURE,
     +   H, W, EB, IC, IR)
C***BEGIN PROLOGUE  DU11LS
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DLLSIA
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (U11LS-S, DU11LS-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C   **** Double Precision version of U11LS ****
C
C       This routine performs a QR factorization of A
C       using Householder transformations. Row and
C       column pivots are chosen to reduce the growth
C       of round-off and to help detect possible rank
C       deficiency.
C
C***SEE ALSO  DLLSIA
C***ROUTINES CALLED  DAXPY, DDOT, DNRM2, DSCAL, DSWAP, IDAMAX, ISWAP,
C                    XERMSG
C***REVISION HISTORY  (YYMMDD)
C   810801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   891009  Removed unreferenced variable.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900328  Added TYPE section.  (WRB)
C***END PROLOGUE  DU11LS
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION DDOT,DNRM2
      DIMENSION A(MDA,*),UB(*),DB(*),H(*),W(*),EB(*)
      INTEGER IC(*),IR(*)
C
C        INITIALIZATION
C
C***FIRST EXECUTABLE STATEMENT  DU11LS
      J=0
      KRANK=N
      DO 10 I=1,N
      IC(I)=I
   10 CONTINUE
      DO 12 I=1,M
      IR(I)=I
   12 CONTINUE
C
C        DETERMINE REL AND ABS ERROR VECTORS
C
C
C
C        CALCULATE COL LENGTH
C
      DO 30 I=1,N
      H(I)=DNRM2(M,A(1,I),1)
      W(I)=H(I)
   30 CONTINUE
C
C         INITIALIZE ERROR BOUNDS
C
      DO  40 I=1,N
      EB(I)=MAX(DB(I),UB(I)*H(I))
      UB(I)=EB(I)
      DB(I)=0.0D0
   40 CONTINUE
C
C          DISCARD SELF DEPENDENT COLUMNS
C
      I=1
   50 IF(EB(I).GE.H(I)) GO TO 60
      IF(I.EQ.KRANK) GO TO 70
      I=I+1
      GO TO 50
C
C          MATRIX REDUCTION
C
   60 CONTINUE
      KK=KRANK
      KRANK=KRANK-1
      IF(MODE.EQ.0) RETURN
      IF(I.GT.NP) GO TO  64
      CALL XERMSG ('SLATEC', 'DU11LS',
     +   'FIRST NP COLUMNS ARE LINEARLY DEPENDENT', 8, 0)
      KRANK=I-1
      RETURN
   64 CONTINUE
      IF(I.GT.KRANK) GO TO 70
      CALL DSWAP(1,EB(I),1,EB(KK),1)
      CALL DSWAP(1,UB(I),1,UB(KK),1)
      CALL DSWAP(1,W(I),1,W(KK),1)
      CALL DSWAP(1,H(I),1,H(KK),1)
      CALL ISWAP(1,IC(I),1,IC(KK),1)
      CALL DSWAP(M,A(1,I),1,A(1,KK),1)
      GO TO 50
C
C           TEST FOR ZERO RANK
C
   70 IF(KRANK.GT.0) GO TO 80
      KRANK=0
      KSURE=0
      RETURN
   80 CONTINUE
C
C        M A I N    L O O P
C
  110 CONTINUE
      J=J+1
      JP1=J+1
      JM1=J-1
      KZ=KRANK
      IF(J.LE.NP) KZ=J
C
C        EACH COL HAS MM=M-J+1 COMPONENTS
C
      MM=M-J+1
C
C         UB DETERMINES COLUMN PIVOT
C
  115 IMIN=J
      IF(H(J).EQ.0.D0) GO TO 170
      RMIN=UB(J)/H(J)
      DO 120 I=J,KZ
      IF(UB(I).GE.H(I)*RMIN) GO TO 120
      RMIN=UB(I)/H(I)
      IMIN=I
  120 CONTINUE
C
C     TEST FOR RANK DEFICIENCY
C
      IF(RMIN.LT.1.0D0) GO TO 200
      TT=(EB(IMIN)+ABS(DB(IMIN)))/H(IMIN)
      IF(TT.GE.1.0D0) GO TO 170
C     COMPUTE EXACT UB
      DO 125 I=1,JM1
      W(I)=A(I,IMIN)
  125 CONTINUE
      L=JM1
  130 W(L)=W(L)/A(L,L)
      IF(L.EQ.1) GO TO 150
      LM1=L-1
      DO 140 I=L,JM1
      W(LM1)=W(LM1)-A(LM1,I)*W(I)
  140 CONTINUE
      L=LM1
      GO TO 130
  150 TT=EB(IMIN)
      DO 160 I=1,JM1
      TT=TT+ABS(W(I))*EB(I)
  160 CONTINUE
      UB(IMIN)=TT
      IF(UB(IMIN)/H(IMIN).GE.1.0D0) GO TO 170
      GO TO 200
C
C        MATRIX REDUCTION
C
  170 CONTINUE
      KK=KRANK
      KRANK=KRANK-1
      KZ=KRANK
      IF(MODE.EQ.0) RETURN
      IF(J.GT.NP) GO TO 172
      CALL XERMSG ('SLATEC', 'DU11LS',
     +   'FIRST NP COLUMNS ARE LINEARLY DEPENDENT', 8, 0)
      KRANK=J-1
      RETURN
  172 CONTINUE
      IF(IMIN.GT.KRANK) GO TO 180
      CALL ISWAP(1,IC(IMIN),1,IC(KK),1)
      CALL DSWAP(M,A(1,IMIN),1,A(1,KK),1)
      CALL DSWAP(1,EB(IMIN),1,EB(KK),1)
      CALL DSWAP(1,UB(IMIN),1,UB(KK),1)
      CALL DSWAP(1,DB(IMIN),1,DB(KK),1)
      CALL DSWAP(1,W(IMIN),1,W(KK),1)
      CALL DSWAP(1,H(IMIN),1,H(KK),1)
  180 IF(J.GT.KRANK) GO TO 300
      GO TO 115
C
C        COLUMN PIVOT
C
  200 IF(IMIN.EQ.J) GO TO 230
      CALL DSWAP(1,H(J),1,H(IMIN),1)
      CALL DSWAP(M,A(1,J),1,A(1,IMIN),1)
      CALL DSWAP(1,EB(J),1,EB(IMIN),1)
      CALL DSWAP(1,UB(J),1,UB(IMIN),1)
      CALL DSWAP(1,DB(J),1,DB(IMIN),1)
      CALL DSWAP(1,W(J),1,W(IMIN),1)
      CALL ISWAP(1,IC(J),1,IC(IMIN),1)
C
C        ROW PIVOT
C
  230 CONTINUE
      JMAX=IDAMAX(MM,A(J,J),1)
      JMAX=JMAX+J-1
      IF(JMAX.EQ.J) GO TO 240
      CALL DSWAP(N,A(J,1),MDA,A(JMAX,1),MDA)
      CALL ISWAP(1,IR(J),1,IR(JMAX),1)
  240 CONTINUE
C
C     APPLY HOUSEHOLDER TRANSFORMATION
C
      TN=DNRM2(MM,A(J,J),1)
      IF(TN.EQ.0.0D0) GO TO 170
      IF(A(J,J).NE.0.0D0) TN=SIGN(TN,A(J,J))
      CALL DSCAL(MM,1.0D0/TN,A(J,J),1)
      A(J,J)=A(J,J)+1.0D0
      IF(J.EQ.N) GO TO 250
      DO 248 I=JP1,N
      BB=-DDOT(MM,A(J,J),1,A(J,I),1)/A(J,J)
      CALL DAXPY(MM,BB,A(J,J),1,A(J,I),1)
      IF(I.LE.NP) GO TO 248
      IF(H(I).EQ.0.0D0) GO TO 248
      TT=1.0D0-(ABS(A(J,I))/H(I))**2
      TT=MAX(TT,0.0D0)
      T=TT
      TT=1.0D0+.05D0*TT*(H(I)/W(I))**2
      IF(TT.EQ.1.0D0) GO TO 244
      H(I)=H(I)*SQRT(T)
      GO TO 246
  244 CONTINUE
      H(I)=DNRM2(M-J,A(J+1,I),1)
      W(I)=H(I)
  246 CONTINUE
  248 CONTINUE
  250 CONTINUE
      H(J)=A(J,J)
      A(J,J)=-TN
C
C
C          UPDATE UB, DB
C
      UB(J)=UB(J)/ABS(A(J,J))
      DB(J)=(SIGN(EB(J),DB(J))+DB(J))/A(J,J)
      IF(J.EQ.KRANK) GO TO 300
      DO 260 I=JP1,KRANK
      UB(I)=UB(I)+ABS(A(J,I))*UB(J)
      DB(I)=DB(I)-A(J,I)*DB(J)
  260 CONTINUE
      GO TO 110
C
C        E N D    M A I N    L O O P
C
  300 CONTINUE
C
C        COMPUTE KSURE
C
      KM1=KRANK-1
      DO 318 I=1,KM1
      IS=0
      KMI=KRANK-I
      DO 315 II=1,KMI
      IF(UB(II).LE.UB(II+1)) GO TO 315
      IS=1
      TEMP=UB(II)
      UB(II)=UB(II+1)
      UB(II+1)=TEMP
  315 CONTINUE
      IF(IS.EQ.0) GO TO 320
  318 CONTINUE
  320 CONTINUE
      KSURE=0
      SUM=0.0D0
      DO 328 I=1,KRANK
      R2=UB(I)*UB(I)
      IF(R2+SUM.GE.1.0D0) GO TO 330
      SUM=SUM+R2
      KSURE=KSURE+1
  328 CONTINUE
  330 CONTINUE
C
C     IF SYSTEM IS OF REDUCED RANK AND MODE = 2
C     COMPLETE THE DECOMPOSITION FOR SHORTEST LEAST SQUARES SOLUTION
C
      IF(KRANK.EQ.N .OR. MODE.LT.2) GO TO 360
      NMK=N-KRANK
      KP1=KRANK+1
      I=KRANK
  340 TN=DNRM2(NMK,A(I,KP1),MDA)/A(I,I)
      TN=A(I,I)*SQRT(1.0D0+TN*TN)
      CALL DSCAL(NMK,1.0D0/TN,A(I,KP1),MDA)
      W(I)=A(I,I)/TN+1.0D0
      A(I,I)=-TN
      IF(I.EQ.1) GO TO 350
      IM1=I-1
      DO 345 II=1,IM1
      TT=-DDOT(NMK,A(II,KP1),MDA,A(I,KP1),MDA)/W(I)
      TT=TT-A(II,I)
      CALL DAXPY(NMK,TT,A(I,KP1),MDA,A(II,KP1),MDA)
      A(II,I)=A(II,I)+TT*W(I)
  345 CONTINUE
      I=I-1
      GO TO 340
  350 CONTINUE
  360 CONTINUE
      RETURN
      END