C     ALGORITHM 615, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 2,
C     JUN., 1984, P. 202-206.
C      PROGRAM SUBL1 (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
C****************************************************************
C      *THIS CODE SOLVES THE L1 BEST SUBSET PROBLEM*
C      *SUBMITTED TO ACM TRANS.SPRING 1982*
C       *LAST REVISED MAR. 12, 1984*
C****************************************************************
      DIMENSION X(300,20), Y(300)
      DIMENSION ZL(20), BVAL(210), IDEX(210), ISTAT(20)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER  I,IBASE,IDEX,IFAULT,IK,IK1,IMT,INCOL,INDEX,INPROB,IPAR
      INTEGER  ISTAT,ITER,J,JMIN,K,L,LABEL,LEVEL,M,MININ,MMAX
      INTEGER  N,NAME,NMAX,NPROB
C
      DOUBLE PRECISION ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
C++   REAL             ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
      DOUBLE PRECISION SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
C++   REAL             SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
C
C
      LOGICAL DIRECT,INTL
      INTEGER UNITGO,UNITIN
C
C-----UNITIN    IS THE INPUT UNIT NUMBER
C-----UNITGO    IS THE OUTPUT UNIT NUMBER
C
      DIMENSION IMT(72), NAME(40)
C
C-----IMT       IS USED TO STORE THE INPUT FORMAT
C-----NAME      SAVES THE PROBLEM NAME
C
      UNITIN=5
      UNITGO=6
      NMAX=300
      MMAX=30
      NPROB=0
C
C-----K        STORES THE NUMBER OF PARAMETERS IN THE FULL MODEL
C
   10 READ (UNITIN,50,END=45) K,NAME
      WRITE(UNITGO,50)K,NAME
      IF (K.EQ.0) STOP
      WRITE (UNITGO,60) NAME
C
C-----N          STORES THE NUMBER OF OBSERVATIONS
C-----MININ      STORES THE MINIMUM NUMBER OF PARAMETERS CONSIDERED
C-----POPT       DEVIATION FROM OPTIMALITY ALLOWED
C
      READ (UNITIN,130) N,MININ,POPT
   51 FORMAT (4X,22HNUMBER OF OBSERVATION=,I5/4X,21HNUMBER OF PARAMETERS
     1=,I5/4X,40HMINIMUM NUMBER OF PARAMETERS CONSIDERED=,I5/4X,45HPERCE
     2NTAGE DEVIATION FROM OPTIMALITY ALLOWED=,F6.2/8X,2H**,
     3 23HBEST SUBSET LAV PROGRAM)
      WRITE (UNITGO,51) N,K,MININ,POPT
      READ (UNITIN,150) IMT
      DO 20 I=1,N
         READ (UNITIN,IMT) Y(I), (X(I,J),J=1,K)
C
C        WRITE (UNITGO,IMT) Y(I), (X(I,J),J=1,K)
C
   20 CONTINUE
C
      READ (UNITIN,70) (ISTAT(I),I=1,K)
C
      ITER=0
C
C******CALL THE ROUTINE TO FIND THE BEST SUBSETS
C
C
      CALL KBEST (X,Y,K,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL,IDEX,
     1ISTAT,ZL)
C
      WRITE (UNITGO,80) IFAULT
C
C     WRITE THE FINAL BEST SUBSET SOLUTION
C
      J=K*(K+1)/2
      JMIN=(MININ-1)*(MININ)/2
      J=J-JMIN
      DO 40 I=MININ,K
         WRITE (UNITGO,110) I
         WRITE (UNITGO,100) ZL(I)
         DO 30 L=1,I
            WRITE (UNITGO,90) IDEX(J),BVAL(J)
            J=J-1
   30    CONTINUE
   40 CONTINUE
C
          NPROB=NPROB+1
          WRITE (UNITGO,160)ITER
      GO TO 10
C
   45 STOP
C
   50 FORMAT (I5,3X,40A1)
   60 FORMAT (3X,5H     ,14HPROBLEM TITLE ,2X,40A1)
   70 FORMAT (10I2)
   80 FORMAT (8H IFAULT=,I3)
   90 FORMAT (3X/6H BETA(,I3,1H),F15.3)
  100 FORMAT (4X,14(1H*),23HSUM OF ABSOLUTE VALUES=,F15.3)
  110 FORMAT (//3X,36HBEST RESULTS FOR LAV SUBSET OF SIZE=,I3)
  130 FORMAT (I5,I2,F6.2)
  150 FORMAT (72A1)
  160 FORMAT (//5X,16HITERATION COUNT=,I7)
C
      END
      SUBROUTINE KBEST (X,Y,M,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL
     1,IDEX,ISTAT,ZL)
C
C***********************************************************************
C
C
C     THE PURPOSE OF THIS PROGRAM IS TO DETERMINE THE BEST SUBSET OF
C     PARAMETERS TO FIT A LINEAR REGRESSION UNDER AN LEAST ABSOLUTE
C     VALUE CRITERION. THIS PROGRAM UTILIZES THE SIMPLEX METHOD OF
C     LINEAR PROGRAMMING WITHIN A BRANCH-AND-BOUND ALGORITHM TO
C     SOLVE THE BEST SUBSET PROBLEM.
C
C
C     THE ALGORITHM IS BASED ON THE PUBLICATION:
C        ARMSTRONG,R.D. AND M.T. KUNG "AN ALGORITHM TO SELECT THE BEST
C        SUBSET FOR A LEAST ABSOLUTE VALUE REGRESSION PROBLEM,"
C        OPTIMIZATION IN STATISTICS, TIMS STUDIES OF THE MANAGEMENT
C        SCIENCES.
C
C   FORMAL PARAMETERS
C
C   X      REAL ARRAY        INPUT:  VALUES OF INDEPENDENT VARIABLES
C          (NMAX,MMAX)                SUCH THAT EACH ROW CORREPSONDS TO
C                                     AN OBSERVATION
C
C   Y      REAL ARRAY        INPUT:  VALUES OF THE DEPENDENT VARIABLES
C            (NMAX)
C
C   M      INTEGER           INPUT:  NUMBER OF DEPENDENT VARIABLES
C
C   N      INTEGER           INPUT:  NUMBER OF OBSERVATIONS
C
C   ITER   INTEGER           OUTPUT: NUMBER OF ITERATIONS
C
C   IFAULT INTEGER           OUTPUT: FAILURE INDICATOR
C                              =0     NORMAL TERMINATION
C                              =1     OBSERVATION MATRIX DOES NOT HAVE
C                                     FULL ROW RANK (RANK M)
C                              =2     PROBLEM SIZE OUT OF RANGE
C                              =3     NO PIVOT ELEMENT FOUND IMPLYING
C                                     NEAR SINGULAR BASIS
C
C   POPT    REAL             INPUT:  PERCENTAGE DEVIATION FROM
C                                    OPTIMALITY ALLOWED
C
C   MININ   INTEGER          INPUT:  MINIMUM NUMBER OF PARAMETERS IN THE
C                                    MODEL.  BEST SUBSET OF SIZE MININ T
C                                    M IS OBTAINED.
C
C   NMAX    INTEGER          INPUT:  DIMENSION OF ROWS IN X (ALSO Y)
C
C   MMAX    INTEGER          INPUT:  DIMENSION OF COLUMNS IN X
C
C   BVAL    REAL ARRAY       OUTPUT: ARRAY OF OPTIMAL BETA VALUES FOR
C                                    EACH SUBSET. THE BETA VALUES FOR
C                                    THE SUBSET OF SIZE M ARE STORED
C                                    IN POSITIONS BVAL(1),BVAL(2),...,
C                                    BVAL(M), FOR THE SUBSET OF SIZE
C                                    M-1 THE VALUES ARE STORED IN
C                                    POSITIONS BVAL(M+1),BVAL(M+2),...,
C                                    BVAL(2M-1). IN GENERAL, THE BETA
C                                    VALUES FOR THE OPTIMAL SUBSET OF
C                                    SIZE K ARE STORED IN POSITIONS
C                                    BVAL(L),...,BVAL(L-K+1) WHERE
C                                       L=(M*(M+1)-K*(K+1))/2 + 1
C
C  IDEX     INTEGER ARRAY    OUTPUT: BETA INDEX SET FOR THE OPTIMAL
C        (((MMAX+1)*MMAX)/2)         SUBSET. THIS ARRAY IS A  PARALLEL
C                                    ARRAY FOR BVAL; I.E., IF BVAL(J)=2.
C                                    AND IDEX(J)=7  THEN BETA(7)=2.7 IN
C                                    THE ASSOCIATED OPTIMAL SUBSET.
C
C  ISTAT    INTEGER ARRAY    INPUT:  PARAMETER STATUS ARRAY.
C            (MMAX)
C
C                                        1 IF BETA(J) IS REQUIRED
C                                          IN EVERY MODEL
C                     ISTAT(J)=
C
C                                        0 OTHERWISE
C
C  ZL       REAL ARRAY       OUTPUT: BEST OBJECTIVE VALUE FOR EACH SUBSE
C            (MMAX)
C                                    ZL(J) GIVES THE BEST OBJECTIVE VALU
C                                    FOR THE SUBSET WITH M-J+1 PARAMETER
C***********************************************************************
C
C    IMPLEMENTATION NOTES:
C
C        1. THE ROUTINE USES TWO MACHINE DEPENDENT VALUES ACU AND BIG.
C           (THESE ARE SET AT THE BEGINNING OF THIS SUBROUTINE)
C           A) ACU IS USED TO TEST FOR "ZERO". ACU SHOULD BE SET TO
C              APPROXIMATELY 100 * THE RELATIVE MACHINE ACCURACY OF THE
C              SYSTEM IN USE.
C           B) BIG IS USED TO INITIALIZE THE ARRAY ZL (DESCRIBED ABOVE).
C              BIG SHOULD BE SET TO THE LARGEST FLOATING POINT VALUE
C              ASSIGNABLE.
C
C        2. BOTH SINGLE AND DOUBLE PRECISION VERSIONS ARE SUPPLIED. THIS
C           VERSION IS IN DOUBLE PRECISION. TO GET A SINGLE PRECISION CO
C           ALL STATEMENTS HAVING A "C++" IN COLUMNS 1-3 SHOULD BE MODIF
C           BY DELETING THE "C++". THE CORRESPONDING DOUBLE PRECISION
C           STATEMENTS SHOULD BE EITHER DELETED OR COMMENTED OUT.
C
C        3. ARRAY DIMENSIONS
C           THE CODE IS CURRENTLY DIMENSIONED TO SOLVE PROBLEMS WITH UP
C           20 PARAMETERS AND 300 OBSERVATIONS. THE DIMENSION SIZES ARE
C           DETERMINED AS FOLLOWS
C               20 MMAX
C              300 NMAX
C              210 (NMAX+1)*MMAX/2
C             6000 NMAX*MMAX
C               10 MAXIMUM VALUE OF IK .
C
C            IF THE DIMENSION SIZES ARE CHANGED THE COMMON AND DIMENSION
C            STATEMENTS WILL NEED TO BE MODIFIED IN EACH SUBROUTINE.
C
C
C***********************************************************************
C
C     SUBROUTINES:
C
C         CALBET   - BACK-SOLVES A SYSYTEM OF EQUATIONS
C
C         CALCPI   - FORWARD-SOLVES A SYSYTEM OF EQUATIONS
C
C         KBEST    - THE DRIVER
C
C         L1NORM   - SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
C                    PARAMETERS INCLUDED IN THE MODEL
C
C         PHASE2   - SOLVES THE CURRENT REGRESSION PROBLEM WITH A PRIMAL
C                    ALGORITHM
C
C         SETUP    - DETERMINES THE FORM OF THE CURRENT PROBLEM BY
C                    CHOOSING THE PARAMETER TO LEAVE BASED ON A PENALTY
C
C         UPDATE   - UPDATES LU DECOMPOSITION MATRIX
C
C
C***********************************************************************
C
C
C     DESCRIPTION OF VARIABLES:
C
C     IK:LENGTH OF CANDIDATE LIST
C     BETA:ESTIMATES TO THE CURRENT SUBPROBLEM
C     SAD:THE MINIMUM TOTAL ABSOLUTE DEVIATION
C     IBASE:THE INDEX ARRAY OF COLUMNS OF THE BASIS
C
C        N,M,NMAX,MMAX,X,Y, ARE NOT CHANGED IN THE SUBROUTINE;
C          SAD IS UPDATED AT EACH ITERATION
C
C     LU:THE LU DECOMPOSITION OF THE CURRENT BASIS
C     INDEX:THE INDEX ARRAY OF ROWS OF LU
C     TOT:THE CURRENT RHS OF THE DUAL PROBLEM
C     SIGMA:INDICATOR ARRAY TO SPECIFY WHETHER A NONBASIC DUAL VARIABLE
C           IS AT UPPER OR LOWER BOUND(+1 IMPLIES UPPER;-1 IMPLIES
C           LOWER)
C     INEXT:A LOCAL ARRAY FOR THE SORT ROUTINE
C     RHS: A LOCAL ARRAY FOR THE CALBET ROUTINE
C
C     IPAR(I)=-K IF THE K-TH PARAMETER (BETA) IS OUT OF MODEL AT LEVEL I
C         = K IF THE K-TH PARAMETER IS IN
C
C     LABEL(I)  SAVES THE INDICES OF THE BASIC OBSERVATIONS IN THE
C            PREDECESSOR PATH
C     ZSAVE(I)  SAVES THE OBJECTIVE VALUES ON THE PREDECESSOR PATH
C
C**********************************************************************
C
      DIMENSION X(300,20), Y(300), ZL(20), BVAL(210), IDEX(210), ISTAT(2
     10)
C
      DIMENSION BETA(20), TOT(20), PI(20), REPP(20)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,IDEX,IFAULT,IK,IK1,INCOL,INDEX,INPROB,IPAR
      INTEGER ISTAT,ITER,J,JINM,K2,K3,KKK1,L,LABEL,LEVEL,M,MININ
      INTEGER MMAX,N,NFREE,NMAX,NUMIN
C
      DOUBLE PRECISION ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
C++   REAL             ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
      DOUBLE PRECISION REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
C++   REAL             REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
      DOUBLE PRECISION ZERO,ZL,ZLOW,ZSAVE
C++   REAL             ZERO,ZL,ZLOW,ZSAVE
C
      LOGICAL DIRECT,INTL
C
      DATA ZERO/0.0D0/
C++   DATA ZERO/0.0E0/
C
C     ASSIGN VALUES TO MACHINE DEPENDENT CONSTANTS
C    SEE - FOX, P.A., A.D. HALL AND N.L. SCHRYER, "FRAMEWORK FOR
C          A PORTABLE LIBRARY: ALGORITHM 528," ACM TRANSACTIONS
C          ON MATHEMATICAL SOFTWARE, VOL 4, NO 2, JUNE 1978.
C
      ACU = SQRT(D1MACH(4))
       BIG = D1MACH(2)
C++    ACU = R1MACH(1)
C++    ACU = 100.*ACU
C++    BIG = R1MACH(2)
C
C
      IFAULT=0
      IK=5
      IK1=IK-1
      POPT1=(100.-POPT)/100.
C
C     TEST PROBLEM SIZE
C
      IF (N.LE.NMAX.AND.M.LE.MMAX) GO TO 10
      IFAULT=2
      RETURN
   10 CONTINUE
      DO 20 I=1,M
         INDEX(I)=I
         TOT(I)=ZERO
   20 INCOL(I)=I
      CALL L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
      IF (IFAULT.GT.0) RETURN
C
C     INITIALIZATION
C
C     SAVE THE INITIAL SOLUTION
C
      JINM=0
      DO 30 I=1,M
         BSAVE(I)=BETA(I)
         BVAL(I)=BETA(I)
         PIE(I)=PI(I)
         LABEL(I)=IBASE(I)
         IDEX(I)=I
         INPROB(I)=I
         SAVTOT(I)=TOT(I)
         ZL(I)=BIG
         JINM=JINM+ISTAT(I)
   30 CONTINUE
      DO 40 I=1,N
         SIG(I)=SIGMA(I)
   40 CONTINUE
C
      LEVEL=0
      DIRECT=.TRUE.
C
      IF (JINM.GT.MININ) MININ=JINM
C
C     INITIALIZE NUMBER OF FREE PARAMETERS IN MODEL
C
      NFREE=M-JINM
C
      ZL(M)=SAD
      ZSAVE(M)=SAD
C
      IF (MININ.EQ.M) RETURN
C
C       NUMIN GIVES THE NUMBER OF PARAMETERS IN THE MODEL
C
      NUMIN=M
C
C     START THE MAIN LOOP
C
   50 NUMIN=NUMIN-1
      LEVEL=LEVEL+1
      IF (NUMIN.LT.MININ.OR.NFREE.EQ.0) GO TO 120
      K2=(M*(M+1)-(NUMIN+1)*(NUMIN+2))/2
      KKK1=NUMIN+1
      ZLOW=ZL(NUMIN)*POPT1
      IF (DIRECT) GO TO 80
      DO 60 I=1,KKK1
         J=K2+I
         INDEX(I)=I
         IBASE(I)=LABEL(J)
         INCOL(I)=INPROB(J)
         PI(I)=PIE(J)
         BETA(I)=BSAVE(J)
         TOT(I)=SAVTOT(J)
   60 CONTINUE
C
C     LOAD IN THE BOUNDS FROM THE IMMEDIATE PREDECESSOR
C
      K3=(M-KKK1)*N
      DO 70 I=1,N
         K3=K3+1
         SIGMA(I)=SIG(K3)
   70 CONTINUE
      SAD=ZSAVE(KKK1)
C
C     SET UP A NEW PROBLEM
C
   80 CALL SETUP (X,KKK1,IFAULT,ISTAT,BETA,TOT,PI,REPP)
      IF (IFAULT.GT.0) RETURN
      NFREE=NFREE-1
C
C     CALL PRIMAL L.P. CODE
C
      CALL PHASE2 (X,Y,NUMIN,N,ITER,IFAULT,BETA,TOT,PI,REPP)
C
C      SAVE THE SOLUTION DATA FOR LATER RECALL
C
      K2=K2+KKK1+1
      J=K2
      DO 90 I=1,NUMIN
         L=K2+I
         PIE(L)=PI(I)
         LABEL(L)=IBASE(I)
         BSAVE(L)=BETA(I)
         INPROB(L)=INCOL(I)
         SAVTOT(L)=TOT(I)
   90 CONTINUE
C
C     SAVE THE OBJECTIVE VALUE
C
      ZSAVE(NUMIN)=SAD
C
C     SAVE THE NONBASIC BOUNDS
C
      K3=(M-NUMIN)*N
      DO 100 I=1,N
         K3=K3+1
         SIG(K3)=SIGMA(I)
  100 CONTINUE
C
      DIRECT=.TRUE.
C
C     CHECK FOR A BETTER SOLUTION
C
      IF (SADT.GE.ZLOW) GO TO 50
C
C     SAVE THE BETTER SOLUTION
C
      ZL(NUMIN)=SAD
      K2=J
      DO 110 I=1,NUMIN
         L=K2+I
         IDEX(L)=INCOL(I)
         BVAL(L)=BETA(I)
  110 CONTINUE
C
C     GO TO TOP OF LOOP
C
      GO TO 50
  120 NUMIN=NUMIN+1
C
      DIRECT=.FALSE.
C
C     MUST WORK BACK UP THE TREE
C
      LEVEL=LEVEL-1
  130 IF (IPAR(LEVEL).GT.0) GO TO 140
      IPAR(LEVEL)=-IPAR(LEVEL)
      J=IPAR(LEVEL)
      ISTAT(J)=1
      NUMIN=NUMIN+1
      GO TO 50
  140 J=IPAR(LEVEL)
      ISTAT(J)=0
      NFREE=NFREE+1
      LEVEL=LEVEL-1
      IF (LEVEL.GT.0) GO TO 130
      RETURN
C
      END
      SUBROUTINE SETUP (X,M,IFAULT,ISTAT,BETA,TOT,PI,REPP)
C
C     THIS SUBROUTINE DETERMINES THE FORM OF THE CURRENT SUBPROBLEM
C
      DIMENSION X(300,20), ISTAT(20), RHS(20), BETA(20), TOT(20), PI(20)
     1, REPP(20)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,IFAULT,IFIXI,II,IK,IK1,INCOL,INDEX,INPROB
      INTEGER IOUT,IPAR,ISAVE,ISTAT,KKK,L,LABEL,LEAVE,LEVEL,LL,M
C
      DOUBLE PRECISION ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
C++   REAL             ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
      DOUBLE PRECISION RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
C++   REAL             RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
      DOUBLE PRECISION SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
C++   REAL             SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
C
      LOGICAL DIRECT,INTL
C
      DATA ONE/1.0D0/,ZERO/0.0D0/
C++   DATA ONE/1.0E0/,ZERO/0.0E0/
C
C     IF THE IMMEDIATE PREDECESSOR IS IN MEMORY GO ON
C
      IF (DIRECT) GO TO 10
C
C     RECONSTRUCT THE LU DECOMPOSITION
C
      KKK=1
      CALL UPDATE (KKK,IFAULT,X,1,M)
      IF (IFAULT.GT.0) RETURN
C
C     RECALCULATE THE BASIC PI VALUES
C
      KKK=0
C
      CALL CALCPI (KKK,PI,TOT,M)
C
C     DETERMINE PARAMETER TO LEAVE BASED ON PENALTY
C
   10 PENALT=BIG
      KKK=0
      DO 20 L=1,M
         RHS(L)=ZERO
   20 CONTINUE
C
      DO 40 I=1,M
         II=INDEX(I)
         LL=INCOL(II)
         IF (ISTAT(LL).EQ.1) GO TO 40
C++      RHO=SIGN(1.,BETA(II))
         RHO=DSIGN(1.D0,BETA(II))
         RHS(II)=ONE
         KKK=I-1
         CALL CALCPI (KKK,REPP,RHS,M)
         RHS(II)=ZERO
         TEST=BIG
         DO 30 L=1,M
C++         IF ( ABS(REPP(L)).LE.ACU) GO TO 30
            IF (DABS(REPP(L)).LE.ACU) GO TO 30
C++         SREPP=SIGN(1.,REPP(L))
            SREPP=DSIGN(1.D0,REPP(L))
C++         RATIO=ABS((1.-RHO*SREPP*PI(L))/REPP(L))
            RATIO=DABS((1.-RHO*SREPP*PI(L))/REPP(L))
            IF (RATIO.GE.TEST) GO TO 30
            SIDE=RHO*SREPP
            TEST=RATIO
            IOUT=L
   30    CONTINUE
C
C     SEE IF THE PENALTY IS LESS
C
C++      IF (TEST*ABS(BETA(II)).GE.PENALT) GO TO 40
         IF (TEST*DABS(BETA(II)).GE.PENALT) GO TO 40
         LEAVE=IOUT
         IFIXI=II
         BESIDE=SIDE
C++      PENALT=TEST*ABS(BETA(II))
         PENALT=TEST*DABS(BETA(II))
   40 CONTINUE
C
C     UPDATE THE OBJECTIVE VALUE
C
      SAD=SAD+PENALT
C
C     PARAMETER INCOL(IFIXI) WILL LEAVE THE MODEL
C     PI(IBASE(LEAVE)) WILL LEAVE THE BASIS AND GO TO BOUND
C
C     SWITCH UNWANTED PARAMETER AND PI TO THE END OF LIST
C
      ISAVE=INCOL(IFIXI)
      INCOL(IFIXI)=INCOL(M)
      INCOL(M)=ISAVE
      TOT(IFIXI)=TOT(M)
C
C     LABEL THIS NODE AND FLAG THE OUTGOING PARAMETER
C
      IPAR(LEVEL)=-ISAVE
      ISTAT(ISAVE)=-1
C
      ISAVE=IBASE(LEAVE)
      IBASE(LEAVE)=IBASE(M)
      IBASE(M)=ISAVE
C
C     PLACE THE OUTGOING PI AT THE CORRECT BOUND
C
      SIGMA(ISAVE)=BESIDE
C
C     FORM A NEW BASIS
C
      M=M-1
C
      DO 50 I=1,M
   50 INDEX(I)=I
C
      KKK=1
      CALL UPDATE (KKK,IFAULT,X,1,M)
      IF (IFAULT.GT.0) RETURN
C
C     UPDATE THE TOT ARRAY
C
      DO 60 I=1,M
         LL=INCOL(I)
         TOT(I)=TOT(I)-BESIDE*X(ISAVE,LL)
   60 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE PHASE2 (X,Y,M,N,ITER,IFAULT,BETA,TOT,PI,REPP)
C
C     SOLVES LP WITH PRIMAL ALGORITHM
C
      DIMENSION X(300,20), Y(300), BETA(20), TOT(20), PI(20), REPP(20),
     1RHS(20),Z(10),KAN(10)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INDEX,INPROB,IOUT,IPAR
      INTEGER ITER,J,K,KAN,KKK,KROW,L,LABEL,LEVEL,LL,M,N
C
      DOUBLE PRECISION ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
C++   REAL             ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
      DOUBLE PRECISION RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
C++   REAL             RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
      DOUBLE PRECISION TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
C++   REAL             TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
C
      LOGICAL DIRECT,INTL
C
      DATA ZERO/0.0D0/,ONE/1.0D0/
C++   DATA ZERO/0.0E0/,ONE/1.0E0/
C
C     Z(I) STORES THE REDUCED COST VALUES ON THE LIST
C     KAN(I) STORES THE ROW NUMBER ON THE LIST
C     IK IS THE LENGTH OF THE LIST
C
C     CALCULATE THE SIMPLEX MULTIPLIERS
C
      DO 10 I=1,M
         K=IBASE(I)
   10 PI(I)=Y(K)
      KKK=0
      CALL CALBET (KKK,BETA,PI,M)
C
C     RECALCULATE THE BASIC PI VALUES
C
      KKK=0
      CALL CALCPI (KKK,PI,TOT,M)
C
C     STORE THE OBJECTIVE VALUE USED FOR TERMINATION CHECK
C
      SADT=SAD
      IF (SAD.GE.ZLOW) RETURN
C
C     GENERATE THE CANDIDATE LIST
C
   20 DO 30 I=1,IK
         Z(I)=ACU
   30 CONTINUE
C
C     CALCULATE THE REDUCED COSTS
C
      DO 70 J=1,N
         IF (SIGMA(J).EQ.ZERO) GO TO 70
         ZC=ZERO
         DO 40 I=1,M
            LL=INCOL(I)
   40    ZC=ZC+BETA(I)*X(J,LL)
         ZZ=(ZC-Y(J))*SIGMA(J)
         IF (ZZ.LE.Z(IK)) GO TO 70
C
C     RANK THE VALUES IN Z(I), IN DESCENDING ORDER
C
         DO 60 K=1,IK1
            IF (ZZ.LE.Z(K)) GO TO 60
            I=IK
   50       Z(I)=Z(I-1)
            KAN(I)=KAN(I-1)
            I=I-1
            IF (I.GT.K) GO TO 50
            Z(K)=ZZ
            KAN(K)=J
            GO TO 70
   60    CONTINUE
         Z(IK)=ZZ
         KAN(IK)=J
   70 CONTINUE
      IF (Z(1).LE.ACU) RETURN
      KROW=1
      ZZ=Z(1)
      IIN=KAN(1)
      ZC=ZZ*SIGMA(IIN)
      GO TO 100
   80 KROW=KROW+1
      IF (KROW.GT.IK) GO TO 20
      IF (Z(KROW).LE.ACU) GO TO 20
C
C     DETERMINE THE ENTERING VARIABLE FROM THE CANDIDATE LIST
C
      IIN=KAN(KROW)
      ZC=ZERO
      DO 90 I=1,M
         LL=INCOL(I)
   90 ZC=ZC+BETA(I)*X(IIN,LL)
      ZC=ZC-Y(IIN)
      ZZ=ZC*SIGMA(IIN)
      IF (ZZ.LE.ACU) GO TO 80
  100 BEST=ZZ
      RHO=SIGMA(IIN)
C
C     RHO: SIGN OF THE INCOMING REDUCED COST
C     IIN: INCOMING CANDIDATE
C
C     FIND THE REPRESENTATION OF THE ENTERING VARIABLE
C
      DO 110 I=1,M
         LL=INCOL(I)
         RHS(I)=X(IIN,LL)
  110 CONTINUE
      KKK=0
      CALL CALCPI (KKK,REPP,RHS,M)
C
C     CALCULATE THE MIN RATIO TEST TO DETERMINE THE LEAVING VARIABLE
C
      TEST=2.0
      DO 120 I=1,M
C++      IF (ABS(REPP(I)).LE.ACU) GO TO 120
         IF (DABS(REPP(I)).LE.ACU) GO TO 120
C++      SREPP=SIGN(1.0,REPP(I))
         SREPP=DSIGN(1.0D0,REPP(I))
C++      RATIO=ABS((1.-RHO*SREPP*PI(I))/REPP(I))
         RATIO=DABS((1.-RHO*SREPP*PI(I))/REPP(I))
         IF (RATIO.GE.TEST) GO TO 120
         TEST=RATIO
         IOUT=I
  120 CONTINUE
C
C     PERFORM FATHOMING TEST BEFORE PIVOT
C
      SADT=SAD+TEST*ZZ
      IF (SADT.GE.ZLOW) RETURN
      SAD=SADT
C
      DO 130 I=1,M
  130 PI(I)=PI(I)+TEST*REPP(I)*RHO
      IF (TEST.LT.2.0) GO TO 150
C
      SIGMA(IIN)=-SIGMA(IIN)
      XTWO=SIGMA(IIN)+SIGMA(IIN)
      DO 140 I=1,M
         LL=INCOL(I)
         TOT(I)=TOT(I)-XTWO*X(IIN,LL)
  140 CONTINUE
      GO TO 80
  150 K=IBASE(IOUT)
C++   SIGMA(K)=SIGN(1.0,PI(IOUT))
      SIGMA(K)=DSIGN(1.0D0,PI(IOUT))
      PI(IOUT)=SIGMA(IIN)*(1.0-TEST)
      IBASE(IOUT)=IIN
C
      DO 160 I=1,M
         LL=INCOL(I)
         TOT(I)=TOT(I)+SIGMA(IIN)*X(IIN,LL)-SIGMA(K)*X(K,LL)
  160 CONTINUE
C
      KKK=IOUT
      CALL UPDATE (IOUT,IFAULT,X,1,M)
C
      DO 170 I=1,M
         RHS(I)=ZERO
  170 CONTINUE
      RHS(KKK)=ONE
      KKK=KKK-1
      CALL CALBET (KKK,REPP,RHS,M)
      DO 180 L=1,M
  180 BETA(L)=BETA(L)-ZC*REPP(L)
      SIGMA(IIN)=ZERO
      ITER=ITER+1
      GO TO 80
C
      END
      SUBROUTINE L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
C
C     THIS ROUTINE  SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
C     THE PARAMETERS INCLUDED IN THE MODEL
C
      DIMENSION X(300,20), Y(300), D(300), DELTA(300), INEXT(300), PRICE
     1(300), RHS(20), BETA(20), TOT(20), PI(20)
C
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INEXT,INPROB,IOUT,IPAR
      INTEGER INDEX,IPT,ITER,J,K,K1,KKK,KOUNT,L,LABEL,LEVEL,M,M1,N
C
      DOUBLE PRECISION ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
C++   REAL             ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
      DOUBLE PRECISION ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
C++   REAL             ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
      DOUBLE PRECISION SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
C++   REAL             SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
      DOUBLE PRECISION ZSAVE,ZZZ
C++   REAL             ZSAVE,ZZZ
C
      LOGICAL DIRECT,INTL
C
      DATA ONE/1.0D0/,ZERO/0.0D0/
C++   DATA ONE/1.0E0/,ZERO/0.0E0/
C
C        INITIAL SETTINGS
C
      ITER=0
      AONE=ONE+ACU
      AHALF=.5+ACU
      BONE=ONE-ACU
      M1=M-1
C
C     SET UP INITIAL LU DECOMPOSITION
C
      INTL=.TRUE.
      K=1
      CALL UPDATE (K,IFAULT,X,N,M)
      IF (IFAULT.NE.0) RETURN
      INTL=.FALSE.
      DO 10 I=1,M
         K1=IBASE(I)
         RHS(I)=Y(K1)
   10 CONTINUE
      KKK=0
      CALL CALBET (KKK,BETA,RHS,M)
C
C        CALCULATE INITIAL D, TOT AND SIGMA VECTORS
C
      DO 30 J=1,N
         VAL=ZERO
         DO 20 I=1,M
            VAL=VAL+BETA(I)*X(J,I)
   20    CONTINUE
         D(J)=Y(J)-VAL
C++      SIGMA(J)=SIGN(1.,D(J))
         SIGMA(J)=DSIGN(1.D0,D(J))
   30 CONTINUE
      DO 40 J=1,M
         RHS(J)=ZERO
         KKK=IBASE(J)
         SIGMA(KKK)=ZERO
   40 CONTINUE
      DO 60 J=1,N
         DO 50 I=1,M
   50    TOT(I)=TOT(I)-SIGMA(J)*X(J,I)
   60 CONTINUE
C
C        MAIN ITERATIVE LOOP BEGINS
C
   70 CONTINUE
C
      KKK=0
      CALL CALCPI (KKK,PI,TOT,M)
C
      T=AONE
      K=0
      DO 80 J=1,M
C++      IF (ABS(PI(J)).LT.T) GO TO 80
         IF (DABS(PI(J)).LT.T) GO TO 80
         K=J
C++      T=ABS(PI(J))
         T=DABS(PI(J))
C++      RHO=-SIGN(1.,PI(J))
         RHO=-DSIGN(1.D0,PI(J))
   80 CONTINUE
      IF (K.EQ.0) GO TO 220
      KKK=K-1
      RHS(K)=ONE
      CALL CALBET (KKK,BETA,RHS,M)
      RHS(K)=ZERO
      DO 100 I=1,N
         DELTA(I)=ZERO
         IF (SIGMA(I).EQ.ZERO) GO TO 100
         DO 90 J=1,M
            DELTA(I)=DELTA(I)+BETA(J)*X(I,J)
   90    CONTINUE
         DELTA(I)=RHO*DELTA(I)
  100 CONTINUE
C
C        PERFORM PARTIAL SORT OF RATIOS
C
      T=T*.5
      KOUNT=0
      RATIO=BIG
      SUM=AHALF
      SUBT=ZERO
      DO 160 I=1,N
         IF (DELTA(I)*SIGMA(I).LE.ACU) GO TO 160
         TEST=D(I)/DELTA(I)
         IF (TEST.GE.RATIO) GO TO 160
C++      SUM=SUM+ABS(DELTA(I))
         SUM=SUM+DABS(DELTA(I))
         IF (SUM-SUBT.GE.T) GO TO 110
C
C     INSERT I IN LIST
C
         KOUNT=KOUNT+1
         PRICE(KOUNT)=TEST
         INEXT(KOUNT)=I
         GO TO 160
C
C     UPDATE SUM AND KICK IIN OUT OF THE LIST
C
  110    SUM=SUM-SUBT
         RATIO=TEST
         IPT=0
         KKK=0
C
C     IDENTIFY A NEW IIN
C
  120    KKK=KKK+1
         IF (KKK.GT.KOUNT) GO TO 130
         IF (PRICE(KKK).LE.RATIO) GO TO 120
         RATIO=PRICE(KKK)
         IPT=KKK
         GO TO 120
  130    IF (IPT.EQ.0) GO TO 150
C
C     SWITCH VALUES
C
         KKK=INEXT(IPT)
C++      SUBT=ABS(DELTA(KKK))
         SUBT=DABS(DELTA(KKK))
         IF (SUM-SUBT.LT.T) GO TO 140
         PRICE(IPT)=PRICE(KOUNT)
         INEXT(IPT)=INEXT(KOUNT)
         KOUNT=KOUNT-1
         GO TO 110
  140    IIN=INEXT(IPT)
         INEXT(IPT)=I
         PRICE(IPT)=TEST
         GO TO 160
  150    IIN=I
C++      SUBT=ABS(DELTA(I))
         SUBT=DABS(DELTA(I))
  160 CONTINUE
C
C        UPDATE BASIC INDICATORS
C
      IF (KOUNT.EQ.0) GO TO 190
      DO 180 J=1,KOUNT
         KKK=INEXT(J)
         ZZZ=SIGMA(KKK)
         DO 170 L=1,M
            SUBT=ZZZ*X(KKK,L)
            TOT(L)=TOT(L)+SUBT+SUBT
  170    CONTINUE
         SIGMA(KKK)=-SIGMA(KKK)
  180 CONTINUE
  190 CONTINUE
      IOUT=IBASE(K)
      DELTA(IOUT)=RHO
      DO 200 L=1,M
  200 TOT(L)=TOT(L)+RHO*X(IOUT,L)+SIGMA(IIN)*X(IIN,L)
      SIGMA(IOUT)=-RHO
      IBASE(K)=IIN
      CALL UPDATE (K,IFAULT,X,1,M)
      DO 210 J=1,N
  210 D(J)=D(J)-RATIO*DELTA(J)
      SIGMA(IIN)=ZERO
      ITER=ITER+1
      GO TO 70
C
C        CALCULATE OPTIMAL BETA AND SUM OF ABSOLUTE DEVIATIONS
C
  220 DO 230 I=1,M
         K1=IBASE(I)
         RHS(I)=Y(K1)
  230 CONTINUE
      KKK=0
      CALL CALBET (KKK,BETA,RHS,M)
      SAD=ZERO
      DO 240 J=1,N
C++  240 SAD=SAD+ABS(D(J))
 240  SAD=SAD+DABS(D(J))
      RETURN
C
      END
      SUBROUTINE CALBET (KKK,BETA,RHS,M)
C
C     SUBROUTINE CALBET FOR BACK-SOLVING SYSTEM OF EQUATIONS
C
      DIMENSION BETA(20), RHS(20)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2
      INTEGER KK,KKK,KKK1,LABEL,LEVEL,M,M1
C
      DOUBLE PRECISION ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
C++   REAL             ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
      DOUBLE PRECISION SIG,SIGMA,ZERO,ZLOW,ZSAVE
C++   REAL             SIG,SIGMA,ZERO,ZLOW,ZSAVE
C
      LOGICAL DIRECT,INTL
C
      DATA ZERO/0.0D0/
C++   DATA ZERO/0.0E0/
C
      IF (M.GT.1) GO TO 10
      K=INDEX(1)
      BETA(K)=RHS(1)/LU(K,1)
      RETURN
   10 CONTINUE
      M1=M-1
      IF (KKK.EQ.0) GO TO 30
      DO 20 I=1,KKK
         K=INDEX(I)
         BETA(K)=ZERO
   20 CONTINUE
   30 KKK=KKK+1
      K=INDEX(KKK)
      BETA(K)=RHS(KKK)/LU(K,KKK)
      IF (KKK.EQ.M) GO TO 60
      KKK1=KKK+1
      DO 50 II=KKK1,M
         K=INDEX(II)
         BETA(K)=RHS(II)
         II1=II-1
         DO 40 I=KKK,II1
            KK=INDEX(I)
            BETA(K)=BETA(K)-LU(KK,II)*BETA(KK)
   40    CONTINUE
         BETA(K)=BETA(K)/LU(K,II)
   50 CONTINUE
   60 CONTINUE
      DO 70 II=1,M1
         K1=M-II
         K=INDEX(K1)
      DO 70 I=1,II
         KK=M-I+1
         K2=INDEX(KK)
         BETA(K)=BETA(K)-LU(K2,K1)*BETA(K2)
   70 CONTINUE
      RETURN
C
      END
      SUBROUTINE UPDATE (KKK,IFAULT,X,N,M)
C
C     SUBROUTINE UPDATE FOR UPDATING LU DECOMPOSITION MATRIX
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,ICOL,IFAULT,II,II1,IK,IK1,INCOL,INDEX,INPROB
      INTEGER IPAR,IROW,ISAVE,J,K,KK,KKK,LABEL,LEVEL,LL,M,N
C
      DOUBLE PRECISION ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
C++   REAL             ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
      DOUBLE PRECISION SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
C++   REAL             SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
C
      LOGICAL DIRECT,INTL
C
      DIMENSION X(300,20)
C
      IROW=0
      DO 90 II=KKK,M
         IF (INTL) GO TO 10
         IROW=IBASE(II)
         GO TO 20
   10    IROW=IROW+1
         IBASE(II)=IROW
         IF (IROW.GT.N) IFAULT=1
         IF (IFAULT.NE.0) RETURN
   20    DO 30 I=1,M
            LL=INCOL(I)
            LU(I,II)=X(IROW,LL)
   30    CONTINUE
C
C     SET UP REPRESENTATION OF INCOMING ROW
C
         IF (II.EQ.1) GO TO 60
         II1=II-1
         DO 50 ICOL=1,II1
            K=INDEX(ICOL)
            SUBT=LU(K,II)
            J=ICOL+1
            DO 40 I=J,M
               K=INDEX(I)
               LU(K,II)=LU(K,II)-SUBT*LU(K,ICOL)
   40       CONTINUE
   50    CONTINUE
C
C     FIND MAXIMUM ENTRY
C
   60    PIVOT=ACU
         KK=0
         DO 70 I=II,M
            K=INDEX(I)
            IF (DABS(LU(K,II)).LE.PIVOT) GO TO 70
C++         IF (ABS(LU(K,II)).LE.PIVOT) GO TO 70
C++         PIVOT=ABS(LU(K,II))
            PIVOT=DABS(LU(K,II))
            KK=I
   70    CONTINUE
         IF (KK.EQ.0) GO TO 10
C
C     SWITCH ORDER
C
         ISAVE=INDEX(KK)
         INDEX(KK)=INDEX(II)
         INDEX(II)=ISAVE
C
C     PUT IN COLUMNS OF LU ONE AT A TIME
C
         IF (II.EQ.M) GO TO 90
         J=II+1
         DO 80 I=J,M
            K=INDEX(I)
            LU(K,II)=LU(K,II)/LU(ISAVE,II)
   80    CONTINUE
   90 CONTINUE
      KKK=IROW
      RETURN
C
      END
      SUBROUTINE CALCPI (KKK,PI,TOT,M)
C
C     THIS ROUTINE FORWARD SOLVES A LINEAR SYSTEM
C
      DIMENSION TOT(20), PI(20)
C
      COMMON /BLOCK/
     -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
     -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
     -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
C
      INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2,KKK
      INTEGER KKK1,LABEL,LEVEL,M,M1
C
      DOUBLE PRECISION ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
C++   REAL             ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
      DOUBLE PRECISION SIGMA,TOT,ZERO,ZLOW,ZSAVE
C++   REAL             SIGMA,TOT,ZERO,ZLOW,ZSAVE
C
      LOGICAL DIRECT,INTL
C
C++   DATA ZERO/0.0E0/
      DATA ZERO/0.0D0/
      M1=M-1
C
      IF (M.GT.1) GO TO 10
      K=INDEX(M)
      PI(M)=TOT(K)/LU(K,M)
      RETURN
   10 IF (KKK.EQ.0) GO TO 30
      DO 20 I=1,KKK
         PI(I)=ZERO
   20 CONTINUE
   30 CONTINUE
      KKK=KKK+1
      K=INDEX(KKK)
      PI(KKK)=TOT(K)
      IF (KKK.EQ.M) GO TO 60
      KKK1=KKK+1
      DO 50 II=KKK1,M
         K=INDEX(II)
         PI(II)=TOT(K)
         II1=II-1
         DO 40 I=KKK,II1
   40    PI(II)=PI(II)-LU(K,I)*PI(I)
   50 CONTINUE
   60 CONTINUE
      K=INDEX(M)
      PI(M)=PI(M)/LU(K,M)
      DO 80 II=1,M1
         K1=M-II
         K=INDEX(K1)
         DO 70 I=1,II
            K2=M-I+1
            PI(K1)=PI(K1)-LU(K,K2)*PI(K2)
   70    CONTINUE
         PI(K1)=PI(K1)/LU(K,K1)
   80 CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION D1MACH(I)
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C     AS FOLLOWS, FOR EXAMPLE
C
C          D = D1MACH(I)
C
C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  DOUBLE-PRECISION MACHINE CONSTANTS
C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C  D1MACH( 5) = LOG10(B)
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA SMALL(1) / ZC00800000 /
C     DATA SMALL(2) / Z000000000 /
C
C     DATA LARGE(1) / ZDFFFFFFFF /
C     DATA LARGE(2) / ZFFFFFFFFF /
C
C     DATA RIGHT(1) / ZCC5800000 /
C     DATA RIGHT(2) / Z000000000 /
C
C     DATA DIVER(1) / ZCC6800000 /
C     DATA DIVER(2) / Z000000000 /
C
C     DATA LOG10(1) / ZD00E730E7 /
C     DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O0000000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O0007777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O7770000000000000 /
C
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O7777777777777777 /
C
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA SMALL(1) / 00604000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 37767777777777777777B /
C     DATA LARGE(2) / 37167777777777777777B /
C
C     DATA RIGHT(1) / 15604000000000000000B /
C     DATA RIGHT(2) / 15000000000000000000B /
C
C     DATA DIVER(1) / 15614000000000000000B /
C     DATA DIVER(2) / 15010000000000000000B /
C
C     DATA LOG10(1) / 17164642023241175717B /
C     DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA SMALL(1) / 200004000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C
C     DATA LARGE(1) / 577777777777777777777B /
C     DATA LARGE(2) / 000007777777777777777B /
C
C     DATA RIGHT(1) / 377214000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C
C     DATA DIVER(1) / 377224000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC DMACH(5)
C
C     DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C     DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C     DATA LOG10/40423K,42023K,50237K,74776K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
C     DATA LARGE(1),LARGE(2) / [37777777, [37777577 /
C     DATA RIGHT(1),RIGHT(2) / [20000000, [00000333 /
C     DATA DIVER(1),DIVER(2) / [20000000, [00000334 /
C     DATA LOG10(1),LOG10(2) / [23210115, [10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
C     DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
C     DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
C     DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
C     DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /
C
C      MACHINE CONSTANTS FOR THE HP 2100
C      THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C      DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
C      DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C      DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
C      DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
C      DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
C
C
C      MACHINE CONSTANTS FOR THE HP 2100
C      FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C      DATA SMALL(1), SMALL(2) /  40000B,       0 /
C      DATA SMALL(3), SMALL(4) /       0,       1 /
C      DATA LARGE(1), LARGE(2) /  77777B, 177777B /
C      DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C      DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
C      DATA RIGHT(3), RIGHT(4) /       0,    225B /
C      DATA DIVER(1), DIVER(2) /  40000B,       0 /
C      DATA DIVER(3), DIVER(4) /       0,    227B /
C      DATA LOG10(1), LOG10(2) /  46420B,  46502B /
C      DATA LOG10(3), LOG10(4) /  76747B, 176377B /
C
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
C     DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /
C     DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C     DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /
C     DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /
C     DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C     DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /
C     DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /
C     DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /
C     DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /
C     DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
C
C     DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /
C     DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /
C     DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /
C     DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /
C     DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /    8388608,           0 /
C     DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
C     DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
C     DATA DIVER(1),DIVER(2) /  620756992,           0 /
C     DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /
C
C     DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /
C     DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /
C     DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /
C     DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /
C     DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /    128,      0 /
C     DATA SMALL(3),SMALL(4) /      0,      0 /
C
C     DATA LARGE(1),LARGE(2) /  32767,     -1 /
C     DATA LARGE(3),LARGE(4) /     -1,     -1 /
C
C     DATA RIGHT(1),RIGHT(2) /   9344,      0 /
C     DATA RIGHT(3),RIGHT(4) /      0,      0 /
C
C     DATA DIVER(1),DIVER(2) /   9472,      0 /
C     DATA DIVER(3),DIVER(4) /      0,      0 /
C
C     DATA LOG10(1),LOG10(2) /  16282,   8346 /
C     DATA LOG10(3),LOG10(4) / -31493, -12296 /
C
C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
C     DATA SMALL(3),SMALL(4) / O000000, O000000 /
C
C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
C     DATA LARGE(3),LARGE(4) / O177777, O177777 /
C
C     DATA RIGHT(1),RIGHT(2) / O022200, O000000 /
C     DATA RIGHT(3),RIGHT(4) / O000000, O000000 /
C
C     DATA DIVER(1),DIVER(2) / O022400, O000000 /
C     DATA DIVER(3),DIVER(4) / O000000, O000000 /
C
C     DATA LOG10(1),LOG10(2) / O037632, O020232 /
C     DATA LOG10(3),LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C     DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
C
C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572/
C
C
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805665541 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB /
C
C***FIRST EXECUTABLE STATEMENT  D1MACH
C
      D1MACH = DMACH(I)
      RETURN
C
      END
      REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE  R1MACH
C***DATE WRITTEN   790101   (YYMMDD)
C***REVISION DATE  820801   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C     R1MACH 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          A = R1MACH(I)
C
C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF A ABOVE IS
C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C  SINGLE-PRECISION MACHINE CONSTANTS
C  R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C  R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C  R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C  R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C  R1MACH(5) = LOG10(B)
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C     DATA RMACH(1) / Z400800000 /
C     DATA RMACH(2) / Z5FFFFFFFF /
C     DATA RMACH(3) / Z4E9800000 /
C     DATA RMACH(4) / Z4EA800000 /
C     DATA RMACH(5) / Z500E730E8 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
C
C     DATA RMACH(1) / O1771000000000000 /
C     DATA RMACH(2) / O0777777777777777 /
C     DATA RMACH(3) / O1311000000000000 /
C     DATA RMACH(4) / O1301000000000000 /
C     DATA RMACH(5) / O1157163034761675 /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C     DATA RMACH(1) / 00014000000000000000B /
C     DATA RMACH(2) / 37767777777777777777B /
C     DATA RMACH(3) / 16404000000000000000B /
C     DATA RMACH(4) / 16414000000000000000B /
C     DATA RMACH(5) / 17164642023241175720B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1
C
C     DATA RMACH(1) / 200004000000000000000B /
C     DATA RMACH(2) / 577777777777777777777B /
C     DATA RMACH(3) / 377214000000000000000B /
C     DATA RMACH(4) / 377224000000000000000B /
C     DATA RMACH(5) / 377774642023241175720B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC RMACH(5)
C
C     DATA SMALL/20K,0/,LARGE/77777K,177777K/
C     DATA RIGHT/35420K,0/,DIVER/36020K,0/
C     DATA LOG10/40423K,42023K/
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
C     DATA LARGE(1),LARGE(2) / [37777777, [00000177 /
C     DATA RIGHT(1),RIGHT(2) / [20000000, [00000352 /
C     DATA DIVER(1),DIVER(2) / [20000000, [00000353 /
C     DATA LOG10(1),LOG10(2) / [23210115, [00000377 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
C
C     DATA RMACH(1) / O402400000000 /
C     DATA RMACH(2) / O376777777777 /
C     DATA RMACH(3) / O714400000000 /
C     DATA RMACH(4) / O716400000000 /
C     DATA RMACH(5) / O776464202324 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C
C     3 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE91), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
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
C     DATA RMACH(1) / Z00100000 /
C     DATA RMACH(2) / Z7FFFFFFF /
C     DATA RMACH(3) / Z3B100000 /
C     DATA RMACH(4) / Z3C100000 /
C     DATA RMACH(5) / Z41134413 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
C
C     DATA RMACH(1) / "000400000000 /
C     DATA RMACH(2) / "377777777777 /
C     DATA RMACH(3) / "146400000000 /
C     DATA RMACH(4) / "147400000000 /
C     DATA RMACH(5) / "177464202324 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1) /    8388608 /
C     DATA LARGE(1) / 2147483647 /
C     DATA RIGHT(1) /  880803840 /
C     DATA DIVER(1) /  889192448 /
C     DATA LOG10(1) / 1067065499 /
C
C     DATA RMACH(1) / O00040000000 /
C     DATA RMACH(2) / O17777777777 /
C     DATA RMACH(3) / O06440000000 /
C     DATA RMACH(4) / O06500000000 /
C     DATA RMACH(5) / O07746420233 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1),SMALL(2) /   128,     0 /
C     DATA LARGE(1),LARGE(2) / 32767,    -1 /
C     DATA RIGHT(1),RIGHT(2) / 13440,     0 /
C     DATA DIVER(1),DIVER(2) / 13568,     0 /
C     DATA LOG10(1),LOG10(2) / 16282,  8347 /
C
C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
C     DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
C     DATA DIVER(1),DIVER(2) / O032400, O000000 /
C     DATA LOG10(1),LOG10(2) / O037632, O020233 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
C
C     DATA RMACH(1) / O000400000000 /
C     DATA RMACH(2) / O377777777777 /
C     DATA RMACH(3) / O146400000000 /
C     DATA RMACH(4) / O147400000000 /
C     DATA RMACH(5) / O177464202324 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C    (EXPRESSED IN INTEGER AND HEXADECIMAL)
C  ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
C  *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C     DATA SMALL(1) /       128 /
C     DATA LARGE(1) /    -32769 /
C     DATA RIGHT(1) /     13440 /
C     DATA DIVER(1) /     13568 /
C     DATA LOG10(1) / 547045274 /
C
C     DATA SMALL(1) / Z00000080 /
C     DATA LARGE(1) / ZFFFF7FFF /
C     DATA RIGHT(1) / Z00003480 /
C     DATA DIVER(1) / Z00003500 /
C     DATA LOG10(1) / Z209B3F9A /
C
C***FIRST EXECUTABLE STATEMENT  R1MACH
C
      R1MACH = RMACH(I)
      RETURN
C
      END
    3   DL1BL3= BLISS VOL II P.294... ALL DATA
   60 1  5.0
  (1X,F6.3,F3.0,2F6.3)
    125  1   300   402
    164  1   156   258
   -065  1   401   503
    120  1   137   239
    034  1   214   316
    010  1   084   186
   -079  1   137   239
   -533  1   031   133  X
   -611  1   227   329  X
    076  1   092   194
    319  1   427   228
    465  1   487   288
    313  1   487   288
    193  1   487   288
    129  1   589   390
    209  1   571   372
    284  1   446   247
    173  1   570   371
    264  1   472   273
    132  1   476   277
    324  1   696   321
    367  1   729   354
    321  1   509   134
    375  1   559   184
    345  1   679   304
    341  1   583   208
    327  1   742   367
    256  1   781   406
    256  1   739   364
    214  1   865   490
    501  1   723   223
    318  1   940   440
    317  1   903   403
    349  1   910   410
    402  1   684   184
    374  1   904   404
    340  1   887   387
    598  1   593   093
    444  1   640   140
    543  1   512   012
    559  1   832   156
    705  1   817   141
    585  1   881   205
    588  1   848   172
    532  1   857   181
    593  1   808   132
    593  1   829   153
    559  1   872   196
    611  1   805   129
    561  1   904   228
    561  1  1047   246
    733  1   942   141
    553  1  1194   393
    593  1  1090   289
    525  1  1000   199
    580  1  1054   253
    711  1   956   155
    546  1  1101   300
    632  1  1060   259
    733  1   893   092
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
    3   DL1BL3= BLISS VOL II P.294... ALL DATA
        PROBLEM TITLE   DL1BL3= BLISS VOL II P.294... ALL DATA
    NUMBER OF OBSERVATION=   60
    NUMBER OF PARAMETERS=    3
    MINIMUM NUMBER OF PARAMETERS CONSIDERED=    1
    PERCENTAGE DEVIATION FROM OPTIMALITY ALLOWED=  5.00
        **BEST SUBSET LAV PROGRAM
 IFAULT=  0
 
 
   BEST RESULTS FOR LAV SUBSET OF SIZE=  1
    **************SUM OF ABSOLUTE VALUES=          7.659
 
 BETA(  2)           .550
 
 
   BEST RESULTS FOR LAV SUBSET OF SIZE=  2
    **************SUM OF ABSOLUTE VALUES=          5.407
 
 BETA(  2)           .818
 
 BETA(  3)          -.782
 
 
   BEST RESULTS FOR LAV SUBSET OF SIZE=  3
    **************SUM OF ABSOLUTE VALUES=          3.877
 
 BETA(  3)         -1.116
 
 BETA(  2)           .628
 
 BETA(  1)           .235
 
 
     ITERATION COUNT=    106
    0