C     documentation: see constr.txt

C     modification by w.h. 2/98:
      SUBROUTINE GETPR(DRELPR)
C     COMPUTE MACHINE PRECISION=DRELPR
      DOUBLE PRECISION DRELPR, SUM, ONE, HALF
      DATA HALF /0.5D0/, ONE /1.D0/
      DRELPR = ONE
   50 SUM=ONE+DRELPR
C     in old version the latter calculation was inside the following condition:
C  50 IF (ONE+DRELPR.EQ.ONE) GO TO 60
C     this had the effect that the calculated precision was that of the 
C     processor which may be bigger than that of the memory storage
      IF (SUM.EQ.ONE) GO TO 60
      DRELPR = DRELPR*HALF
      GO TO 50
   60 DRELPR = DRELPR + DRELPR
      RETURN
C     DRELPR=2.22045e-16 on IBM and Origin
      END

C     ALGORITHM 587, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL. 8, NO. 3,
C     SEP., 1982, P.323.
      SUBROUTINE LSEI(W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, RNORML, LSEI  10
     * MODE, WS, IP)
      DOUBLE PRECISION W(MDW,1), PRGOPT(1), X(1), WS(1), RNORME, RNORML
      INTEGER IP(3)
      DOUBLE PRECISION DUMMY, ENORM, DRELPR, FNORM, GAM, HALF, ONE, RB,
     *RN, RNMAX, SIZE, SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME, ZERO
      DOUBLE PRECISION DASUM, DDOT, DNRM2, DSQRT, DABS, DMAX1
      LOGICAL COV
      DATA ZERO /0.D0/, DRELPR /0.D0/, HALF /0.5D0/, ONE /1.D0/
C
C     CHECK THAT ENOUGH STORAGE WAS ALLOCATED IN WS(*) AND IP(*).
      IF (.NOT.(IP(1).GT.0)) GO TO 20
      LCHK = 2*(ME+N) + MAX0(MA+MG,N) + (MG+2)*(N+7)
      IF (.NOT.(IP(1).LT.LCHK)) GO TO 10
      MODE = 4
      NERR = 2
      IOPT = 1

      CALL XERRWV(67HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR WS(*),
     *NEED LW=I1 BELOW, 67, NERR, IOPT, 1, LCHK, 0,
     * 0, DUMMY, DUMMY)
      RETURN
   10 CONTINUE
   20 IF (.NOT.(IP(2).GT.0)) GO TO 40
      LCHK = MG + 2*N + 2
      IF (.NOT.(IP(2).LT.LCHK)) GO TO 30
      MODE = 4
      NERR = 2
      IOPT = 1
      CALL XERRWV(68HLSEI( ), INSUFFICIENT STORAGE ALLOCATED FOR IP(*),
     *NEED LIP=I1 BELOW, 68, NERR, IOPT, 1, LCHK, 0,
     * 0, DUMMY, DUMMY)
      RETURN
   30 CONTINUE
C
C     COMPUTE MACHINE PRECISION=DRELPR ONLY WHEN NECESSARY.
   40 IF (.NOT.(DRELPR.EQ.ZERO)) GO TO 70
C     modification by w.h. 2/98, see hist9 for old version
      CALL GETPR(DRELPR)
C
C     COMPUTE NUMBER OF POSSIBLE RIGHT MULTIPLYING HOUSEHOLDER
C     TRANSFORMATIONS.
   70 M = ME + MA + MG
      MODE = 0
      IF (N.LE.0 .OR. M.LE.0) RETURN
      IF (.NOT.(MDW.LT.M)) GO TO 80
      NERR = 1
      IOPT = 1
      CALL XERROR(36HLSEI( ), MDW.LT.ME+MA+MG IS AN ERROR, 36, NERR,
     * IOPT)
      MODE = 4
      RETURN
   80 NP1 = N + 1
      KRANKE = MIN0(ME,N)
      N1 = 2*KRANKE + 1
      N2 = N1 + N
C
C     PROCESS-OPTION-VECTOR
      ASSIGN 90 TO IGO990
      GO TO 480
   90 IF (.NOT.(COV .AND. MDW.LT.N)) GO TO 100
      NERR = 2
      IOPT = 1
      CALL XERROR(
     * 54HLSEI( ), MDW.LT.N, WHEN COV MATRIX NEEDED, IS AN ERROR, 54,
     * NERR, IOPT)
      MODE = 4
      RETURN
  100 L = KRANKE
C
C     COMPUTE NORM OF EQUALITY CONSTRAINT MATRIX AND RT SIDE.
      ENORM = ZERO
      DO 110 J=1,N
        ENORM = DMAX1(ENORM,DASUM(ME,W(1,J),1))
  110 CONTINUE
      FNORM = DASUM(ME,W(1,NP1),1)
      IF (.NOT.(L.GT.0)) GO TO 190
      SNMAX = ZERO
      RNMAX = ZERO
      DO 180 I=1,L
C
C     COMPUTE MAXIMUM RATIO OF VECTOR LENGTHS. PARTITION
C     IS AT COL. I.
        DO 150 K=I,ME
          SN = DDOT(N-I+1,W(K,I),MDW,W(K,I),MDW)
          RN = DDOT(I-1,W(K,1),MDW,W(K,1),MDW)
          IF (.NOT.(RN.EQ.ZERO .AND. SN.GT.SNMAX)) GO TO 120
          SNMAX = SN
          IMAX = K
          GO TO 140
  120     IF (.NOT.(K.EQ.I .OR. (SN*RNMAX.GT.RN*SNMAX))) GO TO 130
          SNMAX = SN
          RNMAX = RN
          IMAX = K
  130     CONTINUE
  140     CONTINUE
  150   CONTINUE
C
C     INTERCHANGE ROWS IF NECESSARY.
        IF (I.NE.IMAX) CALL DSWAP(NP1, W(I,1), MDW, W(IMAX,1), MDW)
        IF (.NOT.(SNMAX.GT.TAU**2*RNMAX)) GO TO 160
C
C     ELIMINATE ELEMS I+1,...,N IN ROW I.
        CALL H12(1, I, I+1, N, W(I,1), MDW, WS(I), W(I+1,1), MDW, 1,
     *   M-I)
        GO TO 170
  160   KRANKE = I - 1
        GO TO 200
  170   CONTINUE
  180 CONTINUE
  190 CONTINUE
  200 CONTINUE
C
C     SAVE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL DCOPY(KRANKE, W, MDW+1, WS(KRANKE+1), 1)
C
C     USE HOUSEHOLDER TRANS FROM LEFT TO ACHIEVE KRANKE BY KRANKE UPPER
C     TRIANGULAR FORM.
      IF (.NOT.(KRANKE.GT.0 .AND. KRANKE.LT.ME)) GO TO 220
      DO 210 KK=1,KRANKE
        K = KRANKE + 1 - KK
C
C     APPLY TRANFORMATION TO MATRIX COLS. 1,...,K-1.
        CALL H12(1, K, KRANKE+1, ME, W(1,K), 1, UP, W, 1, MDW, K-1)
C
C     APPLY TO RT SIDE VECTOR.
        CALL H12(2, K, KRANKE+1, ME, W(1,K), 1, UP, W(1,NP1), 1, 1, 1)
  210 CONTINUE
  220 IF (.NOT.(KRANKE.GT.0)) GO TO 240
C
C     SOLVE FOR VARIABLES 1,...,KRANKE IN NEW COORDINATES.
      CALL DCOPY(KRANKE, W(1,NP1), 1, X, 1)
      DO 230 I=1,KRANKE
        X(I) = (X(I)-DDOT(I-1,W(I,1),MDW,X,1))/W(I,I)
  230 CONTINUE
C
C     COMPUTE RESIDUALS FOR REDUCED PROBLEM.
  240 MEP1 = ME + 1
      RNORML = ZERO
      IF (.NOT.(ME.LT.M)) GO TO 270
      DO 260 I=MEP1,M
        W(I,NP1) = W(I,NP1) - DDOT(KRANKE,W(I,1),MDW,X,1)
        SN = DDOT(KRANKE,W(I,1),MDW,W(I,1),MDW)
        RN = DDOT(N-KRANKE,W(I,KRANKE+1),MDW,W(I,KRANKE+1),MDW)
        IF (.NOT.(RN.LE.TAU**2*SN .AND. KRANKE.LT.N)) GO TO 250
        W(I,KRANKE+1) = ZERO
        CALL DCOPY(N-KRANKE, W(I,KRANKE+1), 0, W(I,KRANKE+1), MDW)
  250   CONTINUE
  260 CONTINUE
C
C     COMPUTE EQUAL. CONSTRAINT EQUAS. RESIDUAL LENGTH.
  270 RNORME = DNRM2(ME-KRANKE,W(KRANKE+1,NP1),1)
C
C     MOVE REDUCED PROBLEM DATA UPWARD IF KRANKE.LT.ME.
      IF (.NOT.(KRANKE.LT.ME)) GO TO 290
      DO 280 J=1,NP1
        CALL DCOPY(M-ME, W(ME+1,J), 1, W(KRANKE+1,J), 1)
  280 CONTINUE
C
C     COMPUTE SOLN OF REDUCED PROBLEM.
  290 CALL LSI(W(KRANKE+1,KRANKE+1), MDW, MA, MG, N-KRANKE, PRGOPT,
     * X(KRANKE+1), RNORML, MODE, WS(N2), IP(2))
      IF (.NOT.(ME.GT.0)) GO TO 330
C
C     TEST FOR CONSISTENCY OF EQUALITY CONSTRAINTS.
      MDEQC = 0
      XNRME = DASUM(KRANKE,W(1,NP1),1)
      IF (RNORME.GT.TAU*(ENORM*XNRME+FNORM)) MDEQC = 1
      MODE = MODE + MDEQC
C
C     CHECK IF SOLN TO EQUAL. CONSTRAINTS SATISFIES INEQUAL.
C     CONSTRAINTS WHEN THERE ARE NO DEGREES OF FREEDOM LEFT.
      IF (.NOT.(KRANKE.EQ.N .AND. MG.GT.0)) GO TO 320
      XNORM = DASUM(N,X,1)
      MAPKE1 = MA + KRANKE + 1
      MEND = MA + KRANKE + MG
      DO 310 I=MAPKE1,MEND
        SIZE = DASUM(N,W(I,1),MDW)*XNORM + DABS(W(I,NP1))
        IF (.NOT.(W(I,NP1).GT.TAU*SIZE)) GO TO 300
        MODE = MODE + 2
        GO TO 450
  300   CONTINUE
  310 CONTINUE
  320 CONTINUE
  330 IF (.NOT.(KRANKE.GT.0)) GO TO 420
C
C     REPLACE DIAG. TERMS OF LOWER TRAP. MATRIX.
      CALL DCOPY(KRANKE, WS(KRANKE+1), 1, W, MDW+1)
C
C     REAPPLY TRANS TO PUT SOLN IN ORIGINAL COORDINATES.
      DO 340 II=1,KRANKE
        I = KRANKE + 1 - II
        CALL H12(2, I, I+1, N, W(I,1), MDW, WS(I), X, 1, 1, 1)
  340 CONTINUE
C
C     COMPUTE COV MATRIX OF EQUAL. CONSTRAINED PROBLEM.
      IF (.NOT.(COV)) GO TO 410
      DO 400 JJ=1,KRANKE
        J = KRANKE + 1 - JJ
        IF (.NOT.(J.LT.N)) GO TO 390
        RB = WS(J)*W(J,J)
        IF (RB.NE.ZERO) RB = ONE/RB
        JP1 = J + 1
        DO 350 I=JP1,N
          W(I,J) = DDOT(N-J,W(I,JP1),MDW,W(J,JP1),MDW)*RB
  350   CONTINUE
        GAM = DDOT(N-J,W(JP1,J),1,W(J,JP1),MDW)*RB
        GAM = HALF*GAM
        CALL DAXPY(N-J, GAM, W(J,JP1), MDW, W(JP1,J), 1)
        DO 370 I=JP1,N
          DO 360 K=I,N
            W(I,K) = W(I,K) + W(J,I)*W(K,J) + W(I,J)*W(J,K)
            W(K,I) = W(I,K)
  360     CONTINUE
  370   CONTINUE
        UJ = WS(J)
        VJ = GAM*UJ
        W(J,J) = UJ*VJ + UJ*VJ
        DO 380 I=JP1,N
          W(J,I) = UJ*W(I,J) + VJ*W(J,I)
  380   CONTINUE
        CALL DCOPY(N-J, W(J,JP1), MDW, W(JP1,J), 1)
  390   CONTINUE
  400 CONTINUE
  410 CONTINUE
C
C     APPLY THE SCALING TO THE COVARIANCE MATRIX.
  420 IF (.NOT.(COV)) GO TO 440
      DO 430 I=1,N
        L = N1 + I
        CALL DSCAL(N, WS(L-1), W(I,1), MDW)
        CALL DSCAL(N, WS(L-1), W(1,I), 1)
  430 CONTINUE
  440 CONTINUE
  450 CONTINUE
C
C     RESCALE SOLN. VECTOR.
      IF (.NOT.(MODE.LE.1)) GO TO 470
      DO 460 J=1,N
        L = N1 + J
        X(J) = X(J)*WS(L-1)
  460 CONTINUE
  470 IP(1) = KRANKE
      IP(3) = IP(3) + 2*KRANKE + N
      RETURN
  480 CONTINUE
C     TO PROCESS-OPTION-VECTOR
C
C     THE NOMINAL TOLERANCE USED IN THE CODE
C     FOR THE EQUALITY CONSTRAINT EQUATIONS.
      TAU = DSQRT(DRELPR)
C
C     THE NOMINAL COLUMN SCALING USED IN THE CODE IS
C     THE IDENTITY SCALING.
      WS(N1) = ONE
      CALL DCOPY(N, WS(N1), 0, WS(N1), 1)
C
C     NO COVARIANCE MATRIX IS NOMINALLY COMPUTED.
      COV = .FALSE.
C
C     DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE.
      NOPT = 1000
      NTIMES = 0
C
C     DEFINE BOUND FOR POSITIVE VALUES OF LINK.
      NLINK = 100000
      LAST = 1
      LINK = PRGOPT(1)
      IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 490
      NERR = 3
      IOPT = 1
      CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
      MODE = 4
      RETURN
  490 IF (.NOT.(LINK.GT.1)) GO TO 540
      NTIMES = NTIMES + 1
      IF (.NOT.(NTIMES.GT.NOPT)) GO TO 500
      NERR = 3
      IOPT = 1
      CALL XERROR(
     * 52HLSEI( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 52,
     * NERR, IOPT)
      MODE = 4
      RETURN
  500 KEY = PRGOPT(LAST+1)
      IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO
      IF (.NOT.(KEY.EQ.2 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 520
      DO 510 J=1,N
        T = DNRM2(M,W(1,J),1)
        IF (T.NE.ZERO) T = ONE/T
        L = N1 + J
        WS(L-1) = T
  510 CONTINUE
  520 IF (KEY.EQ.3) CALL DCOPY(N, PRGOPT(LAST+2), 1, WS(N1), 1)
      IF (KEY.EQ.4) TAU = DMAX1(DRELPR,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 530
      NERR = 3
      IOPT = 1
      CALL XERROR(38HLSEI( ) THE OPTION VECTOR IS UNDEFINED, 38, NERR,
     * IOPT)
      MODE = 4
      RETURN
  530 LAST = LINK
      LINK = NEXT
      GO TO 490
  540 DO 550 J=1,N
        L = N1 + J
        CALL DSCAL(M, WS(L-1), W(1,J), 1)
  550 CONTINUE
      GO TO 560
  560 GO TO IGO990, (90)
      END

      SUBROUTINE LSI(W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, IP) LSI   10
      DOUBLE PRECISION W(MDW,1), PRGOPT(1), RNORM, WS(1), X(1)
CCC
CCC    error: subscript out of bounds for IP(2) in original version
CCC
C      INTEGER IP(1)
      INTEGER IP(2)
      DOUBLE PRECISION ANORM, DRELPR, FAC, GAM, HALF, ONE, RB, TAU, TOL,
     * XNORM, ZERO
      DOUBLE PRECISION DASUM, DDOT, DSQRT, DMAX1
      LOGICAL COV
C
      DATA ZERO /0.D0/, DRELPR /0.D0/, ONE /1.D0/, HALF /.5D0/
C
C     COMPUTE MACHINE PRECISION=DRELPR ONLY WHEN NECESSARY.
      IF (.NOT.(DRELPR.EQ.ZERO)) GO TO 30
C     modification by w.h. 2/98, see hist9 for old version
      CALL GETPR(DRELPR)
   30 MODE = 0
      RNORM = ZERO
      M = MA + MG
      NP1 = N + 1
      KRANK = 0
      IF (N.LE.0 .OR. M.LE.0) GO TO 70
      ASSIGN 40 TO IGO994
      GO TO 500
C
C     PROCESS-OPTION-VECTOR
C
C     COMPUTE MATRIX NORM OF LEAST SQUARES EQUAS.
   40 ANORM = ZERO
      DO 50 J=1,N
        ANORM = DMAX1(ANORM,DASUM(MA,W(1,J),1))
   50 CONTINUE
C
C     SET TOL FOR HFTI( ) RANK TEST.
      TAU = TOL*ANORM
C
C     COMPUTE HOUSEHOLDER ORTHOGONAL DECOMP OF MATRIX.
      IF (N.GT.0) WS(1) = ZERO
      CALL DCOPY(N, WS, 0, WS, 1)
      CALL DCOPY(MA, W(1,NP1), 1, WS, 1)
      K = MAX0(M,N)
      MINMAN = MIN0(MA,N)
      N1 = K + 1
      N2 = N1 + N
      CALL HFTI(W, MDW, MA, N, WS, 1, 1, TAU, KRANK, RNORM, WS(N2),
     * WS(N1), IP)
      FAC = ONE
      GAM=MA-KRANK
      IF (KRANK.LT.MA) FAC = RNORM**2/GAM
      ASSIGN 60 TO IGO990
      GO TO 80
C
C     REDUCE-TO-LPDP-AND-SOLVE
   60 CONTINUE
   70 IP(1) = KRANK
      IP(2) = N + MAX0(M,N) + (MG+2)*(N+7)
      RETURN
   80 CONTINUE
C
C     TO REDUCE-TO-LPDP-AND-SOLVE
      MAP1 = MA + 1
C
C     COMPUTE INEQ. RT-HAND SIDE FOR LPDP.
      IF (.NOT.(MA.LT.M)) GO TO 260
      IF (.NOT.(MINMAN.GT.0)) GO TO 160
      DO 90 I=MAP1,M
        W(I,NP1) = W(I,NP1) - DDOT(N,W(I,1),MDW,WS,1)
   90 CONTINUE
      DO 100 I=1,MINMAN
        J = IP(I)
C
C     APPLY PERMUTATIONS TO COLS OF INEQ. CONSTRAINT MATRIX.
        CALL DSWAP(MG, W(MAP1,I), 1, W(MAP1,J), 1)
  100 CONTINUE
C
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO CONSTRAINT MATRIX.
      IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 120
      DO 110 II=1,KRANK
        I = KRANK + 1 - II
        L = N1 + I
        CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), W(MAP1,1),
     *   MDW, 1, MG)
  110 CONTINUE
C
C     COMPUTE PERMUTED INEQ. CONSTR. MATRIX TIMES R-INVERSE.
  120 DO 150 I=MAP1,M
        IF (.NOT.(0.LT.KRANK)) GO TO 140
        DO 130 J=1,KRANK
          W(I,J) = (W(I,J)-DDOT(J-1,W(1,J),1,W(I,1),MDW))/W(J,J)
  130   CONTINUE
  140   CONTINUE
  150 CONTINUE
C
C     SOLVE THE REDUCED PROBLEM WITH LPDP ALGORITHM,
C     THE LEAST PROJECTED DISTANCE PROBLEM.
  160 CALL LPDP(W(MAP1,1), MDW, MG, KRANK, N-KRANK, PRGOPT, X, XNORM,
     * MDLPDP, WS(N2), IP(N+1))
      IF (.NOT.(MDLPDP.EQ.1)) GO TO 240
      IF (.NOT.(KRANK.GT.0)) GO TO 180
C
C     COMPUTE SOLN IN ORIGINAL COORDINATES.
      DO 170 II=1,KRANK
        I = KRANK + 1 - II
        X(I) = (X(I)-DDOT(II-1,W(I,I+1),MDW,X(I+1),1))/W(I,I)
  170 CONTINUE
C
C     APPLY HOUSEHOLDER TRANS. TO SOLN VECTOR.
  180 IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.N)) GO TO 200
      DO 190 I=1,KRANK
        L = N1 + I
        CALL H12(2, I, KRANK+1, N, W(I,1), MDW, WS(L-1), X, 1, 1, 1)
  190 CONTINUE
  200 IF (.NOT.(MINMAN.GT.0)) GO TO 230
C
C     REPERMUTE VARIABLES TO THEIR INPUT ORDER.
      DO 210 II=1,MINMAN
        I = MINMAN + 1 - II
        J = IP(I)
        CALL DSWAP(1, X(I), 1, X(J), 1)
  210 CONTINUE
C
C     VARIABLES ARE NOW IN ORIG. COORDINATES.
C     ADD SOLN OF UNSCONSTRAINED PROB.
      DO 220 I=1,N
        X(I) = X(I) + WS(I)
  220 CONTINUE
C
C     COMPUTE THE RESIDUAL VECTOR NORM.
      RNORM = DSQRT(RNORM**2+XNORM**2)
  230 GO TO 250
  240 MODE = 2
  250 GO TO 270
  260 CALL DCOPY(N, WS, 1, X, 1)
  270 IF (.NOT.(COV .AND. KRANK.GT.0)) GO TO 490
C
C     COMPUTE COVARIANCE MATRIX BASED ON THE ORTHOGONAL DECOMP.
C     FROM HFTI( ).
C
      KRM1 = KRANK - 1
      KRP1 = KRANK + 1
C
C     COPY DIAG. TERMS TO WORKING ARRAY.
      CALL DCOPY(KRANK, W, MDW+1, WS(N2), 1)
C
C     RECIPROCATE DIAG. TERMS.
      DO 280 J=1,KRANK
        W(J,J) = ONE/W(J,J)
  280 CONTINUE
      IF (.NOT.(KRANK.GT.1)) GO TO 310
C
C     INVERT THE UPPER TRIANGULAR QR FACTOR ON ITSELF.
      DO 300 I=1,KRM1
        IP1 = I + 1
        DO 290 J=IP1,KRANK
          W(I,J) = -DDOT(J-I,W(I,I),MDW,W(I,J),1)*W(J,J)
  290   CONTINUE
  300 CONTINUE
C
C     COMPUTE THE INVERTED FACTOR TIMES ITS TRANSPOSE.
  310 DO 330 I=1,KRANK
        DO 320 J=I,KRANK
          W(I,J) = DDOT(KRANK+1-J,W(I,J),MDW,W(J,J),MDW)
  320   CONTINUE
  330 CONTINUE
      IF (.NOT.(KRANK.LT.N)) GO TO 450
C
C     ZERO OUT LOWER TRAPEZOIDAL PART.
C     COPY UPPER TRI. TO LOWER TRI. PART.
      DO 340 J=1,KRANK
        CALL DCOPY(J, W(1,J), 1, W(J,1), MDW)
  340 CONTINUE
      DO 350 I=KRP1,N
        W(I,1) = ZERO
        CALL DCOPY(I, W(I,1), 0, W(I,1), MDW)
  350 CONTINUE
C
C     APPLY RIGHT SIDE TRANSFORMATIONS TO LOWER TRI.
      N3 = N2 + KRP1
      DO 430 I=1,KRANK
        L = N1 + I
        K = N2 + I
        RB = WS(L-1)*WS(K-1)
        IF (.NOT.(RB.LT.ZERO)) GO TO 420
C
C     IF RB.GE.ZERO, TRANSFORMATION CAN BE REGARDED AS ZERO.
        RB = ONE/RB
C
C     STORE UNSCALED RANK-ONE HOUSEHOLDER UPDATE IN WORK ARRAY.
        WS(N3) = ZERO
        CALL DCOPY(N, WS(N3), 0, WS(N3), 1)
        L = N1 + I
        K = N3 + I
        WS(K-1) = WS(L-1)
        DO 360 J=KRP1,N
          K = N3 + J
          WS(K-1) = W(I,J)
  360   CONTINUE
        DO 370 J=1,N
          L = N3 + I
          K = N3 + J
          WS(J) = DDOT(J-I,W(J,I),MDW,WS(L-1),1) + DDOT(N-J+1,W(J,J),1,
     *     WS(K-1),1)
          WS(J) = WS(J)*RB
  370   CONTINUE
        L = N3 + I
        GAM = DDOT(N-I+1,WS(L-1),1,WS(I),1)*RB
        GAM = GAM*HALF
        CALL DAXPY(N-I+1, GAM, WS(L-1), 1, WS(I), 1)
        DO 410 J=I,N
          IF (.NOT.(I.GT.1)) GO TO 390
          IM1 = I - 1
          K = N3 + J
          DO 380 L=1,IM1
            W(J,L) = W(J,L) + WS(K-1)*WS(L)
  380     CONTINUE
  390     K = N3 + J
          DO 400 L=I,J
            IL = N3 + L
            W(J,L) = W(J,L) + WS(J)*WS(IL-1) + WS(L)*WS(K-1)
  400     CONTINUE
  410   CONTINUE
  420   CONTINUE
  430 CONTINUE
C
C     COPY LOWER TRI. TO UPPER TRI. TO SYMMETRIZE THE COVARIANCE MATRIX.
      DO 440 I=1,N
        CALL DCOPY(I, W(I,1), MDW, W(1,I), 1)
  440 CONTINUE
C
C     REPERMUTE ROWS AND COLS.
  450 DO 470 II=1,MINMAN
        I = MINMAN + 1 - II
        K = IP(I)
        IF (.NOT.(I.NE.K)) GO TO 460
        CALL DSWAP(1, W(I,I), 1, W(K,K), 1)
        CALL DSWAP(I-1, W(1,I), 1, W(1,K), 1)
        CALL DSWAP(K-I-1, W(I,I+1), MDW, W(I+1,K), 1)
        CALL DSWAP(N-K, W(I,K+1), MDW, W(K,K+1), MDW)
  460   CONTINUE
  470 CONTINUE
C
C     PUT IN NORMALIZED RESIDUAL SUM OF SQUARES SCALE FACTOR
C     AND SYMMETRIZE THE RESULTING COVARIANCE MARIX.
      DO 480 J=1,N
        CALL DSCAL(J, FAC, W(1,J), 1)
        CALL DCOPY(J, W(1,J), 1, W(J,1), MDW)
  480 CONTINUE
  490 GO TO 540
  500 CONTINUE
C
C     TO PROCESS-OPTION-VECTOR
C
C     THE NOMINAL TOLERANCE USED IN THE CODE,
      TOL = DSQRT(DRELPR)
      COV = .FALSE.
      LAST = 1
      LINK = PRGOPT(1)
  510 IF (.NOT.(LINK.GT.1)) GO TO 520
      KEY = PRGOPT(LAST+1)
      IF (KEY.EQ.1) COV = PRGOPT(LAST+2).NE.ZERO
      IF (KEY.EQ.5) TOL = DMAX1(DRELPR,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      LAST = LINK
      LINK = NEXT
      GO TO 510
  520 GO TO 530
  530 GO TO IGO994, (40)
  540 GO TO IGO990, (60)
      END
      SUBROUTINE LPDP(A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, IS)LPDP  10
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START EDITING AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/DNRM2/DNRM2/,/DDOT/DDOT/,
C     /DCOPY/DCOPY/,/DSCAL/DSCAL/,/DABS(/DABS(/, ABS/, DABS/,/D0/D0/
C
C     DIMENSION A(MDA,N+1),PRGOPT(*),X(N),WS((M+2)*(N+7)),IS(M+N+1),
C     WHERE N=N1+N2.  THIS IS A SLIGHT OVERESTIMATE FOR WS(*).
C
C     WRITTEN BY R. J. HANSON AND K. H. HASKELL, SANDIA LABS
C     REVISED OCT. 1, 1981.
C
C     DETERMINE AN N1-VECTOR W, AND
C               AN N2-VECTOR Z
C     WHICH MINIMIZES THE EUCLIDEAN LENGTH OF W
C     SUBJECT TO G*W+H*Z .GE. Y.
C     THIS IS THE LEAST PROJECTED DISTANCE PROBLEM, LPDP.
C     THE MATRICES G AND H ARE OF RESPECTIVE
C     DIMENSIONS M BY N1 AND M BY N2.
C
C     CALLED BY SUBPROGRAM LSI( ).
C
C     THE MATRIX
C                (G H Y)
C
C     OCCUPIES ROWS 1,...,M AND COLS 1,...,N1+N2+1 OF A(*,*).
C
C     THE SOLUTION (W) IS RETURNED IN X(*).
C                  (Z)
C
C     THE VALUE OF MODE INDICATES THE STATUS OF
C     THE COMPUTATION AFTER RETURNING TO THE USER.
C
C          MODE=1  THE SOLUTION WAS SUCCESSFULLY OBTAINED.
C
C          MODE=2  THE INEQUALITIES ARE INCONSISTENT.
C
C     SUBROUTINES CALLED
C
C     WNNLS         SOLVES A NONNEGATIVELY CONSTRAINED LINEAR LEAST
C                   SQUARES PROBLEM WITH LINEAR EQUALITY CONSTRAINTS.
C                   PART OF THIS PACKAGE.
C
C++
C     DDOT,         SUBROUTINES FROM THE BLAS PACKAGE.
C     DSCAL,DNRM2,  SEE TRANS. MATH. SOFT., VOL. 5, NO. 3, P. 308.
C     DCOPY
C
      DOUBLE PRECISION A(MDA,1), PRGOPT(1), WS(1), WNORM, X(1)
CCC
CCC    error: subscript out of bounds for IS(2) in original version
CCC
C      INTEGER IS(1)
      INTEGER IS(2)
      DOUBLE PRECISION FAC, ONE, RNORM, SC, YNORM, ZERO
      DOUBLE PRECISION DDOT, DNRM2, ABS
      DATA ZERO, ONE /0.D0,1.D0/, FAC /0.1D0/
      N = N1 + N2
      MODE = 1
      IF (.NOT.(M.LE.0)) GO TO 20
      IF (.NOT.(N.GT.0)) GO TO 10
      X(1) = ZERO
      CALL DCOPY(N, X, 0, X, 1)
   10 WNORM = ZERO
      RETURN
   20 NP1 = N + 1
C
C     SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE.
      DO 40 I=1,M
        SC = DNRM2(N,A(I,1),MDA)
        IF (.NOT.(SC.NE.ZERO)) GO TO 30
        SC = ONE/SC
        CALL DSCAL(NP1, SC, A(I,1), MDA)
   30   CONTINUE
   40 CONTINUE
C
C     SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO).
      YNORM = DNRM2(M,A(1,NP1),1)
      IF (.NOT.(YNORM.NE.ZERO)) GO TO 50
      SC = ONE/YNORM
      CALL DSCAL(M, SC, A(1,NP1), 1)
C
C     SCALE COLS OF MATRIX H.
   50 J = N1 + 1
   60 IF (.NOT.(J.LE.N)) GO TO 70
      SC = DNRM2(M,A(1,J),1)
      IF (SC.NE.ZERO) SC = ONE/SC
      CALL DSCAL(M, SC, A(1,J), 1)
      X(J) = SC
      J = J + 1
      GO TO 60
   70 IF (.NOT.(N1.GT.0)) GO TO 130
C
C     COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*).
      IW = 0
      DO 80 I=1,M
C
C     MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY.
        CALL DCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1)
        IW = IW + N2
C
C     MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY.
        CALL DCOPY(N1, A(I,1), MDA, WS(IW+1), 1)
        IW = IW + N1
C
C     MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY.
        WS(IW+1) = A(I,NP1)
        IW = IW + 1
   80 CONTINUE
      WS(IW+1) = ZERO
      CALL DCOPY(N, WS(IW+1), 0, WS(IW+1), 1)
      IW = IW + N
      WS(IW+1) = ONE
      IW = IW + 1
C
C     SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U.GE.0.  THE
C     MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR
C     F = TRANSPOSE OF (0,...,0,1).
      IX = IW + 1
      IW = IW + M
C
C     DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ).
      IS(1) = 0
      IS(2) = 0
      CALL WNNLS(WS, NP1, N2, NP1-N2, M, 0, PRGOPT, WS(IX), RNORM,
     * MODEW, IS, WS(IW+1))
C
C     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W.
      SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1)
      IF (.NOT.(ONE+FAC*DABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 110
      SC = ONE/SC
      DO 90 J=1,N1
        X(J) = SC*DDOT(M,A(1,J),1,WS(IX),1)
   90 CONTINUE
C
C     COMPUTE THE VECTOR Q=Y-GW.  OVERWRITE Y WITH THIS VECTOR.
      DO 100 I=1,M
        A(I,NP1) = A(I,NP1) - DDOT(N1,A(I,1),MDA,X,1)
  100 CONTINUE
      GO TO 120
  110 MODE = 2
      RETURN
  120 CONTINUE
  130 IF (.NOT.(N2.GT.0)) GO TO 180
C
C     COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*).
      IW = 0
      DO 140 I=1,M
        CALL DCOPY(N2, A(I,N1+1), MDA, WS(IW+1), 1)
        IW = IW + N2
        WS(IW+1) = A(I,NP1)
        IW = IW + 1
  140 CONTINUE
      WS(IW+1) = ZERO
      CALL DCOPY(N2, WS(IW+1), 0, WS(IW+1), 1)
      IW = IW + N2
      WS(IW+1) = ONE
      IW = IW + 1
      IX = IW + 1
      IW = IW + M
C
C     SOLVE RV=S SUBJECT TO V.GE.0.  THE MATRIX R =(TRANSPOSE
C     OF (H Q)), WHERE Q=Y-GW.  THE (N2+1)-VECTOR S =(TRANSPOSE
C     OF (0,...,0,1)).
C
C     DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF WNNLS( ).
      IS(1) = 0
      IS(2) = 0
      CALL WNNLS(WS, N2+1, 0, N2+1, M, 0, PRGOPT, WS(IX), RNORM, MODEW,
     * IS, WS(IW+1))
C
C     COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z.
      SC = ONE - DDOT(M,A(1,NP1),1,WS(IX),1)
      IF (.NOT.(ONE+FAC*DABS(SC).NE.ONE .AND. RNORM.GT.ZERO)) GO TO 160
      SC = ONE/SC
      DO 150 J=1,N2
        L = N1 + J
        X(L) = SC*DDOT(M,A(1,L),1,WS(IX),1)*X(L)
  150 CONTINUE
      GO TO 170
  160 MODE = 2
      RETURN
  170 CONTINUE
C
C     ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION.
  180 CALL DSCAL(N, YNORM, X, 1)
      WNORM = DNRM2(N1,X,1)
      RETURN
      END
      SUBROUTINE WNNLS(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,    WNN   10
     * IWORK, WORK)
C
C     DIMENSION W(MDW,N+1),PRGOPT(*),X(N),IWORK(M+N),WORK(M+5*N)
C
C     ABSTRACT
C
C     THIS SUBPROGRAM SOLVES A LINEARLY CONSTRAINED LEAST SQUARES
C     PROBLEM.  SUPPOSE THERE ARE GIVEN MATRICES E AND A OF
C     RESPECTIVE DIMENSIONS ME BY N AND MA BY N, AND VECTORS F
C     AND B OF RESPECTIVE LENGTHS ME AND MA.  THIS SUBROUTINE
C     SOLVES THE PROBLEM
C
C               EX = F, (EQUATIONS TO BE EXACTLY SATISFIED)
C
C               AX = B, (EQUATIONS TO BE APPROXIMATELY SATISFIED,
C                        IN THE LEAST SQUARES SENSE)
C
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE
C
C     ANY VALUES ME.GE.0, MA.GE.0 AND 0.LE. L .LE.N ARE PERMITTED.
C
C     THE PROBLEM IS REPOSED AS PROBLEM WNNLS
C
C               (WT*E)X = (WT*F)
C               (   A)    (   B), (LEAST SQUARES)
C               SUBJECT TO COMPONENTS L+1,...,N NONNEGATIVE.
C
C     THE SUBPROGRAM CHOOSES THE HEAVY WEIGHT (OR PENALTY PARAMETER) WT.
C
C     THE PARAMETERS FOR WNNLS ARE
C
C     INPUT..
C
C     W(*,*),MDW,  THE ARRAY W(*,*) IS DOUBLE SUBSCRIPTED WITH FIRST
C     ME,MA,N,L    DIMENSIONING PARAMETER EQUAL TO MDW.  FOR THIS
C                  DISCUSSION LET US CALL M = ME + MA.  THEN MDW
C                  MUST SATISFY MDW.GE.M.  THE CONDITION MDW.LT.M
C                  IS AN ERROR.
C
C                  THE ARRAY W(*,*) CONTAINS THE MATRICES AND VECTORS
C
C                       (E  F)
C                       (A  B)
C
C                  IN ROWS AND COLUMNS 1,...,M AND 1,...,N+1
C                  RESPECTIVELY.  COLUMNS 1,...,L CORRESPOND TO
C                  UNCONSTRAINED VARIABLES X(1),...,X(L).  THE
C                  REMAINING VARIABLES ARE CONSTRAINED TO BE
C                  NONNEGATIVE.  THE CONDITION L.LT.0 .OR. L.GT.N IS
C                  AN ERROR.
C
C     PRGOPT(*)    THIS ARRAY IS THE OPTION VECTOR.
C                  IF THE USER IS SATISFIED WITH THE NOMINAL
C                  SUBPROGRAM FEATURES SET
C
C                  PRGOPT(1)=1 (OR PRGOPT(1)=1.0)
C
C                  OTHERWISE PRGOPT(*) IS A LINKED LIST CONSISTING OF
C                  GROUPS OF DATA OF THE FOLLOWING FORM
C
C                  LINK
C                  KEY
C                  DATA SET
C
C                  THE PARAMETERS LINK AND KEY ARE EACH ONE WORD.
C                  THE DATA SET CAN BE COMPRISED OF SEVERAL WORDS.
C                  THE NUMBER OF ITEMS DEPENDS ON THE VALUE OF KEY.
C                  THE VALUE OF LINK POINTS TO THE FIRST
C                  ENTRY OF THE NEXT GROUP OF DATA WITHIN
C                  PRGOPT(*).  THE EXCEPTION IS WHEN THERE ARE
C                  NO MORE OPTIONS TO CHANGE.  IN THAT
C                  CASE LINK=1 AND THE VALUES KEY AND DATA SET
C                  ARE NOT REFERENCED. THE GENERAL LAYOUT OF
C                  PRGOPT(*) IS AS FOLLOWS.
C
C               ...PRGOPT(1)=LINK1 (LINK TO FIRST ENTRY OF NEXT GROUP)
C               .  PRGOPT(2)=KEY1 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(3)=DATA VALUE (DATA VALUE FOR THIS CHANGE)
C               .       .
C               .       .
C               .       .
C               ...PRGOPT(LINK1)=LINK2 (LINK TO THE FIRST ENTRY OF
C               .                       NEXT GROUP)
C               .  PRGOPT(LINK1+1)=KEY2 (KEY TO THE OPTION CHANGE)
C               .  PRGOPT(LINK1+2)=DATA VALUE
C               ...     .
C               .       .
C               .       .
C               ...PRGOPT(LINK)=1 (NO MORE OPTIONS TO CHANGE)
C
C                  VALUES OF LINK THAT ARE NONPOSITIVE ARE ERRORS.
C                  A VALUE OF LINK.GT.NLINK=100000 IS ALSO AN ERROR.
C                  THIS HELPS PREVENT USING INVALID BUT POSITIVE
C                  VALUES OF LINK THAT WILL PROBABLY EXTEND
C                  BEYOND THE PROGRAM LIMITS OF PRGOPT(*).
C                  UNRECOGNIZED VALUES OF KEY ARE IGNORED.  THE
C                  ORDER OF THE OPTIONS IS ARBITRARY AND ANY NUMBER
C                  OF OPTIONS CAN BE CHANGED WITH THE FOLLOWING
C                  RESTRICTION.  TO PREVENT CYCLING IN THE
C                  PROCESSING OF THE OPTION ARRAY A COUNT OF THE
C                  NUMBER OF OPTIONS CHANGED IS MAINTAINED.
C                  WHENEVER THIS COUNT EXCEEDS NOPT=1000 AN ERROR
C                  MESSAGE IS PRINTED AND THE SUBPROGRAM RETURNS.
C
C                  OPTIONS..
C
C                  KEY=6
C                         SCALE THE NONZERO COLUMNS OF THE
C                  ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  TO HAVE LENGTH ONE.  THE DATA SET FOR
C                  THIS OPTION IS A SINGLE VALUE.  IT MUST
C                  BE NONZERO IF UNIT LENGTH COLUMN SCALING IS
C                  DESIRED.
C
C                  KEY=7
C                         SCALE COLUMNS OF THE ENTIRE DATA MATRIX
C                  (E)
C                  (A)
C                  WITH A USER-PROVIDED DIAGONAL MATRIX.
C                  THE DATA SET FOR THIS OPTION CONSISTS
C                  OF THE N DIAGONAL SCALING FACTORS, ONE FOR
C                  EACH MATRIX COLUMN.
C
C                  KEY=8
C                         CHANGE THE RANK DETERMINATION TOLERANCE FROM
C                  THE NOMINAL VALUE OF DSQRT(EPS).  THIS QUANTITY CAN
C                  BE NO SMALLER THAN EPS, THE ARITHMETIC-
C                  STORAGE PRECISION.  THE QUANTITY USED
C                  HERE IS INTERNALLY RESTRICTED TO BE AT
C                  LEAST EPS.  THE DATA SET FOR THIS OPTION
C                  IS THE NEW TOLERANCE.
C
C                  KEY=9
C                         CHANGE THE BLOW-UP PARAMETER FROM THE
C                  NOMINAL VALUE OF DSQRT(EPS).  THE RECIPROCAL OF
C                  THIS PARAMETER IS USED IN REJECTING SOLUTION
C                  COMPONENTS AS TOO LARGE WHEN A VARIABLE IS
C                  FIRST BROUGHT INTO THE ACTIVE SET.  TOO LARGE
C                  MEANS THAT THE PROPOSED COMPONENT TIMES THE
C                  RECIPROCAL OF THE PARAMETERIS NOT LESS THAN
C                  THE RATIO OF THE NORMS OF THE RIGHT-SIDE
C                  VECTOR AND THE DATA MATRIX.
C                  THIS PARAMETER CAN BE NO SMALLER THAN EPS,
C                  THE ARITHMETIC-STORAGE PRECISION.
C
C                  FOR EXAMPLE, SUPPOSE WE WANT TO PROVIDE
C                  A DIAGONAL MATRIX TO SCALE THE PROBLEM
C                  MATRIX AND CHANGE THE TOLERANCE USED FOR
C                  DETERMINING LINEAR DEPENDENCE OF DROPPED COL
C                  VECTORS.  FOR THESE OPTIONS THE DIMENSIONS OF
C                  PRGOPT(*) MUST BE AT LEAST N+6.  THE FORTRAN
C                  STATEMENTS DEFINING THESE OPTIONS WOULD
C                  BE AS FOLLOWS.
C
C                  PRGOPT(1)=N+3 (LINK TO ENTRY N+3 IN PRGOPT(*))
C                  PRGOPT(2)=7 (USER-PROVIDED SCALING KEY)
C
C                  CALL DCOPY(N,D,1,PRGOPT(3),1) (COPY THE N
C                  SCALING FACTORS FROM A USER ARRAY CALLED D(*)
C                  INTO PRGOPT(3)-PRGOPT(N+2))
C
C                  PRGOPT(N+3)=N+6 (LINK TO ENTRY N+6 OF PRGOPT(*))
C                  PRGOPT(N+4)=8 (LINEAR DEPENDENCE TOLERANCE KEY)
C                  PRGOPT(N+5)=... (NEW VALUE OF THE TOLERANCE)
C
C                  PRGOPT(N+6)=1 (NO MORE OPTIONS TO CHANGE)
C
C     IWORK(1),    THE AMOUNTS OF WORKING STORAGE ACTUALLY ALLOCATED
C     IWORK(2)     FOR THE WORKING ARRAYS WORK(*) AND IWORK(*),
C                  RESPECTIVELY.  THESE QUANTITIES ARE COMPARED WITH
C                  THE ACTUAL AMOUNTS OF STORAGE NEEDED FOR WNNLS( ).
C                  INSUFFICIENT STORAGE ALLOCATED FOR EITHER WORK(*)
C                  OR IWORK(*) IS CONSIDERED AN ERROR.  THIS FEATURE
C                  WAS INCLUDED IN WNNLS( ) BECAUSE MISCALCULATING
C                  THE STORAGE FORMULAS FOR WORK(*) AND IWORK(*)
C                  MIGHT VERY WELL LEAD TO SUBTLE AND HARD-TO-FIND
C                  EXECUTION ERRORS.
C
C                  THE LENGTH OF WORK(*) MUST BE AT LEAST
C
C                  LW = ME+MA+5*N
C                  THIS TEST WILL NOT BE MADE IF IWORK(1).LE.0.
C
C                  THE LENGTH OF IWORK(*) MUST BE AT LEAST
C
C                  LIW = ME+MA+N
C                  THIS TEST WILL NOT BE MADE IF IWORK(2).LE.0.
C
C     OUTPUT..
C
C     X(*)         AN ARRAY DIMENSIONED AT LEAST N, WHICH WILL
C                  CONTAIN THE N COMPONENTS OF THE SOLUTION VECTOR
C                  ON OUTPUT.
C
C     RNORM        THE RESIDUAL NORM OF THE SOLUTION.  THE VALUE OF
C                  RNORM CONTAINS THE RESIDUAL VECTOR LENGTH OF THE
C                  EQUALITY CONSTRAINTS AND LEAST SQUARES EQUATIONS.
C
C     MODE         THE VALUE OF MODE INDICATES THE SUCCESS OR FAILURE
C                  OF THE SUBPROGRAM.
C
C                  MODE = 0  SUBPROGRAM COMPLETED SUCCESSFULLY.
C
C                       = 1  MAX. NUMBER OF ITERATIONS (EQUAL TO
C                            3*(N-L)) EXCEEDED. NEARLY ALL PROBLEMS
C                            SHOULD COMPLETE IN FEWER THAN THIS
C                            NUMBER OF ITERATIONS. AN APPROXIMATE
C                            SOLUTION AND ITS CORRESPONDING RESIDUAL
C                            VECTOR LENGTH ARE IN X(*) AND RNORM.
C
C                       = 2  USAGE ERROR OCCURRED.  THE OFFENDING
C                            CONDITION IS NOTED WITH THE ERROR
C                            PROCESSING SUBPROGRAM, XERROR( ).
C
C     USER-DESIGNATED
C     WORKING ARRAYS..
C
C     WORK(*)      A WORKING ARRAY OF LENGTH AT LEAST
C                  M + 5*N.
C
C     IWORK(*)     AN INTEGER-VALUED WORKING ARRAY OF LENGTH AT LEAST
C                  M+N.
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/, DUMMY/,SNGL(DUMMY)/
C
C     WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES,
C     AND R.J. HANSON, SANDIA LABORATORIES.
C     REVISED FEB.25, 1982.
C
C     SUBROUTINES CALLED BY WNNLS( )
C
C++
C     WNLSM         COMPANION SUBROUTINE TO WNNLS( ), WHERE
C                   MOST OF THE COMPUTATION TAKES PLACE.
C
C     XERROR,XERRWV FROM SLATEC ERROR PROCESSING PACKAGE.
C                   THIS IS DOCUMENTED IN SANDIA TECH. REPT.,
C                   SAND78-1189.
C
C     REFERENCES
C
C     1. SOLVING LEAST SQUARES PROBLEMS, BY C.L. LAWSON
C        AND R.J. HANSON.  PRENTICE-HALL, INC. (1974).
C
C     2. BASIC LINEAR ALGEBRA SUBPROGRAMS FOR FORTRAN USAGE, BY
C        C.L. LAWSON, R.J. HANSON, D.R. KINCAID, AND F.T. KROGH.
C        TOMS, V. 5, NO. 3, P. 308.  ALSO AVAILABLE AS
C        SANDIA TECHNICAL REPORT NO. SAND77-0898.
C
C     3. AN ALGORITHM FOR LINEAR LEAST SQUARES WITH EQUALITY
C        AND NONNEGATIVITY CONSTRAINTS, BY K.H. HASKELL AND
C        R.J. HANSON.  AVAILABLE AS SANDIA TECHNICAL REPORT NO.
C        SAND77-0552, AND MATH. PROGRAMMING, VOL. 21, (1981), P. 98-118.
C
C     4. SLATEC COMMON MATH. LIBRARY ERROR HANDLING
C        PACKAGE.  BY R. E. JONES.  AVAILABLE AS SANDIA
C        TECHNICAL REPORT SAND78-1189.
C
      DOUBLE PRECISION  DUMMY, W(MDW,1), PRGOPT(1), X(1), WORK(1), RNORM
CCC
CCC    error subscript out of bounds for IWORK(2) in original version
CCC
C      INTEGER IWORK(1)
      INTEGER IWORK(2)
C
C
      MODE = 0
      IF (MA+ME.LE.0 .OR. N.LE.0) RETURN
      IF (.NOT.(IWORK(1).GT.0)) GO TO 20
      LW = ME + MA + 5*N
      IF (.NOT.(IWORK(1).LT.LW)) GO TO 10
      NERR = 2
      IOPT = 1
      CALL XERRWV(70HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR WORK(*
     *), NEED LW=I1 BELOW, 70, NERR, IOPT, 1, LW, 0, 0, DUMMY, DUMMY)
      MODE = 2
      RETURN
   10 CONTINUE
   20 IF (.NOT.(IWORK(2).GT.0)) GO TO 40
      LIW = ME + MA + N
      IF (.NOT.(IWORK(2).LT.LIW)) GO TO 30
      NERR = 2
      IOPT = 1
      CALL XERRWV(72HWNNLS( ), INSUFFICIENT STORAGE ALLOCATED FOR IWORK(
     **), NEED LIW=I1 BELOW, 72, NERR, IOPT, 1, LIW, 0, 0, DUMMY, DUMMY)
      MODE = 2
      RETURN
   30 CONTINUE
   40 IF (.NOT.(MDW.LT.ME+MA)) GO TO 50
      NERR = 1
      IOPT = 1
      CALL XERROR(44HWNNLS( ), THE VALUE MDW.LT.ME+MA IS AN ERROR, 44,
     * NERR, IOPT)
      MODE = 2
      RETURN
   50 IF (0.LE.L .AND. L.LE.N) GO TO 60
      NERR = 2
      IOPT = 1
      CALL XERROR(39HWNNLS( ), L.LE.0.AND.L.LE.N IS REQUIRED, 39, NERR,
     * IOPT)
      MODE = 2
      RETURN
C
C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS
C     REQUIRED BY THE MAIN SUBROUTINE WNLSM( ).
C
   60 L1 = N + 1
      L2 = L1 + N
      L3 = L2 + ME + MA
      L4 = L3 + N
      L5 = L4 + N
C
      CALL WNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK,
     * IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3), WORK(L4),
     * WORK(L5))
      RETURN
      END
      SUBROUTINE WNLSM(W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE,   WNL   10
     * IPIVOT, ITYPE, WD, H, SCALE, Z, TEMP, D)
C
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START CHANGES AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/DASUM/DASUM/,/DROTMG/DROTMG/,
C     /DNRM2/DNRM2/,/ DSQRT/ DDSQRT/,/DROTM/DROTM/,/DMAX1/DMAX1/,
C     /DCOPY/DCOPY/,/DSCAL/DSCAL/,/DAXPY/DAXPY/,/D0/D0/,/DSWAP/DSWAP/,
C     /IDAMAX/IDAMAX/,/DRELPR/DRELPR/,/.D-/.D-/                   REMK
C
C     THIS IS A COMPANION SUBPROGRAM TO WNNLS( ).
C     THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE
C     USAGE INSTRUCTIONS.
C
C     WRITTEN BY KAREN H. HASKELL, SANDIA LABORATORIES,
C     WITH THE HELP OF R.J. HANSON, SANDIA LABORATORIES,
C     DECEMBER 1976 - JANUARY 1978.
C     REVISED MAR. 4, 1982.
C
C     IN ADDITION TO THE PARAMETERS DISCUSSED IN THE PROLOGUE TO
C     SUBROUTINE WNNLS, THE FOLLOWING WORK ARRAYS ARE USED IN
C     SUBROUTINE WNLSM  (THEY ARE PASSED THROUGH THE CALLING
C     SEQUENCE FROM WNNLS FOR PURPOSES OF VARIABLE DIMENSIONING).
C     THEIR CONTENTS WILL IN GENERAL BE OF NO INTEREST TO THE USER.
C
C         IPIVOT(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         PIVOTING INFORMATION FOR THE COLS OF W(*,*).
C
C         ITYPE(*)
C            AN ARRAY OF LENGTH M WHICH IS USED TO KEEP TRACK
C         OF THE CLASSIFICATION OF THE EQUATIONS.  ITYPE(I)=0
C         DENOTES EQUATION I AS AN EQUALITY CONSTRAINT.
C         ITYPE(I)=1 DENOTES EQUATION I AS A LEAST SQUARES
C         EQUATION.
C
C         WD(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         DUAL SOLUTION VECTOR.
C
C         H(*)
C            AN ARRAY OF LENGTH N.  UPON COMPLETION IT CONTAINS THE
C         PIVOT SCALARS OF THE HOUSEHOLDER TRANSFORMATIONS PERFORMED
C         IN THE CASE KRANK.LT.L.
C
C         SCALE(*)
C            AN ARRAY OF LENGTH M WHICH IS USED BY THE SUBROUTINE
C         TO STORE THE DIAGONAL MATRIX OF WEIGHTS.
C         THESE ARE USED TO APPLY THE MODIFIED GIVENS
C         TRANSFORMATIONS.
C
C         Z(*),TEMP(*)
C            WORKING ARRAYS OF LENGTH N.
C
C         D(*)
C            AN ARRAY OF LENGTH N THAT CONTAINS THE
C         COLUMN SCALING FOR THE MATRIX (E).
C                                       (A)
C
C     SUBROUTINE WNLSM (W,MDW,MME,MA,N,L,PRGOPT,X,RNORM,MODE,
C    1                  IPIVOT,ITYPE,WD,H,SCALE,Z,TEMP,D)
C++
      DOUBLE PRECISION W(MDW,1), X(1), WD(1), H(1), SCALE(1), DOPE(4)
      DOUBLE PRECISION Z(1), TEMP(1), PRGOPT(1), D(1), SPARAM(5)
      DOUBLE PRECISION ALAMDA, ALPHA, ALSQ, AMAX, BNORM, EANORM
      DOUBLE PRECISION DRELPR, FAC, ONE, BLOWUP
      DOUBLE PRECISION RNORM, SM, T, TAU, TWO, WMAX, ZERO, ZZ, Z2
      DOUBLE PRECISION DMAX1, DSQRT, DNRM2, DASUM
      INTEGER IPIVOT(1), ITYPE(1), IDAMAX, IDOPE(8)
      LOGICAL HITCON, FEASBL, DONE, POS
      DATA ZERO /0.D0/, ONE /1.D0/, TWO /2.D0/, DRELPR /0.D0/
C
C     INITIALIZE-VARIABLES
      ASSIGN 10 TO IGO998
      GO TO 180
C
C     PERFORM INITIAL TRIANGULARIZATION IN THE SUBMATRIX
C     CORRESPONDING TO THE UNCONSTRAINED VARIABLES USING
C     THE PROCEDURE INITIALLY-TRIANGULARIZE.
   10 ASSIGN 20 TO IGO995
      GO TO 280
C
C     PERFORM WNNLS ALGORITHM USING THE FOLLOWING STEPS.
C
C     UNTIL(DONE)
C
C        COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
C        WHEN (HITCON) ADD-CONSTRAINTS
C
C        ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
C
C        FIN
C
C     COMPUTE-FINAL-SOLUTION
C
   20 IF (DONE) GO TO 80
C
      ASSIGN 30 TO IGO991
      GO TO 300
C
C     COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
   30 IF (.NOT.(HITCON)) GO TO 50
      ASSIGN 40 TO IGO986
      GO TO 370
   40 GO TO 70
C
C     WHEN (HITCON) ADD-CONSTRAINTS
C
   50 ASSIGN 60 TO IGO983
      GO TO 640
   60 CONTINUE
C
C     ELSE PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
C
   70 GO TO 20
C
   80 ASSIGN 90 TO IGO980
      GO TO 1000
C
C     COMPUTE-FINAL-SOLUTION
C
   90 RETURN
  100 CONTINUE
C
C     TO PROCESS-OPTION-VECTOR
      FAC = 1.D-4
C
C     THE NOMINAL TOLERANCE USED IN THE CODE,
      TAU = DSQRT(DRELPR)
C
C     THE NOMINAL BLOW-UP FACTOR USED IN THE CODE.
      BLOWUP = TAU
C
C     THE NOMINAL COLUMN SCALING USED IN THE CODE IS
C     THE IDENTITY SCALING.
      D(1) = ONE
      CALL DCOPY(N, D, 0, D, 1)
C
C     DEFINE BOUND FOR NUMBER OF OPTIONS TO CHANGE.
      NOPT = 1000
C
C     DEFINE BOUND FOR POSITIVE VALUE OF LINK.
      NLINK = 100000
      NTIMES = 0
      LAST = 1
      LINK = PRGOPT(1)
      IF (.NOT.(LINK.LE.0 .OR. LINK.GT.NLINK)) GO TO 110
      NERR = 3
      IOPT = 1
      CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR,
     * IOPT)
      MODE = 2
      RETURN
  110 IF (.NOT.(LINK.GT.1)) GO TO 160
      NTIMES = NTIMES + 1
      IF (.NOT.(NTIMES.GT.NOPT)) GO TO 120
      NERR = 3
      IOPT = 1
      CALL XERROR(
     * 53HWNNLS( ). THE LINKS IN THE OPTION VECTOR ARE CYCLING., 53,
     * NERR, IOPT)
      MODE = 2
      RETURN
  120 KEY = PRGOPT(LAST+1)
      IF (.NOT.(KEY.EQ.6 .AND. PRGOPT(LAST+2).NE.ZERO)) GO TO 140
      DO 130 J=1,N
        T = DNRM2(M,W(1,J),1)
        IF (T.NE.ZERO) T = ONE/T
        D(J) = T
  130 CONTINUE
  140 IF (KEY.EQ.7) CALL DCOPY(N, PRGOPT(LAST+2), 1, D, 1)
      IF (KEY.EQ.8) TAU = DMAX1(DRELPR,PRGOPT(LAST+2))
      IF (KEY.EQ.9) BLOWUP = DMAX1(DRELPR,PRGOPT(LAST+2))
      NEXT = PRGOPT(LINK)
      IF (.NOT.(NEXT.LE.0 .OR. NEXT.GT.NLINK)) GO TO 150
      NERR = 3
      IOPT = 1
      CALL XERROR(39HWNNLS( ) THE OPTION VECTOR IS UNDEFINED, 39, NERR,
     * IOPT)
      MODE = 2
      RETURN
  150 LAST = LINK
      LINK = NEXT
      GO TO 110
  160 DO 170 J=1,N
        CALL DSCAL(M, D(J), W(1,J), 1)
  170 CONTINUE
      GO TO 1260
  180 CONTINUE
C
C     TO INITIALIZE-VARIABLES
C
C     DRELPR IS THE PRECISION FOR THE PARTICULAR MACHINE
C     BEING USED.  THIS LOGIC AVOIDS RECOMPUTING IT EVERY ENTRY.
      IF (.NOT.(DRELPR.EQ.ZERO)) GO TO 210
C     modification by w.h. 2/98, see hist9 for old version
      CALL GETPR(DRELPR)
  210 M = MA + MME
      ME = MME
      MEP1 = ME + 1
      ASSIGN 220 TO IGO977
      GO TO 100
C
C     PROCESS-OPTION-VECTOR
  220 DONE = .FALSE.
      ITER = 0
      ITMAX = 3*(N-L)
      MODE = 0
      LP1 = L + 1
      NSOLN = L
      NSP1 = NSOLN + 1
      NP1 = N + 1
      NM1 = N - 1
      L1 = MIN0(M,L)
C
C     COMPUTE SCALE FACTOR TO APPLY TO EQUAL. CONSTRAINT EQUAS.
      DO 230 J=1,N
        WD(J) = DASUM(M,W(1,J),1)
  230 CONTINUE
      IMAX = IDAMAX(N,WD,1)
      EANORM = WD(IMAX)
      BNORM = DASUM(M,W(1,NP1),1)
      ALAMDA = EANORM/(DRELPR*FAC)
C
C     DEFINE SCALING DIAG MATRIX FOR MOD GIVENS USAGE AND
C     CLASSIFY EQUATION TYPES.
      ALSQ = ALAMDA**2
      DO 260 I=1,M
C
C     WHEN EQU I IS HEAVILY WEIGHTED ITYPE(I)=0, ELSE ITYPE(I)=1.
        IF (.NOT.(I.LE.ME)) GO TO 240
        T = ALSQ
        ITEMP = 0
        GO TO 250
  240   T = ONE
        ITEMP = 1
  250   SCALE(I) = T
        ITYPE(I) = ITEMP
  260 CONTINUE
C
C     SET THE SOLN VECTOR X(*) TO ZERO AND THE COL INTERCHANGE
C     MATRIX TO THE IDENTITY.
      X(1) = ZERO
      CALL DCOPY(N, X, 0, X, 1)
      DO 270 I=1,N
        IPIVOT(I) = I
  270 CONTINUE
      GO TO 1230
  280 CONTINUE
C
C     TO INITIALLY-TRIANGULARIZE
C
C     SET FIRST L COMPS. OF DUAL VECTOR TO ZERO BECAUSE
C     THESE CORRESPOND TO THE UNCONSTRAINED VARIABLES.
      IF (.NOT.(L.GT.0)) GO TO 290
      WD(1) = ZERO
      CALL DCOPY(L, WD, 0, WD, 1)
C
C     THE ARRAYS IDOPE(*) AND DOPE(*) ARE USED TO PASS
C     INFORMATION TO WNLIT().  THIS WAS DONE TO AVOID
C     A LONG CALLING SEQUENCE OR THE USE OF COMMON.
  290 IDOPE(1) = ME
      IDOPE(2) = MEP1
      IDOPE(3) = 0
      IDOPE(4) = 1
      IDOPE(5) = NSOLN
      IDOPE(6) = 0
      IDOPE(7) = 1
      IDOPE(8) = L1
C
      DOPE(1) = ALSQ
      DOPE(2) = EANORM
      DOPE(3) = FAC
      DOPE(4) = TAU
      CALL WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM,
     * IDOPE, DOPE, DONE)
      ME = IDOPE(1)
      MEP1 = IDOPE(2)
      KRANK = IDOPE(3)
      KRP1 = IDOPE(4)
      NSOLN = IDOPE(5)
      NIV = IDOPE(6)
      NIV1 = IDOPE(7)
      L1 = IDOPE(8)
      GO TO 1240
  300 CONTINUE
C
C     TO COMPUTE-SEARCH-DIRECTION-AND-FEASIBLE-POINT
C
C     SOLVE THE TRIANGULAR SYSTEM OF CURRENTLY NON-ACTIVE
C     VARIABLES AND STORE THE SOLUTION IN Z(*).
C
C     SOLVE-SYSTEM
      ASSIGN 310 TO IGO958
      GO TO 1110
C
C     INCREMENT ITERATION COUNTER AND CHECK AGAINST MAX. NUMBER
C     OF ITERATIONS.
  310 ITER = ITER + 1
      IF (.NOT.(ITER.GT.ITMAX)) GO TO 320
      MODE = 1
      DONE = .TRUE.
C
C     CHECK TO SEE IF ANY CONSTRAINTS HAVE BECOME ACTIVE.
C     IF SO, CALCULATE AN INTERPOLATION FACTOR SO THAT ALL
C     ACTIVE CONSTRAINTS ARE REMOVED FROM THE BASIS.
  320 ALPHA = TWO
      HITCON = .FALSE.
      IF (.NOT.(L.LT.NSOLN)) GO TO 360
      DO 350 J=LP1,NSOLN
        ZZ = Z(J)
        IF (.NOT.(ZZ.LE.ZERO)) GO TO 340
        T = X(J)/(X(J)-ZZ)
        IF (.NOT.(T.LT.ALPHA)) GO TO 330
        ALPHA = T
        JCON = J
  330   HITCON = .TRUE.
  340   CONTINUE
  350 CONTINUE
  360 GO TO 1220
  370 CONTINUE
C
C     TO ADD-CONSTRAINTS
C
C     USE COMPUTED ALPHA TO INTERPOLATE BETWEEN LAST
C     FEASIBLE SOLUTION X(*) AND CURRENT UNCONSTRAINED
C     (AND INFEASIBLE) SOLUTION Z(*).
      IF (.NOT.(LP1.LE.NSOLN)) GO TO 390
      DO 380 J=LP1,NSOLN
        X(J) = X(J) + ALPHA*(Z(J)-X(J))
  380 CONTINUE
  390 FEASBL = .FALSE.
      GO TO 410
  400 IF (FEASBL) GO TO 610
C
C     REMOVE COL JCON AND SHIFT COLS JCON+1 THROUGH N TO THE
C     LEFT. SWAP COL JCON INTO THE N-TH POSITION.  THIS ACHIEVES
C     UPPER HESSENBERG FORM FOR THE NONACTIVE CONSTRAINTS AND
C     LEAVES AN UPPER HESSENBERG MATRIX TO RETRIANGULARIZE.
  410 DO 420 I=1,M
        T = W(I,JCON)
        CALL DCOPY(N-JCON, W(I,JCON+1), MDW, W(I,JCON), MDW)
        W(I,N) = T
  420 CONTINUE
C
C     UPDATE PERMUTED INDEX VECTOR TO REFLECT THIS SHIFT AND SWAP.
      ITEMP = IPIVOT(JCON)
      IF (.NOT.(JCON.LT.N)) GO TO 440
      DO 430 I=JCON,NM1
        IPIVOT(I) = IPIVOT(I+1)
  430 CONTINUE
  440 IPIVOT(N) = ITEMP
C
C     SIMILARLY REPERMUTE X(*) VECTOR.
      CALL DCOPY(N-JCON, X(JCON+1), 1, X(JCON), 1)
      X(N) = ZERO
      NSP1 = NSOLN
      NSOLN = NSOLN - 1
      NIV1 = NIV
      NIV = NIV - 1
C
C     RETRIANGULARIZE UPPER HESSENBERG MATRIX AFTER ADDING CONSTRAINTS.
      J = JCON
      I = KRANK + JCON - L
  450 IF (.NOT.(J.LE.NSOLN)) GO TO 570
      IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0)) GO TO 470
      ASSIGN 460 TO IGO938
      GO TO 620
C
C     (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0) ZERO-IP1-TO-I-IN-COL-J
  460 GO TO 560
  470 IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1)) GO TO 490
      ASSIGN 480 TO IGO938
      GO TO 620
C
C     (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1) ZERO-IP1-TO-I-IN-COL-J
  480 GO TO 560
  490 IF (.NOT.(ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.0)) GO TO 510
      CALL DSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW)
      CALL DSWAP(1, SCALE(I), 1, SCALE(I+1), 1)
      ITEMP = ITYPE(I+1)
      ITYPE(I+1) = ITYPE(I)
      ITYPE(I) = ITEMP
C
C     SWAPPED ROW WAS FORMERLY A PIVOT ELT., SO IT WILL
C     BE LARGE ENOUGH TO PERFORM ELIM.
      ASSIGN 500 TO IGO938
      GO TO 620
C
C     ZERO-IP1-TO-I-IN-COL-J
  500 GO TO 560
  510 IF (.NOT.(ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.1)) GO TO 550
      T = SCALE(I)*W(I,J)**2/ALSQ
      IF (.NOT.(T.GT.TAU**2*EANORM**2)) GO TO 530
      ASSIGN 520 TO IGO938
      GO TO 620
  520 GO TO 540
  530 CALL DSWAP(NP1, W(I,1), MDW, W(I+1,1), MDW)
      CALL DSWAP(1, SCALE(I), 1, SCALE(I+1), 1)
      ITEMP = ITYPE(I+1)
      ITYPE(I+1) = ITYPE(I)
      ITYPE(I) = ITEMP
      W(I+1,J) = ZERO
  540 CONTINUE
  550 CONTINUE
  560 I = I + 1
      J = J + 1
      GO TO 450
C
C     SEE IF THE REMAINING COEFFS IN THE SOLN SET ARE FEASIBLE.  THEY
C     SHOULD BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.  IF ANY ARE
C     INFEASIBLE IT IS DUE TO ROUNDOFF ERROR.  ANY THAT ARE NON-
C     POSITIVE WILL BE SET TO ZERO AND REMOVED FROM THE SOLN SET.
  570 IF (.NOT.(LP1.LE.NSOLN)) GO TO 590
      DO 580 JCON=LP1,NSOLN
        IF (X(JCON).LE.ZERO) GO TO 600
  580 CONTINUE
  590 FEASBL = .TRUE.
  600 CONTINUE
      GO TO 400
  610 GO TO 1200
  620 CONTINUE
C
C     TO ZERO-IP1-TO-I-IN-COL-J
      IF (.NOT.(W(I+1,J).NE.ZERO)) GO TO 630
      CALL DROTMG(SCALE(I), SCALE(I+1), W(I,J), W(I+1,J), SPARAM)
      W(I+1,J) = ZERO
      CALL DROTM(NP1-J, W(I,J+1), MDW, W(I+1,J+1), MDW, SPARAM)
  630 GO TO 1290
  640 CONTINUE
C
C     TO PERFORM-MULTIPLIER-TEST-AND-DROP-A-CONSTRAINT
      CALL DCOPY(NSOLN, Z, 1, X, 1)
      IF (.NOT.(NSOLN.LT.N)) GO TO 650
      X(NSP1) = ZERO
      CALL DCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1)
  650 I = NIV1
  660 IF (.NOT.(I.LE.ME)) GO TO 690
C
C     RECLASSIFY LEAST SQUARES EQATIONS AS EQUALITIES AS
C     NECESSARY.
      IF (.NOT.(ITYPE(I).EQ.0)) GO TO 670
      I = I + 1
      GO TO 680
  670 CALL DSWAP(NP1, W(I,1), MDW, W(ME,1), MDW)
      CALL DSWAP(1, SCALE(I), 1, SCALE(ME), 1)
      ITEMP = ITYPE(I)
      ITYPE(I) = ITYPE(ME)
      ITYPE(ME) = ITEMP
      MEP1 = ME
      ME = ME - 1
  680 GO TO 660
C
C     FORM INNER PRODUCT VECTOR WD(*) OF DUAL COEFFS.
  690 IF (.NOT.(NSP1.LE.N)) GO TO 730
      DO 720 J=NSP1,N
        SM = ZERO
        IF (.NOT.(NSOLN.LT.M)) GO TO 710
        DO 700 I=NSP1,M
          SM = SM + SCALE(I)*W(I,J)*W(I,NP1)
  700   CONTINUE
  710   WD(J) = SM
  720 CONTINUE
  730 GO TO 750
  740 IF (POS .OR. DONE) GO TO 970
C
C     FIND J SUCH THAT WD(J)=WMAX IS MAXIMUM.  THIS DETERMINES
C     THAT THE INCOMING COL J WILL REDUCE THE RESIDUAL VECTOR
C     AND BE POSITIVE.
  750 WMAX = ZERO
      IWMAX = NSP1
      IF (.NOT.(NSP1.LE.N)) GO TO 780
      DO 770 J=NSP1,N
        IF (.NOT.(WD(J).GT.WMAX)) GO TO 760
        WMAX = WD(J)
        IWMAX = J
  760   CONTINUE
  770 CONTINUE
  780 IF (.NOT.(WMAX.LE.ZERO)) GO TO 790
      DONE = .TRUE.
      GO TO 960
C
C     SET DUAL COEFF TO ZERO FOR INCOMING COL.
  790 WD(IWMAX) = ZERO
C
C     WMAX .GT. ZERO, SO OKAY TO MOVE COL IWMAX TO SOLN SET.
C     PERFORM TRANSFORMATION TO RETRIANGULARIZE, AND TEST
C     FOR NEAR LINEAR DEPENDENCE.
C     SWAP COL IWMAX INTO NSOLN-TH POSITION TO MAINTAIN UPPER
C     HESSENBERG FORM OF ADJACENT COLS, AND ADD NEW COL TO
C     TRIANGULAR DECOMPOSITION.
      NSOLN = NSP1
      NSP1 = NSOLN + 1
      NIV = NIV1
      NIV1 = NIV + 1
      IF (.NOT.(NSOLN.NE.IWMAX)) GO TO 800
      CALL DSWAP(M, W(1,NSOLN), 1, W(1,IWMAX), 1)
      WD(IWMAX) = WD(NSOLN)
      WD(NSOLN) = ZERO
      ITEMP = IPIVOT(NSOLN)
      IPIVOT(NSOLN) = IPIVOT(IWMAX)
      IPIVOT(IWMAX) = ITEMP
C
C     REDUCE COL NSOLN SO THAT THE MATRIX OF NONACTIVE
C     CONSTRAINTS VARIABLES IS TRIANGULAR.
  800 J = M
  810 IF (.NOT.(J.GT.NIV)) GO TO 870
      JM1 = J - 1
      JP = JM1
C
C     WHEN OPERATING NEAR THE ME LINE, TEST TO SEE IF THE PIVOT ELT.
C     IS NEAR ZERO.  IF SO, USE THE LARGEST ELT. ABOVE IT AS THE PIVOT.
C     THIS IS TO MAINTAIN THE SHARP INTERFACE BETWEEN WEIGHTED AND
C     NON-WEIGHTED ROWS IN ALL CASES.
      IF (.NOT.(J.EQ.MEP1)) GO TO 850
      IMAX = ME
      AMAX = SCALE(ME)*W(ME,NSOLN)**2
  820 IF (.NOT.(JP.GE.NIV)) GO TO 840
      T = SCALE(JP)*W(JP,NSOLN)**2
      IF (.NOT.(T.GT.AMAX)) GO TO 830
      IMAX = JP
      AMAX = T
  830 JP = JP - 1
      GO TO 820
  840 JP = IMAX
  850 IF (.NOT.(W(J,NSOLN).NE.ZERO)) GO TO 860
      CALL DROTMG(SCALE(JP), SCALE(J), W(JP,NSOLN), W(J,NSOLN), SPARAM)
      W(J,NSOLN) = ZERO
      CALL DROTM(NP1-NSOLN, W(JP,NSP1), MDW, W(J,NSP1), MDW, SPARAM)
  860 J = JM1
      GO TO 810
C
C     SOLVE FOR Z(NSOLN)=PROPOSED NEW VALUE FOR X(NSOLN).
C     TEST IF THIS IS NONPOSITIVE OR TOO LARGE.
C     IF THIS WAS TRUE OR IF THE PIVOT TERM WAS ZERO REJECT
C     THE COL AS DEPENDENT.
  870 IF (.NOT.(W(NIV,NSOLN).NE.ZERO)) GO TO 890
      ISOL = NIV
      ASSIGN 880 TO IGO897
      GO TO 980
C
C     TEST-PROPOSED-NEW-COMPONENT
  880 GO TO 940
  890 IF (.NOT.(NIV.LE.ME .AND. W(MEP1,NSOLN).NE.ZERO)) GO TO 920
C
C     TRY TO ADD ROW MEP1 AS AN ADDITIONAL EQUALITY CONSTRAINT.
C     CHECK SIZE OF PROPOSED NEW SOLN COMPONENT.
C     REJECT IT IF IT IS TOO LARGE.
      ISOL = MEP1
      ASSIGN 900 TO IGO897
      GO TO 980
C
C     TEST-PROPOSED-NEW-COMPONENT
  900 IF (.NOT.(POS)) GO TO 910
C
C     SWAP ROWS MEP1 AND NIV, AND SCALE FACTORS FOR THESE ROWS.
      CALL DSWAP(NP1, W(MEP1,1), MDW, W(NIV,1), MDW)
      CALL DSWAP(1, SCALE(MEP1), 1, SCALE(NIV), 1)
      ITEMP = ITYPE(MEP1)
      ITYPE(MEP1) = ITYPE(NIV)
      ITYPE(NIV) = ITEMP
      ME = MEP1
      MEP1 = ME + 1
  910 GO TO 930
  920 POS = .FALSE.
  930 CONTINUE
  940 IF (POS) GO TO 950
      NSP1 = NSOLN
      NSOLN = NSOLN - 1
      NIV1 = NIV
      NIV = NIV - 1
  950 CONTINUE
  960 GO TO 740
  970 GO TO 1250
  980 CONTINUE
C
C     TO TEST-PROPOSED-NEW-COMPONENT
      Z2 = W(ISOL,NP1)/W(ISOL,NSOLN)
      Z(NSOLN) = Z2
      POS = Z2.GT.ZERO
      IF (.NOT.(Z2*EANORM.GE.BNORM .AND. POS)) GO TO 990
      POS = .NOT.(BLOWUP*Z2*EANORM.GE.BNORM)
  990 GO TO 1280
 1000 CONTINUE
C     TO COMPUTE-FINAL-SOLUTION
C
C     SOLVE SYSTEM, STORE RESULTS IN X(*).
C
      ASSIGN 1010 TO IGO958
      GO TO 1110
C     SOLVE-SYSTEM
 1010 CALL DCOPY(NSOLN, Z, 1, X, 1)
C
C     APPLY HOUSEHOLDER TRANSFORMATIONS TO X(*) IF KRANK.LT.L
      IF (.NOT.(0.LT.KRANK .AND. KRANK.LT.L)) GO TO 1030
      DO 1020 I=1,KRANK
        CALL H12(2, I, KRP1, L, W(I,1), MDW, H(I), X, 1, 1, 1)
 1020 CONTINUE
C
C     FILL IN TRAILING ZEROES FOR CONSTRAINED VARIABLES NOT IN SOLN.
 1030 IF (.NOT.(NSOLN.LT.N)) GO TO 1040
      X(NSP1) = ZERO
      CALL DCOPY(N-NSOLN, X(NSP1), 0, X(NSP1), 1)
C
C     REPERMUTE SOLN VECTOR TO NATURAL ORDER.
 1040 DO 1070 I=1,N
        J = I
 1050   IF (IPIVOT(J).EQ.I) GO TO 1060
        J = J + 1
        GO TO 1050
 1060   IPIVOT(J) = IPIVOT(I)
        IPIVOT(I) = J
        CALL DSWAP(1, X(J), 1, X(I), 1)
 1070 CONTINUE
C
C     RESCALE THE SOLN USING THE COL SCALING.
      DO 1080 J=1,N
        X(J) = X(J)*D(J)
 1080 CONTINUE
C     IF (.NOT.(NSOLN.LT.M)) GO TO 1100                           REMK
C     DO 1090 I=NSP1,M                                            REMK
      IF (.NOT.(NIV.LT.M)) GO TO 1100
      DO 1090 I = NIV1,M
        T = W(I,NP1)
        IF (I.LE.ME) T = T/ALAMDA
        T = (SCALE(I)*T)*T
        RNORM = RNORM + T
 1090 CONTINUE
 1100 RNORM = DSQRT(RNORM)
      GO TO 1210
C
C     TO SOLVE-SYSTEM
C
 1110 CONTINUE
      IF (.NOT.(DONE)) GO TO 1120
      ISOL = 1
      GO TO 1130
 1120 ISOL = LP1
 1130 IF (.NOT.(NSOLN.GE.ISOL)) GO TO 1190
C
C     COPY RT. HAND SIDE INTO TEMP VECTOR TO USE OVERWRITING METHOD.
      CALL DCOPY(NIV, W(1,NP1), 1, TEMP, 1)
      DO 1180 JJ=ISOL,NSOLN
        J = NSOLN - JJ + ISOL
        IF (.NOT.(J.GT.KRANK)) GO TO 1140
        I = NIV - JJ + ISOL
        GO TO 1150
 1140   I = J
 1150   IF (.NOT.(J.GT.KRANK .AND. J.LE.L)) GO TO 1160
        Z(J) = ZERO
        GO TO 1170
 1160   Z(J) = TEMP(I)/W(I,J)
        CALL DAXPY(I-1, -Z(J), W(1,J), 1, TEMP, 1)
 1170   CONTINUE
 1180 CONTINUE
 1190 GO TO 1270
 1200 GO TO IGO986, (40)
 1210 GO TO IGO980, (90)
 1220 GO TO IGO991, (30)
 1230 GO TO IGO998, (10)
 1240 GO TO IGO995, (20)
 1250 GO TO IGO983, (60)
 1260 GO TO IGO977, (220)
 1270 GO TO IGO958, (310, 1010)
 1280 GO TO IGO897, (880, 900)
 1290 GO TO IGO938, (460, 480, 500, 520)
      END
      SUBROUTINE WNLIT(W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM, WNL   10
     * IDOPE, DOPE, DONE)
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (BEGIN CHANGES AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/DCOPY/DCOPY/,/DROTM/DROTM/,
C     /DSCAL/DSCAL/,/DSQRT/DDSQRT/,                                  REMK
C     /DSWAP/DSWAP/,/DMAX1/DMAX1/,/IDAMAX/IDAMAX/,/.D-/.D-/,/D0/D0/
C
C     THIS IS A COMPANION SUBPROGRAM TO WNNLS( ).
C     THE DOCUMENTATION FOR WNNLS( ) HAS MORE COMPLETE
C     USAGE INSTRUCTIONS.
C
C     NOTE  THE M BY (N+1) MATRIX W( , ) CONTAINS THE RT. HAND SIDE
C           B AS THE (N+1)ST COL.
C
C
C     TRIANGULARIZE L1 BY L1 SUBSYSTEM, WHERE L1=MIN(M,L), WITH
C     COL INTERCHANGES.
C     REVISED MARCH 4, 1982
C
C++
      DOUBLE PRECISION W(MDW,1), H(1), SCALE(1), DOPE(4), SPARAM(5)
C     DOUBLE PRECISION ALSQ, AMAX, EANORM, FAC, FACTOR, HBAR, ONE, RN
      DOUBLE PRECISION ALSQ, EANORM, FACTOR, HBAR, ONE, RN
      DOUBLE PRECISION RNORM, SN, T, TAU, TENM3, ZERO
      DOUBLE PRECISION DMAX1
      INTEGER ITYPE(1), IPIVOT(1), IDOPE(8)
      INTEGER IDAMAX
      LOGICAL INDEP, DONE, RECALC
      DATA TENM3 /1.D-3/, ZERO /0.D0/, ONE /1.D0/
C
      ME = IDOPE(1)
      MEP1 = IDOPE(2)
      KRANK = IDOPE(3)
      KRP1 = IDOPE(4)
      NSOLN = IDOPE(5)
      NIV = IDOPE(6)
      NIV1 = IDOPE(7)
      L1 = IDOPE(8)
C
      ALSQ = DOPE(1)
      EANORM = DOPE(2)
C     FAC = DOPE(3)                                               REMK
      TAU = DOPE(4)
      NP1 = N + 1
      LB = MIN0(M-1,L)
      RECALC = .TRUE.
      RNORM = ZERO
      KRANK = 0
C     WE SET FACTOR=1.D0 SO THAT THE HEAVY WEIGHT ALAMDA WILL BE
C     INCLUDED IN THE TEST FOR COL INDEPENDENCE.
      FACTOR = 1.D0
      I = 1
      IP1 = 2
      LEND = L
   10 IF (.NOT.(I.LE.LB)) GO TO 150
      IF (.NOT.(I.LE.ME)) GO TO 130
C
C     SET IR TO POINT TO THE I-TH ROW.
      IR = I
      MEND = M
      ASSIGN 20 TO IGO996
      GO TO 460
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
   20 ASSIGN 30 TO IGO993
      GO TO 560
C
C     PERFORM-COL-INTERCHANGE
C
C     SET IC TO POINT TO I-TH COL.
   30 IC = I
      ASSIGN 40 TO IGO990
      GO TO 520
C
C     TEST-INDEP-OF-INCOMING-COL
   40 IF (.NOT.(INDEP)) GO TO 110
C
C     ELIMINATE I-TH COL BELOW DIAG. USING MOD. GIVENS TRANSFORMATIONS
C     APPLIED TO (A B).
      J = M
      DO 100 JJ=IP1,M
        JM1 = J - 1
        JP = JM1
C     WHEN OPERATING NEAR THE ME LINE, USE THE LARGEST ELT.        REMK
C     ABOVE IT AS THE PIVOT.                                       REMK
C       IF (.NOT.(J.EQ.MEP1)) GO TO 80                             REMK
C       IMAX = ME                                                  REMK
C       AMAX = SCALE(ME)*W(ME,I)**2                                REMK
C  50   IF (.NOT.(JP.GE.I)) GO TO 70                               REMK
C       T = SCALE(JP)*W(JP,I)**2                                   REMK
C       IF (.NOT.(T.GT.AMAX)) GO TO 60                             REMK
C       IMAX = JP                                                  REMK
C       AMAX = T                                                   REMK
C  60   JP = JP - 1                                                REMK
C       GO TO 50                                                   REMK
C  70   JP = IMAX                                                  REMK
        IF (.NOT. (JJ.EQ.M)) GO TO 70
        IF (.NOT. (I.LT.MEP1)) GO TO 80
        J = MEP1
        JP = I
        T = SCALE(JP)*W(JP,I)**2*TAU**2
        IF (.NOT.(T.GT.SCALE(J)*W(J,I)**2)) GO TO 130
        GO TO 80
   70   IF (.NOT.(J.EQ.MEP1)) GO TO 80
        J = JM1
        JM1 = J - 1
        JP = JM1
   80   IF (.NOT.(W(J,I).NE.ZERO)) GO TO 90
        CALL DROTMG(SCALE(JP), SCALE(J), W(JP,I), W(J,I), SPARAM)
        W(J,I) = ZERO
        CALL DROTM(NP1-I, W(JP,IP1), MDW, W(J,IP1), MDW, SPARAM)
   90   J = JM1
  100 CONTINUE
      GO TO 140
  110 CONTINUE
      IF (.NOT.(LEND.GT.I)) GO TO 130
C
C     COL I IS DEPENDENT. SWAP WITH COL LEND.
      MAX = LEND
C
C     PERFORM-COL-INTERCHANGE
      ASSIGN 120 TO IGO993
      GO TO 560
  120 CONTINUE
      LEND = LEND - 1
C
C     FIND COL IN REMAINING SET WITH LARGEST SS.
      MAX = IDAMAX(LEND-I+1,H(I),1) + I - 1
      HBAR = H(MAX)
      GO TO 30
  130 CONTINUE
      KRANK = I - 1
      GO TO 160
  140 I = IP1
      IP1 = IP1 + 1
      GO TO 10
  150 KRANK = L1
  160 CONTINUE
      KRP1 = KRANK + 1
      IF (.NOT.(KRANK.LT.ME)) GO TO 290
      FACTOR = ALSQ
      DO 170 I=KRP1,ME
        IF (L.GT.0) W(I,1) = ZERO
        CALL DCOPY(L, W(I,1), 0, W(I,1), MDW)
  170 CONTINUE
C
C     DETERMINE THE RANK OF THE REMAINING EQUALITY CONSTRAINT
C     EQUATIONS BY ELIMINATING WITHIN THE BLOCK OF CONSTRAINED
C     VARIABLES.  REMOVE ANY REDUNDANT CONSTRAINTS.
      IR = KRP1
      IF (.NOT. (L.LT.N)) GO TO 245
      LP1 = L + 1
      RECALC = .TRUE.
      LB = MIN0(L+ME-KRANK,N)
      I = LP1
      IP1 = I + 1
  180 IF (.NOT.(I.LE.LB)) GO TO 280
      IR = KRANK + I - L
      LEND = N
      MEND = ME
      ASSIGN 190 TO IGO996
      GO TO 460
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
  190 ASSIGN 200 TO IGO993
      GO TO 560
C
C     PERFORM-COL-INTERCHANGE
C
C     ELIMINATE ELEMENTS IN THE I-TH COL.
  200 J = ME
  210 IF (.NOT.(J.GT.IR)) GO TO 230
      JM1 = J - 1
      IF (.NOT.(W(J,I).NE.ZERO)) GO TO 220
      CALL DROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM)
      W(J,I) = ZERO
      CALL DROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM)
  220 J = JM1
      GO TO 210
C
C     SET IC=I=COL BEING ELIMINATED
  230 IC = I
      ASSIGN 240 TO IGO990
      GO TO 520
C
C     TEST-INDEP-OF-INCOMING-COL
  240 IF (INDEP) GO TO 270
C
C     REMOVE ANY REDUNDANT OR DEPENDENT EQUALITY CONSTRAINTS.
  245 CONTINUE
      JJ = IR
  250 IF (.NOT.(IR.LE.ME)) GO TO 260
      W(IR,1) = ZERO
      CALL DCOPY(N, W(IR,1), 0, W(IR,1), MDW)
      RNORM = RNORM + (SCALE(IR)*W(IR,NP1)/ALSQ)*W(IR,NP1)
      W(IR,NP1) = ZERO
      SCALE(IR) = ONE
C     RECLASSIFY THE ZEROED ROW AS A LEAST SQUARES EQUATION.
      ITYPE(IR) = 1
      IR = IR + 1
      GO TO 250
C
C     REDUCE ME TO REFLECT ANY DISCOVERED DEPENDENT EQUALITY
C     CONSTRAINTS.
  260 CONTINUE
      ME = JJ - 1
      MEP1 = ME + 1
      GO TO 300
  270 I = IP1
      IP1 = IP1 + 1
      GO TO 180
  280 CONTINUE
  290 CONTINUE
  300 CONTINUE
      IF (.NOT.(KRANK.LT.L1)) GO TO 420
C
C     TRY TO DETERMINE THE VARIABLES KRANK+1 THROUGH L1 FROM THE
C     LEAST SQUARES EQUATIONS.  CONTINUE THE TRIANGULARIZATION WITH
C     PIVOT ELEMENT W(MEP1,I).
C
      RECALC = .TRUE.
C
C     SET FACTOR=ALSQ TO REMOVE EFFECT OF HEAVY WEIGHT FROM
C     TEST FOR COL INDEPENDENCE.
      FACTOR = ALSQ
      KK = KRP1
      I = KK
      IP1 = I + 1
  310 IF (.NOT.(I.LE.L1)) GO TO 410
C
C     SET IR TO POINT TO THE MEP1-ST ROW.
      IR = MEP1
      LEND = L
      MEND = M
      ASSIGN 320 TO IGO996
      GO TO 460
C
C     UPDATE-COL-SS-AND-FIND-PIVOT-COL
  320 ASSIGN 330 TO IGO993
      GO TO 560
C
C     PERFORM-COL-INTERCHANGE
C
C     ELIMINATE I-TH COL BELOW THE IR-TH ELEMENT.
  330 IRP1 = IR + 1
      IF (.NOT.(IRP1.LE.M)) GO TO 355
      J = M
      DO 350 JJ=IRP1,M
        JM1 = J - 1
        IF (.NOT.(W(J,I).NE.ZERO)) GO TO 340
        CALL DROTMG(SCALE(JM1), SCALE(J), W(JM1,I), W(J,I), SPARAM)
        W(J,I) = ZERO
        CALL DROTM(NP1-I, W(JM1,IP1), MDW, W(J,IP1), MDW, SPARAM)
  340   J = JM1
  350 CONTINUE
  355 CONTINUE
C
C     TEST IF NEW PIVOT ELEMENT IS NEAR ZERO. IF SO, THE COL IS
C     DEPENDENT.
      T = SCALE(IR)*W(IR,I)**2
      INDEP = T.GT.TAU**2*EANORM**2
      IF (.NOT.INDEP) GO TO 380
C
C     COL TEST PASSED. NOW MUST PASS ROW NORM TEST TO BE CLASSIFIED
C     AS INDEPENDENT.
      RN = ZERO
      DO 370 I1=IR,M
        DO 360 J1=IP1,N
          RN = DMAX1(RN,SCALE(I1)*W(I1,J1)**2)
  360   CONTINUE
  370 CONTINUE
      INDEP = T.GT.TAU**2*RN
C
C     IF INDEPENDENT, SWAP THE IR-TH AND KRP1-ST ROWS TO MAINTAIN THE
C     TRIANGULAR FORM.  UPDATE THE RANK INDICATOR KRANK AND THE
C     EQUALITY CONSTRAINT POINTER ME.
  380 IF (.NOT.(INDEP)) GO TO 390
      CALL DSWAP(NP1, W(KRP1,1), MDW, W(IR,1), MDW)
      CALL DSWAP(1, SCALE(KRP1), 1, SCALE(IR), 1)
C     RECLASSIFY THE LEAST SQ. EQUATION AS AN EQUALITY CONSTRAINT AND
C     RESCALE IT.
      ITYPE(IR) = 0
      T = DSQRT(SCALE(KRP1))
      CALL DSCAL(NP1, T, W(KRP1,1), MDW)
      SCALE(KRP1) = ALSQ
      ME = MEP1
      MEP1 = ME + 1
      KRANK = KRP1
      KRP1 = KRANK + 1
      GO TO 400
  390 GO TO 430
  400 I = IP1
      IP1 = IP1 + 1
      GO TO 310
  410 CONTINUE
  420 CONTINUE
  430 CONTINUE
C
C     IF PSEUDORANK IS LESS THAN L, APPLY HOUSEHOLDER TRANS.
C     FROM RIGHT.
      IF (.NOT.(KRANK.LT.L)) GO TO 450
      DO 440 I=1,KRANK
        J = KRP1 - I
        CALL H12(1, J, KRP1, L, W(J,1), MDW, H(J), W, MDW, 1, J-1)
  440 CONTINUE
  450 NIV = KRANK + NSOLN - L
      NIV1 = NIV + 1
      IF (L.EQ.N) DONE = .TRUE.
C
C  END OF INITIAL TRIANGULARIZATION.
      IDOPE(1) = ME
      IDOPE(2) = MEP1
      IDOPE(3) = KRANK
      IDOPE(4) = KRP1
      IDOPE(5) = NSOLN
      IDOPE(6) = NIV
      IDOPE(7) = NIV1
      IDOPE(8) = L1
      RETURN
  460 CONTINUE
C
C     TO UPDATE-COL-SS-AND-FIND-PIVOT-COL
C
C     THE COL SS VECTOR WILL BE UPDATED AT EACH STEP. WHEN
C     NUMERICALLY NECESSARY, THESE VALUES WILL BE RECOMPUTED.
C
      IF (.NOT.(IR.NE.1 .AND. (.NOT.RECALC))) GO TO 480
C     UPDATE COL SS =SUM OF SQUARES.
      DO 470 J=I,LEND
        H(J) = H(J) - SCALE(IR-1)*W(IR-1,J)**2
  470 CONTINUE
C
C     TEST FOR NUMERICAL ACCURACY.
      MAX = IDAMAX(LEND-I+1,H(I),1) + I - 1
      RECALC = HBAR + TENM3*H(MAX).EQ.HBAR
C
C     IF REQUIRED, RECALCULATE COL SS, USING ROWS IR THROUGH MEND.
  480 IF (.NOT.(RECALC)) GO TO 510
      DO 500 J=I,LEND
        H(J) = ZERO
        DO 490 K=IR,MEND
          H(J) = H(J) + SCALE(K)*W(K,J)**2
  490   CONTINUE
  500 CONTINUE
C
C     FIND COL WITH LARGEST SS.
      MAX = IDAMAX(LEND-I+1,H(I),1) + I - 1
      HBAR = H(MAX)
  510 GO TO 600
  520 CONTINUE
C
C     TO TEST-INDEP-OF-INCOMING-COL
C
C     TEST THE COL IC TO DETERMINE IF IT IS LINEARLY INDEPENDENT
C     OF THE COLS ALREADY IN THE BASIS.  IN THE INIT TRI
C     STEP, WE USUALLY WANT THE HEAVY WEIGHT ALAMDA TO
C     BE INCLUDED IN THE TEST FOR INDEPENDENCE.  IN THIS CASE THE
C     VALUE OF FACTOR WILL HAVE BEEN SET TO 1.D0 BEFORE THIS
C     PROCEDURE IS INVOKED.  IN THE POTENTIALLY RANK DEFICIENT
C     PROBLEM, THE VALUE OF FACTOR WILL HAVE BEEN
C     SET TO ALSQ=ALAMDA**2 TO REMOVE THE EFFECT OF THE HEAVY WEIGHT
C     FROM THE TEST FOR INDEPENDENCE.
C
C     WRITE NEW COL AS PARTITIONED VECTOR
C             (A1)  NUMBER OF COMPONENTS IN SOLN SO FAR = NIV
C             (A2)  M-NIV COMPONENTS
C     AND COMPUTE  SN = INVERSE WEIGHTED LENGTH OF A1
C                  RN = INVERSE WEIGHTED LENGTH OF A2
C     CALL THE COL INDEPENDENT WHEN RN .GT. TAU*SN
      SN = ZERO
      RN = ZERO
      DO 550 J=1,MEND
        T = SCALE(J)
        IF (J.LE.ME) T = T/FACTOR
        T = T*W(J,IC)**2
        IF (.NOT.(J.LT.IR)) GO TO 530
        SN = SN + T
        GO TO 540
  530   RN = RN + T
  540   CONTINUE
  550 CONTINUE
      INDEP = RN.GT.TAU**2*SN
      GO TO 590
  560 CONTINUE
C
C     TO PERFORM-COL-INTERCHANGE
C
      IF (.NOT.(MAX.NE.I)) GO TO 570
C     EXCHANGE ELEMENTS OF PERMUTED INDEX VECTOR AND PERFORM COL
C     INTERCHANGES.
      ITEMP = IPIVOT(I)
      IPIVOT(I) = IPIVOT(MAX)
      IPIVOT(MAX) = ITEMP
      CALL DSWAP(M, W(1,MAX), 1, W(1,I), 1)
      T = H(MAX)
      H(MAX) = H(I)
      H(I) = T
  570 GO TO 580
  580 GO TO IGO993, (30, 200, 330, 120)
  590 GO TO IGO990, (40, 240)
  600 GO TO IGO996, (20, 190, 320)
      END
      SUBROUTINE HFTI(A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H,   HFTI  10
     * G, IP)
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (BEGIN CHANGES AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/ DSQRT/ DDSQRT/,
C     /, DABS/, DABS/,/ABS(/DABS(/,/D0/D0/
C
C     DIMENSION A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N)
C
C     WRITTEN BY C. L. LAWSON AND R. J. HANSON.  FROM THE BOOK SOLVING
C     LEAST SQUARES PROBLEMS, PRENTICE-HALL, INC. (1974). FOR FURTHER
C     ALGORITHMIC DETAILS SEE ALGORITHM HFTI IN CHAPTER 14.
C
C     ABSTRACT
C
C     THIS SUBROUTINE SOLVES A LINEAR LEAST SQUARES PROBLEM OR A SET OF
C     LINEAR LEAST SQUARES PROBLEMS HAVING THE SAME MATRIX BUT DIFFERENT
C     RIGHT-SIDE VECTORS.  THE PROBLEM DATA CONSISTS OF AN M BY N MATRIX
C     A, AN M BY NB MATRIX B, AND AN DABSOLUTE TOLERANCE PARAMETER TAU
C     WHOSE USAGE IS DESCRIBED BELOW.  THE NB COLUMN VECTORS OF B
C     REPRESENT RIGHT-SIDE VECTORS FOR NB DISTINCT LINEAR LEAST SQUARES
C     PROBLEMS.
C
C     THIS SET OF PROBLEMS CAN ALSO BE WRITTEN AS THE MATRIX LEAST
C     SQUARES PROBLEM
C
C                       AX = B,
C
C     WHERE X IS THE N BY NB SOLUTION MATRIX.
C
C     NOTE THAT IF B IS THE M BY M IDENTITY MATRIX, THEN X WILL BE THE
C     PSEUDO-INVERSE OF A.
C
C     THIS SUBROUTINE FIRST TRANSFORMS THE AUGMENTED MATRIX (A B) TO A
C     MATRIX (R C) USING PREMULTIPLYING HOUSEHOLDER TRANSFORMATIONS WITH
C     COLUMN INTERCHANGES.  ALL SUBDIAGONAL ELEMENTS IN THE MATRIX R ARE
C     ZERO AND ITS DIAGONAL ELEMENTS SATISFY
C
C                       DABS(R(I,I)).GE.ABS(R(I+1,I+1)),
C
C                       I = 1,...,L-1, WHERE
C
C                       L = MIN(M,N).
C
C     THE SUBROUTINE WILL COMPUTE AN INTEGER, KRANK, EQUAL TO THE NUMBER
C     OF DIAGONAL TERMS OF R THAT EXCEED TAU IN MAGNITUDE.  THEN A
C     SOLUTION OF MINIMUM EUCLIDEAN LENGTH IS COMPUTED USING THE FIRST
C     KRANK ROWS OF (R C).
C
C     TO BE SPECIFIC WE SUGGEST THAT THE USER CONSIDER AN EASILY
C     COMPUTABLE MATRIX NORM, SUCH AS, THE MAXIMUM OF ALL COLUMN SUMS OF
C     MAGNITUDES.
C
C     NOW IF THE RELATIVE UNCERTAINTY OF B IS EPS, (NORM OF UNCERTAINTY/
C     NORM OF B), IT IS SUGGESTED THAT TAU BE SET APPROXIMATELY EQUAL TO
C     EPS*(NORM OF A).
C
C     THE USER MUST DIMENSION ALL ARRAYS APPEARING IN THE CALL LIST..
C     A(MDA,N),(B(MDB,NB) OR B(M)),RNORM(NB),H(N),G(N),IP(N).  THIS
C     PERMITS THE SOLUTION OF A RANGE OF PROBLEMS IN THE SAME ARRAY
C     SPACE.
C
C     THE ENTIRE SET OF PARAMETERS FOR HFTI ARE
C
C     INPUT..
C
C     A(*,*),MDA,M,N    THE ARRAY A(*,*) INITIALLY CONTAINS THE M BY N
C                       MATRIX A OF THE LEAST SQUARES PROBLEM AX = B.
C                       THE FIRST DIMENSIONING PARAMETER OF THE ARRAY
C                       A(*,*) IS MDA, WHICH MUST SATISFY MDA.GE.M
C                       EITHER M.GE.N OR M.LT.N IS PERMITTED.  THERE
C                       IS NO RESTRICTION ON THE RANK OF A.  THE
C                       CONDITION MDA.LT.M IS CONSIDERED AN ERROR.
C
C     B(*),MDB,NB       IF NB = 0 THE SUBROUTINE WILL PERFORM THE
C                       ORTHOGONAL DECOMPOSITION BUT WILL MAKE NO
C                       REFERENCES TO THE ARRAY B(*).  IF NB.GT.0
C                       THE ARRAY B(*) MUST INITIALLY CONTAIN THE M BY
C                       NB MATRIX B OF THE LEAST SQUARES PROBLEM AX =
C                       B.  IF NB.GE.2 THE ARRAY B(*) MUST BE DOUBLY
C                       SUBSCRIPTED WITH FIRST DIMENSIONING PARAMETER
C                       MDB.GE.MAX(M,N).  IF NB = 1 THE ARRAY B(*) MAY
C                       BE EITHER DOUBLY OR SINGLY SUBSCRIPTED.  IN
C                       THE LATTER CASE THE VALUE OF MDB IS ARBITRARY
C                       BUT IT SHOULD BE SET TO SOME VALID INTEGER
C                       VALUE SUCH AS MDB = M.
C
C                       THE CONDITION OF NB.GT.1.AND.MDB.LT. MAX(M,N)
C                       IS CONSIDERED AN ERROR.
C
C     TAU               DABSOLUTE TOLERANCE PARAMETER PROVIDED BY USER
C                       FOR PSEUDORANK DETERMINATION.
C
C     H(*),G(*),IP(*)   ARRAYS OF WORKING SPACE USED BY HFTI.
C
C     OUTPUT..
C
C     A(*,*)            THE CONTENTS OF THE ARRAY A(*,*) WILL BE
C                       MODIFIED BY THE SUBROUTINE.  THESE CONTENTS
C                       ARE NOT GENERALLY REQUIRED BY THE USER.
C
C     B(*)              ON RETURN THE ARRAY B(*) WILL CONTAIN THE N BY
C                       NB SOLUTION MATRIX X.
C
C     KRANK             SET BY THE SUBROUTINE TO INDICATE THE
C                       PSEUDORANK OF A.
C
C     RNORM(*)          ON RETURN, RNORM(J) WILL CONTAIN THE EUCLIDEAN
C                       NORM OF THE RESIDUAL VECTOR FOR THE PROBLEM
C                       DEFINED BY THE J-TH COLUMN VECTOR OF THE ARRAY
C                       B(*,*) FOR J = 1,...,NB.
C
C     H(*),G(*)         ON RETURN THESE ARRAYS RESPECTIVELY CONTAIN
C                       ELEMENTS OF THE PRE- AND POST-MULTIPLYING
C                       HOUSEHOLDER TRANSFORMATIONS USED TO COMPUTE
C                       THE MINIMUM EUCLIDEAN LENGTH SOLUTION.
C
C     IP(*)             ARRAY IN WHICH THE SUBROUTINE RECORDS INDICES
C                       DESCRIBING THE PERMUTATION OF COLUMN VECTORS.
C                       THE CONTENTS OF ARRAYS H(*),G(*) AND IP(*)
C                       ARE NOT GENERALLY REQUIRED BY THE USER.
C
C++
      DOUBLE PRECISION A(MDA,N), B(MDB,1), H(N), G(N), RNORM(NB), TAU
      DOUBLE PRECISION FACTOR, HMAX, SM1, ZERO, SM, TMP
      DOUBLE PRECISION DIFF, DSQRT, ABS
      INTEGER   IP(N)
      ZERO = 0.D0
      FACTOR = 0.001D0
C
      K = 0
      LDIAG = MIN0(M,N)
      IF (LDIAG.LE.0) GO TO 310
      IF (.NOT.MDA.LT.M) GO TO 10
      NERR = 2
      IOPT = 2
      CALL XERROR(31HHFTI MDA.LT.M.. PROBABLE ERROR., 31, NERR, IOPT)
      RETURN
   10 CONTINUE
C
      IF (.NOT.(NB.GT.1 .AND. MAX0(M,N).GT.MDB)) GO TO 20
      NERR = 2
      IOPT = 2
      CALL XERROR(49HHFTI MDB.LT.MAX(M,N).AND.NB.GT.1. PROBABLE ERROR.,
     * 49, NERR, IOPT)
      RETURN
   20 CONTINUE
C
      DO 100 J=1,LDIAG
        IF (J.EQ.1) GO TO 40
C
C     UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX
C    ..
        LMAX = J
        DO 30 L=J,N
          H(L) = H(L) - A(J-1,L)**2
          IF (H(L).GT.H(LMAX)) LMAX = L
   30   CONTINUE
        IF (DIFF(HMAX+FACTOR*H(LMAX),HMAX)) 40, 40, 70
C
C     COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX
C    ..
   40   LMAX = J
        DO 60 L=J,N
          H(L) = ZERO
          DO 50 I=J,M
            H(L) = H(L) + A(I,L)**2
   50     CONTINUE
          IF (H(L).GT.H(LMAX)) LMAX = L
   60   CONTINUE
        HMAX = H(LMAX)
C    ..
C     LMAX HAS BEEN DETERMINED
C
C     DO COLUMN INTERCHANGES IF NEEDED.
C    ..
   70   CONTINUE
        IP(J) = LMAX
        IF (IP(J).EQ.J) GO TO 90
        DO 80 I=1,M
          TMP = A(I,J)
          A(I,J) = A(I,LMAX)
          A(I,LMAX) = TMP
   80   CONTINUE
        H(LMAX) = H(J)
   90 JCOL = MIN0(J+1,N)
C
C     COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B.
C    ..
        CALL H12(1, J, J+1, M, A(1,J), 1, H(J), A(1,JCOL), 1, MDA, N-J)
        CALL H12(2, J, J+1, M, A(1,J), 1, H(J), B, 1, MDB, NB)
  100 CONTINUE
C
C     DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU.
C    ..
      DO 110 J=1,LDIAG
        IF (DABS(A(J,J)).LE.TAU) GO TO 120
  110 CONTINUE
      K = LDIAG
      GO TO 130
  120 K = J - 1
  130 KP1 = K + 1
C
C     COMPUTE THE NORMS OF THE RESIDUAL VECTORS.
C
      IF (NB.LE.0) GO TO 170
      DO 160 JB=1,NB
        TMP = ZERO
        IF (KP1.GT.M) GO TO 150
        DO 140 I=KP1,M
          TMP = TMP + B(I,JB)**2
  140   CONTINUE
  150   RNORM(JB) = DSQRT(TMP)
  160 CONTINUE
  170 CONTINUE
C                                           SPECIAL FOR PSEUDORANK = 0
      IF (K.GT.0) GO TO 200
      IF (NB.LE.0) GO TO 310
      DO 190 JB=1,NB
        DO 180 I=1,N
          B(I,JB) = ZERO
  180   CONTINUE
  190 CONTINUE
      GO TO 310
C
C     IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER
C     DECOMPOSITION OF FIRST K ROWS.
C    ..
  200 IF (K.EQ.N) GO TO 220
      DO 210 II=1,K
        I = KP1 - II
        CALL H12(1, I, KP1, N, A(I,1), MDA, G(I), A, MDA, 1, I-1)
  210 CONTINUE
  220 CONTINUE
C
C
      IF (NB.LE.0) GO TO 310
      DO 300 JB=1,NB
C
C     SOLVE THE K BY K TRIANGULAR SYSTEM.
C    ..
        DO 250 L=1,K
          SM = ZERO
          I = KP1 - L
          IF (I.EQ.K) GO TO 240
          IP1 = I + 1
          DO 230 J=IP1,K
            SM = SM + A(I,J)*B(J,JB)
  230     CONTINUE
  240     SM1 = SM
          B(I,JB) = (B(I,JB)-SM1)/A(I,I)
  250   CONTINUE
C
C     COMPLETE COMPUTATION OF SOLUTION VECTOR.
C    ..
        IF (K.EQ.N) GO TO 280
        DO 260 J=KP1,N
          B(J,JB) = ZERO
  260   CONTINUE
        DO 270 I=1,K
          CALL H12(2, I, KP1, N, A(I,1), MDA, G(I), B(1,JB), 1, MDB, 1)
  270   CONTINUE
C
C      RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE
C      COLUMN INTERCHANGES.
C    ..
  280   DO 290 JJ=1,LDIAG
          J = LDIAG + 1 - JJ
          IF (IP(J).EQ.J) GO TO 290
          L = IP(J)
          TMP = B(L,JB)
          B(L,JB) = B(J,JB)
          B(J,JB) = TMP
  290   CONTINUE
  300 CONTINUE
C    ..
C     THE SOLUTION VECTORS, X, ARE NOW
C     IN THE FIRST  N  ROWS OF THE ARRAY B(,).
C
  310 KRANK = K
      RETURN
      END
      SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)          H 12  10
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (START CHANGES AT LINE WITH C++ IN COLS. 1-3.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/,/DDOT/DDOT/,/DABS,/DABS,/,
C     /DSWAP/DSWAP/,/DSQRT/DDSQRT/,/DABS(/ DABS(/,/DMAX1/DMAX1/,
C     /DAXPY/DAXPY/,/D0/D0/
C
C
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
C
C     MODIFIED AT SANDIA LABS., MAY 1977, TO --
C
C     1)  REMOVE DOUBLE PRECISION ACCUMULATION, AND
C     2)  INCLUDE USAGE OF THE BASIC LINEAR ALGEBRA PACKAGE FOR
C         VECTORS LONGER THAN A PARTICULAR THRESHOLD.
C
C     CONSTRUCTION AND/OR APPLICATION OF A SINGLE
C     HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B
C
C     MODE    = 1 OR 2   TO SELECT ALGORITHM  H1  OR  H2 .
C     LPIVOT IS THE INDEX OF THE PIVOT ELEMENT.
C     L1,M   IF L1 .LE. M   THE TRANSFORMATION WILL BE CONSTRUCTED TO
C            ZERO ELEMENTS INDEXED FROM L1 THROUGH M.   IF L1 GT. M
C            THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION.
C     U(),IUE,UP    ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR.
C                   IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS.
C                                       ON EXIT FROM H1 U() AND UP
C                   CONTAIN QUANTITIES DEFINING THE VECTOR U OF THE
C                   HOUSEHOLDER TRANSFORMATION.   ON ENTRY TO H2 U()
C                   AND UP SHOULD CONTAIN QUANTITIES PREVIOUSLY COMPUTED
C                   BY H1.  THESE WILL NOT BE MODIFIED BY H2.
C     C()    ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE
C            REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER
C            TRANSFORMATION IS TO BE APPLIED.  ON EXIT C() CONTAINS THE
C            SET OF TRANSFORMED VECTORS.
C     ICE    STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C().
C     ICV    STORAGE INCREMENT BETWEEN VECTORS IN C().
C     NCV    NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0
C            NO OPERATIONS WILL BE DONE ON C().
C
C     SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV)
C++
      DOUBLE PRECISION U(IUE,M), C(1), UP
      DOUBLE PRECISION B, CL, CLINV, ONE, SM, UL1M1
      DOUBLE PRECISION DABS, DMAX1, DSQRT, DDOT
      ONE=1.D0
C
      IF (0.GE.LPIVOT.OR.LPIVOT.GE.L1.OR.L1.GT.M) RETURN
      CL=DABS(U(1,LPIVOT))
      IF (MODE.EQ.2) GO TO 60
C                            ****** CONSTRUCT THE TRANSFORMATION. ******
          DO 10 J=L1,M
   10     CL=DMAX1(DABS(U(1,J)),CL)
      IF (CL) 130,130,20
   20 CLINV=ONE/CL
      SM=(U(1,LPIVOT)*CLINV)**2
          DO 30 J=L1,M
   30     SM=SM+(U(1,J)*CLINV)**2
      CL=CL*DSQRT(SM)
      IF (U(1,LPIVOT)) 50,50,40
   40 CL=-CL
   50 UP=U(1,LPIVOT)-CL
      U(1,LPIVOT)=CL
      GO TO 70
C            ****** APPLY THE TRANSFORMATION  I+U*(U**T)/B  TO C. ******
C
   60 IF (CL) 130,130,70
   70 IF (NCV.LE.0) RETURN
      B=UP*U(1,LPIVOT)
C                       B  MUST BE NONPOSITIVE HERE.  IF B = 0., RETURN.
C
      IF (B) 80,130,130
   80 B=ONE/B
      MML1P2=M-L1+2
      IF (MML1P2.GT.20) GO TO 140
      I2=1-ICV+ICE*(LPIVOT-1)
      INCR=ICE*(L1-LPIVOT)
          DO 120 J=1,NCV
          I2=I2+ICV
          I3=I2+INCR
          I4=I3
          SM=C(I2)*UP
              DO 90 I=L1,M
              SM=SM+C(I3)*U(1,I)
   90         I3=I3+ICE
          IF (SM) 100,120,100
  100     SM=SM*B
          C(I2)=C(I2)+SM*UP
              DO 110 I=L1,M
              C(I4)=C(I4)+SM*U(1,I)
  110         I4=I4+ICE
  120     CONTINUE
  130 RETURN
  140 CONTINUE
      L1M1=L1-1
      KL1=1+(L1M1-1)*ICE
      KL2=KL1
      KLP=1+(LPIVOT-1)*ICE
      UL1M1=U(1,L1M1)
      U(1,L1M1)=UP
      IF (LPIVOT.EQ.L1M1) GO TO 150
      CALL DSWAP(NCV,C(KL1),ICV,C(KLP),ICV)
  150 CONTINUE
          DO 160 J=1,NCV
          SM=DDOT(MML1P2,U(1,L1M1),IUE,C(KL1),ICE)
          SM=SM*B
          CALL DAXPY (MML1P2,SM,U(1,L1M1),IUE,C(KL1),ICE)
          KL1=KL1+ICV
  160 CONTINUE
      U(1,L1M1)=UL1M1
      IF (LPIVOT.EQ.L1M1) RETURN
      KL1=KL2
      CALL DSWAP(NCV,C(KL1),ICV,C(KLP),ICV)
      RETURN
      END
      DOUBLE PRECISION FUNCTION DIFF(X,Y)                               DIFF  10
C
C     THE EDITING REQUIRED TO CONVERT THIS SUBROUTINE FROM SINGLE TO
C     DOUBLE PRECISION INVOLVES THE FOLLOWING CHARACTER STRING CHANGES.
C     USE AN EDITING COMMAND (CHANGE) /STRING-1/(TO)STRING-2/.
C     (APPLY CHANGES TO ENTIRE PROGRAM UNIT.)
C     /REAL (12 BLANKS)/DOUBLE PRECISION/
C
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974
      DOUBLE PRECISION  X, Y
      DIFF=X-Y
      RETURN
      END
C     !!!sed-endpoint!!! sed-conversion to dconstr.f should end here

C     modification by w.h. 09/02: 
C     the following 1778 lines of the original file are now in blas.f
     
C modification by w.h. 5/98: 
C XERROR and XERRWV are the only routines called from above.
C we shorten them:
C just print out messg and exit.
C see hist/hist10/constr.f for previous version

      SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL)                        XER   10
C
C     (old) ABSTRACT
C        XERROR PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER
C        DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE
C        OF THE LIBRARY ERROR CONTROL FLAG, KONTRL.
C        (SEE SUBROUTINE XSETF FOR DETAILS.)
C
C     DESCRIPTION OF PARAMETERS
C      --INPUT--
C        MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED, CONTAINING
C                NO MORE THAN 72 CHARACTERS.
C        NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG.
C        further parameters ignored
C
      call myErr(MESSG,NMESSG)
      END

      SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2)      XER   10
C     modifikation 22.10.98
      double precision R1,R2
C
C     (old) ABSTRACT
C        XERRWV PROCESSES A DIAGNOSTIC MESSAGE, IN A MANNER
C        DETERMINED BY THE VALUE OF LEVEL AND THE CURRENT VALUE
C        OF THE LIBRARY ERROR CONTROL FLAG, KONTRL.
C        (SEE SUBROUTINE XSETF FOR DETAILS.)
C        IN ADDITION, UP TO TWO INTEGER VALUES AND TWO REAL
C        VALUES MAY BE PRINTED ALONG WITH THE MESSAGE.
C
C     (old) DESCRIPTION OF PARAMETERS
C      --INPUT--
C        MESSG - THE HOLLERITH MESSAGE TO BE PROCESSED.
C        NMESSG- THE ACTUAL NUMBER OF CHARACTERS IN MESSG.
C        further parameters ignored
C END OF ABSTRACT
C     modifikation 04.12.01: wrap to odessa's XERR, now in solvlin.cc

      call xerr(MESSG, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
C      call myErr(MESSG,NMESSG)
      END

      INTEGER FUNCTION I1MACH(I)
C***BEGIN PROLOGUE  I1MACH
C***REVISION DATE  811015   (YYMMDD)
C***CATEGORY NO.  Q
C***KEYWORDS  MACHINE CONSTANTS,INTEGER
C***DATE WRITTEN   1975
C***AUTHOR FOX P.A.,HALL A.D.,SCHRYER N.L. (BELL LABS)
C***PURPOSE
C RETURNS INTEGER MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   THESE MACHINE CONSTANT ROUTINES MUST BE ACTIVATED FOR
C   A PARTICULAR ENVIRONMENT.
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     I1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C     SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C     AS FOLLOWS, FOR EXAMPLE
C
C          K = I1MACH(I)
C
C     WHERE I=1,...,16.  THE (OUTPUT) VALUE OF K ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  I/O UNIT NUMBERS.
C    I1MACH( 1) = THE STANDARD INPUT UNIT.
C    I1MACH( 2) = THE STANDARD OUTPUT UNIT.
C    I1MACH( 3) = THE STANDARD PUNCH UNIT.
C    I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT.
C
C  WORDS.
C    I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT.
C    I1MACH( 6) = THE NUMBER OF CHARACTERS PER INTEGER STORAGE UNIT.
C
C  INTEGERS.
C    ASSUME INTEGERS ARE REPRESENTED IN THE S-DIGIT, BASE-A FORM
C
C               SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
C
C               WHERE 0 .LE. X(I) .LT. A FOR I=0,...,S-1.
C    I1MACH( 7) = A, THE BASE.
C    I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS.
C    I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE.
C
C  FLOATING-POINT NUMBERS.
C    ASSUME FLOATING-POINT NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C    BASE-B FORM
C               SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C               WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T,
C               0 .LT. X(1), AND EMIN .LE. E .LE. EMAX.
C    I1MACH(10) = B, THE BASE.
C
C  SINGLE-PRECISION
C    I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C    I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C    I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C  DOUBLE-PRECISION
C    I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C    I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C    I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C  REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C  I1MACH(1) - I1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C  WITH THE LOCAL OPERATING SYSTEM.
C
C***REFERENCES
C  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A PORTABLE LIBRARY*,
C  ACM TRANSACTION ON MATHEMATICAL SOFTWARE, VOL. 4, NO. 2,
C  JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  I1MACH
C
      INTEGER IMACH(16),OUTPUT
C
C     EQUIVALENCE (IMACH(4),OUTPUT)
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA IMACH( 1) /    7 /
C     DATA IMACH( 2) /    2 /
C     DATA IMACH( 3) /    2 /
C     DATA IMACH( 4) /    2 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    4 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   33 /
C     DATA IMACH( 9) / Z1FFFFFFFF /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   24 /
C     DATA IMACH(12) / -256 /
C     DATA IMACH(13) /  255 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) / -256 /
C     DATA IMACH(16) /  255 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA IMACH( 1) /   5 /
C     DATA IMACH( 2) /   6 /
C     DATA IMACH( 3) /   7 /
C     DATA IMACH( 4) /   6 /
C     DATA IMACH( 5) /  48 /
C     DATA IMACH( 6) /   6 /
C     DATA IMACH( 7) /   2 /
C     DATA IMACH( 8) /  39 /
C     DATA IMACH( 9) / O0007777777777777 /
C     DATA IMACH(10) /   8 /
C     DATA IMACH(11) /  13 /
C     DATA IMACH(12) / -50 /
C     DATA IMACH(13) /  76 /
C     DATA IMACH(14) /  26 /
C     DATA IMACH(15) / -50 /
C     DATA IMACH(16) /  76 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA IMACH( 1) /   5 /
C     DATA IMACH( 2) /   6 /
C     DATA IMACH( 3) /   7 /
C     DATA IMACH( 4) /   6 /
C     DATA IMACH( 5) /  48 /
C     DATA IMACH( 6) /   6 /
C     DATA IMACH( 7) /   2 /
C     DATA IMACH( 8) /  39 /
C     DATA IMACH( 9) / O0007777777777777 /
C     DATA IMACH(10) /   8 /
C     DATA IMACH(11) /  13 /
C     DATA IMACH(12) / -50 /
C     DATA IMACH(13) /  76 /
C     DATA IMACH(14) /  26 /
C     DATA IMACH(15) / -32754 /
C     DATA IMACH(16) /  32780 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    7 /
C     DATA IMACH( 4) /6LOUTPUT/
C     DATA IMACH( 5) /   60 /
C     DATA IMACH( 6) /   10 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   48 /
C     DATA IMACH( 9) / 00007777777777777777B /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   48 /
C     DATA IMACH(12) / -974 /
C     DATA IMACH(13) / 1070 /
C     DATA IMACH(14) /   96 /
C     DATA IMACH(15) / -927 /
C     DATA IMACH(16) / 1070 /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA IMACH( 1) /   100 /
C     DATA IMACH( 2) /   101 /
C     DATA IMACH( 3) /   102 /
C     DATA IMACH( 4) /   101 /
C     DATA IMACH( 5) /    64 /
C     DATA IMACH( 6) /     8 /
C     DATA IMACH( 7) /     2 /
C     DATA IMACH( 8) /    63 /
C     DATA IMACH( 9) /  777777777777777777777B /
C     DATA IMACH(10) /     2 /
C     DATA IMACH(11) /    48 /
C     DATA IMACH(12) / -8192 /
C     DATA IMACH(13) /  8191 /
C     DATA IMACH(14) /    96 /
C     DATA IMACH(15) / -8192 /
C     DATA IMACH(16) /  8191 /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     DATA IMACH( 1) /   11 /
C     DATA IMACH( 2) /   12 /
C     DATA IMACH( 3) /    8 /
C     DATA IMACH( 4) /   10 /
C     DATA IMACH( 5) /   16 /
C     DATA IMACH( 6) /    2 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   15 /
C     DATA IMACH( 9) /32767 /
C     DATA IMACH(10) /   16 /
C     DATA IMACH(11) /    6 /
C     DATA IMACH(12) /  -64 /
C     DATA IMACH(13) /   63 /
C     DATA IMACH(14) /   14 /
C     DATA IMACH(15) /  -64 /
C     DATA IMACH(16) /   63 /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA IMACH( 1) /       5 /
C     DATA IMACH( 2) /       6 /
C     DATA IMACH( 3) /       0 /
C     DATA IMACH( 4) /       6 /
C     DATA IMACH( 5) /      24 /
C     DATA IMACH( 6) /       3 /
C     DATA IMACH( 7) /       2 /
C     DATA IMACH( 8) /      23 /
C     DATA IMACH( 9) / 8388607 /
C     DATA IMACH(10) /       2 /
C     DATA IMACH(11) /      23 /
C     DATA IMACH(12) /    -127 /
C     DATA IMACH(13) /     127 /
C     DATA IMACH(14) /      38 /
C     DATA IMACH(15) /    -127 /
C     DATA IMACH(16) /     127 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /   43 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    6 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -127 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   63 /
C     DATA IMACH(15) / -127 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA IMACH(1) /      5/
C     DATA IMACH(2) /      6 /
C     DATA IMACH(3) /      4 /
C     DATA IMACH(4) /      1 /
C     DATA IMACH(5) /     16 /
C     DATA IMACH(6) /      2 /
C     DATA IMACH(7) /      2 /
C     DATA IMACH(8) /     15 /
C     DATA IMACH(9) /  32767 /
C     DATA IMACH(10)/      2 /
C     DATA IMACH(11)/     23 /
C     DATA IMACH(12)/   -128 /
C     DATA IMACH(13)/    127 /
C     DATA IMACH(14)/     39 /
C     DATA IMACH(15)/   -128 /
C     DATA IMACH(16)/    127 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA IMACH(1) /      5 /
C     DATA IMACH(2) /      6 /
C     DATA IMACH(3) /      4 /
C     DATA IMACH(4) /      1 /
C     DATA IMACH(5) /     16 /
C     DATA IMACH(6) /      2 /
C     DATA IMACH(7) /      2 /
C     DATA IMACH(8) /     15 /
C     DATA IMACH(9) /  32767 /
C     DATA IMACH(10)/      2 /
C     DATA IMACH(11)/     23 /
C     DATA IMACH(12)/   -128 /
C     DATA IMACH(13)/    127 /
C     DATA IMACH(14)/     55 /
C     DATA IMACH(15)/   -128 /
C     DATA IMACH(16)/    127 /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
CCC
CCC assuming IBM 360 will work on IBM RS76000
CCC
C
      DATA IMACH( 1) /   5 /
      DATA IMACH( 2) /   6 /
      DATA IMACH( 3) /   7 /
      DATA IMACH( 4) /   6 /
      DATA IMACH( 5) /  32 /
      DATA IMACH( 6) /   4 /
      DATA IMACH( 7) /   2 /
      DATA IMACH( 8) /  31 /
      DATA IMACH( 9) / 2147483647 /
CCC was Z7FFFFFFF /
      DATA IMACH(10) /  16 /
      DATA IMACH(11) /   6 /
      DATA IMACH(12) / -64 /
      DATA IMACH(13) /  63 /
      DATA IMACH(14) /  14 /
      DATA IMACH(15) / -64 /
      DATA IMACH(16) /  63 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    5 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / "377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   54 /
C     DATA IMACH(15) / -101 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    5 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / "377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   62 /
C     DATA IMACH(15) / -128 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN S SUPPORTING
C     32-BIT INTEGER ARITHMETIC.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   32 /
C     DATA IMACH( 6) /    4 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   31 /
C     DATA IMACH( 9) / 2147483647 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   24 /
C     DATA IMACH(12) / -127 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   56 /
C     DATA IMACH(15) / -127 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN S SUPPORTING
C     16-BIT INTEGER ARITHMETIC.
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    5 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   16 /
C     DATA IMACH( 6) /    2 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   15 /
C     DATA IMACH( 9) / 32767 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   24 /
C     DATA IMACH(12) / -127 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   56 /
C     DATA IMACH(15) / -127 /
C     DATA IMACH(16) /  127 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    1 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    4 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) /-1024 /
C     DATA IMACH(16) / 1023 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
C
C     DATA IMACH( 1) /    5 /
C     DATA IMACH( 2) /    6 /
C     DATA IMACH( 3) /    7 /
C     DATA IMACH( 4) /    6 /
C     DATA IMACH( 5) /   36 /
C     DATA IMACH( 6) /    6 /
C     DATA IMACH( 7) /    2 /
C     DATA IMACH( 8) /   35 /
C     DATA IMACH( 9) / O377777777777 /
C     DATA IMACH(10) /    2 /
C     DATA IMACH(11) /   27 /
C     DATA IMACH(12) / -128 /
C     DATA IMACH(13) /  127 /
C     DATA IMACH(14) /   60 /
C     DATA IMACH(15) /-1024/
C     DATA IMACH(16) / 1023 /
C
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C
C     DATA IMACH(1) /    5 /
C     DATA IMACH(2) /    6 /
C     DATA IMACH(3) /    5 /
C     DATA IMACH(4) /    6 /
C     DATA IMACH(5) /   32 /
C     DATA IMACH(6) /    4 /
C     DATA IMACH(7) /    2 /
C     DATA IMACH(8) /   31 /
C     DATA IMACH(9) /2147483647 /
C     DATA IMACH(10)/    2 /
C     DATA IMACH(11)/   24 /
C     DATA IMACH(12)/ -127 /
C     DATA IMACH(13)/  127 /
C     DATA IMACH(14)/   56 /
C     DATA IMACH(15)/ -127 /
C     DATA IMACH(16)/  127 /
C
C***FIRST EXECUTABLE STATEMENT  I1MACH
C
      IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
C
      I1MACH=IMACH(I)
      RETURN
C
   10 CONTINUE
C modification by w.h. 5/98
      CALL myErr(39H1ERROR    1 IN I1MACH - I OUT OF BOUNDS, 39)
      END

C     main program: s. constr.main
      REAL FUNCTION RAN(K)                                              RAN   10
C
C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266
C      BY PIKE AND HILL (MODIFIED BY HANSSON)
C      COLLECTED ALG. FROM CACM.
C
C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
C      29 BITS.
C
      INTEGER IY,K
      DATA IY/100001/
C
      IF(K.GT.0) IY = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      RAN = FLOAT(IY) / 2796203.0E0
      RETURN
C     ---------- LAST CARD OF RAN ----------
      END
CCC
CCC the following lines were not commented out in the original version
CCC
C
C    3    1    2    0    2                                               DATA  10
C   1.     1000.     1000.                                               DATA  20
C    3    1    2    1    2                                               DATA  30
C   1.     1000.     1000.                                               DATA  40
C    3    1    2    2    2                                               DATA  50
C1000.     1000.     1000.                                               DATA  60
C    4    2    2    2    3                                               DATA  70
C1000.     1000.     1000.                                               DATA  80
C    5    3    2    2    4                                               DATA  90
C10000.    1000.     1000.                                               DATA 100
C    0    0    0    0   -1                                               DATA 110
