home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / CLIPPER / MISC / ACM615.ZIP / ACM615.FOR
Encoding:
Text File  |  1986-07-23  |  53.5 KB  |  1,819 lines
  1. C     ALGORITHM 615, COLLECTED ALGORITHMS FROM ACM.
  2. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 2,
  3. C     JUN., 1984, P. 202-206.
  4. C      PROGRAM SUBL1 (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
  5. C****************************************************************
  6. C      *THIS CODE SOLVES THE L1 BEST SUBSET PROBLEM*
  7. C      *SUBMITTED TO ACM TRANS.SPRING 1982*
  8. C       *LAST REVISED MAR. 12, 1984*
  9. C****************************************************************
  10.       DIMENSION X(300,20), Y(300)
  11.       DIMENSION ZL(20), BVAL(210), IDEX(210), ISTAT(20)
  12. C
  13.       COMMON /BLOCK/
  14.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  15.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  16.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  17. C
  18.       INTEGER  I,IBASE,IDEX,IFAULT,IK,IK1,IMT,INCOL,INDEX,INPROB,IPAR
  19.       INTEGER  ISTAT,ITER,J,JMIN,K,L,LABEL,LEVEL,M,MININ,MMAX
  20.       INTEGER  N,NAME,NMAX,NPROB
  21. C
  22.       DOUBLE PRECISION ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
  23. C++   REAL             ACU,BIG,BSAVE,BVAL,LU,PIE,POPT,SAD,SADT,SAVTOT
  24.       DOUBLE PRECISION SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
  25. C++   REAL             SIG,SIGMA,X,Y,ZL,ZLOW,ZSAVE
  26. C
  27. C
  28.       LOGICAL DIRECT,INTL
  29.       INTEGER UNITGO,UNITIN
  30. C
  31. C-----UNITIN    IS THE INPUT UNIT NUMBER
  32. C-----UNITGO    IS THE OUTPUT UNIT NUMBER
  33. C
  34.       DIMENSION IMT(72), NAME(40)
  35. C
  36. C-----IMT       IS USED TO STORE THE INPUT FORMAT
  37. C-----NAME      SAVES THE PROBLEM NAME
  38. C
  39.       UNITIN=5
  40.       UNITGO=6
  41.       NMAX=300
  42.       MMAX=30
  43.       NPROB=0
  44. C
  45. C-----K        STORES THE NUMBER OF PARAMETERS IN THE FULL MODEL
  46. C
  47.    10 READ (UNITIN,50,END=45) K,NAME
  48.       WRITE(UNITGO,50)K,NAME
  49.       IF (K.EQ.0) STOP
  50.       WRITE (UNITGO,60) NAME
  51. C
  52. C-----N          STORES THE NUMBER OF OBSERVATIONS
  53. C-----MININ      STORES THE MINIMUM NUMBER OF PARAMETERS CONSIDERED
  54. C-----POPT       DEVIATION FROM OPTIMALITY ALLOWED
  55. C
  56.       READ (UNITIN,130) N,MININ,POPT
  57.    51 FORMAT (4X,22HNUMBER OF OBSERVATION=,I5/4X,21HNUMBER OF PARAMETERS
  58.      1=,I5/4X,40HMINIMUM NUMBER OF PARAMETERS CONSIDERED=,I5/4X,45HPERCE
  59.      2NTAGE DEVIATION FROM OPTIMALITY ALLOWED=,F6.2/8X,2H**,
  60.      3 23HBEST SUBSET LAV PROGRAM)
  61.       WRITE (UNITGO,51) N,K,MININ,POPT
  62.       READ (UNITIN,150) IMT
  63.       DO 20 I=1,N
  64.          READ (UNITIN,IMT) Y(I), (X(I,J),J=1,K)
  65. C
  66. C        WRITE (UNITGO,IMT) Y(I), (X(I,J),J=1,K)
  67. C
  68.    20 CONTINUE
  69. C
  70.       READ (UNITIN,70) (ISTAT(I),I=1,K)
  71. C
  72.       ITER=0
  73. C
  74. C******CALL THE ROUTINE TO FIND THE BEST SUBSETS
  75. C
  76. C
  77.       CALL KBEST (X,Y,K,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL,IDEX,
  78.      1ISTAT,ZL)
  79. C
  80.       WRITE (UNITGO,80) IFAULT
  81. C
  82. C     WRITE THE FINAL BEST SUBSET SOLUTION
  83. C
  84.       J=K*(K+1)/2
  85.       JMIN=(MININ-1)*(MININ)/2
  86.       J=J-JMIN
  87.       DO 40 I=MININ,K
  88.          WRITE (UNITGO,110) I
  89.          WRITE (UNITGO,100) ZL(I)
  90.          DO 30 L=1,I
  91.             WRITE (UNITGO,90) IDEX(J),BVAL(J)
  92.             J=J-1
  93.    30    CONTINUE
  94.    40 CONTINUE
  95. C
  96.           NPROB=NPROB+1
  97.           WRITE (UNITGO,160)ITER
  98.       GO TO 10
  99. C
  100.    45 STOP
  101. C
  102.    50 FORMAT (I5,3X,40A1)
  103.    60 FORMAT (3X,5H     ,14HPROBLEM TITLE ,2X,40A1)
  104.    70 FORMAT (10I2)
  105.    80 FORMAT (8H IFAULT=,I3)
  106.    90 FORMAT (3X/6H BETA(,I3,1H),F15.3)
  107.   100 FORMAT (4X,14(1H*),23HSUM OF ABSOLUTE VALUES=,F15.3)
  108.   110 FORMAT (//3X,36HBEST RESULTS FOR LAV SUBSET OF SIZE=,I3)
  109.   130 FORMAT (I5,I2,F6.2)
  110.   150 FORMAT (72A1)
  111.   160 FORMAT (//5X,16HITERATION COUNT=,I7)
  112. C
  113.       END
  114.       SUBROUTINE KBEST (X,Y,M,N,ITER,IFAULT,POPT,MININ,NMAX,MMAX,BVAL
  115.      1,IDEX,ISTAT,ZL)
  116. C
  117. C***********************************************************************
  118. C
  119. C
  120. C     THE PURPOSE OF THIS PROGRAM IS TO DETERMINE THE BEST SUBSET OF
  121. C     PARAMETERS TO FIT A LINEAR REGRESSION UNDER AN LEAST ABSOLUTE
  122. C     VALUE CRITERION. THIS PROGRAM UTILIZES THE SIMPLEX METHOD OF
  123. C     LINEAR PROGRAMMING WITHIN A BRANCH-AND-BOUND ALGORITHM TO
  124. C     SOLVE THE BEST SUBSET PROBLEM.
  125. C
  126. C
  127. C     THE ALGORITHM IS BASED ON THE PUBLICATION:
  128. C        ARMSTRONG,R.D. AND M.T. KUNG "AN ALGORITHM TO SELECT THE BEST
  129. C        SUBSET FOR A LEAST ABSOLUTE VALUE REGRESSION PROBLEM,"
  130. C        OPTIMIZATION IN STATISTICS, TIMS STUDIES OF THE MANAGEMENT
  131. C        SCIENCES.
  132. C
  133. C   FORMAL PARAMETERS
  134. C
  135. C   X      REAL ARRAY        INPUT:  VALUES OF INDEPENDENT VARIABLES
  136. C          (NMAX,MMAX)                SUCH THAT EACH ROW CORREPSONDS TO
  137. C                                     AN OBSERVATION
  138. C
  139. C   Y      REAL ARRAY        INPUT:  VALUES OF THE DEPENDENT VARIABLES
  140. C            (NMAX)
  141. C
  142. C   M      INTEGER           INPUT:  NUMBER OF DEPENDENT VARIABLES
  143. C
  144. C   N      INTEGER           INPUT:  NUMBER OF OBSERVATIONS
  145. C
  146. C   ITER   INTEGER           OUTPUT: NUMBER OF ITERATIONS
  147. C
  148. C   IFAULT INTEGER           OUTPUT: FAILURE INDICATOR
  149. C                              =0     NORMAL TERMINATION
  150. C                              =1     OBSERVATION MATRIX DOES NOT HAVE
  151. C                                     FULL ROW RANK (RANK M)
  152. C                              =2     PROBLEM SIZE OUT OF RANGE
  153. C                              =3     NO PIVOT ELEMENT FOUND IMPLYING
  154. C                                     NEAR SINGULAR BASIS
  155. C
  156. C   POPT    REAL             INPUT:  PERCENTAGE DEVIATION FROM
  157. C                                    OPTIMALITY ALLOWED
  158. C
  159. C   MININ   INTEGER          INPUT:  MINIMUM NUMBER OF PARAMETERS IN THE
  160. C                                    MODEL.  BEST SUBSET OF SIZE MININ T
  161. C                                    M IS OBTAINED.
  162. C
  163. C   NMAX    INTEGER          INPUT:  DIMENSION OF ROWS IN X (ALSO Y)
  164. C
  165. C   MMAX    INTEGER          INPUT:  DIMENSION OF COLUMNS IN X
  166. C
  167. C   BVAL    REAL ARRAY       OUTPUT: ARRAY OF OPTIMAL BETA VALUES FOR
  168. C                                    EACH SUBSET. THE BETA VALUES FOR
  169. C                                    THE SUBSET OF SIZE M ARE STORED
  170. C                                    IN POSITIONS BVAL(1),BVAL(2),...,
  171. C                                    BVAL(M), FOR THE SUBSET OF SIZE
  172. C                                    M-1 THE VALUES ARE STORED IN
  173. C                                    POSITIONS BVAL(M+1),BVAL(M+2),...,
  174. C                                    BVAL(2M-1). IN GENERAL, THE BETA
  175. C                                    VALUES FOR THE OPTIMAL SUBSET OF
  176. C                                    SIZE K ARE STORED IN POSITIONS
  177. C                                    BVAL(L),...,BVAL(L-K+1) WHERE
  178. C                                       L=(M*(M+1)-K*(K+1))/2 + 1
  179. C
  180. C  IDEX     INTEGER ARRAY    OUTPUT: BETA INDEX SET FOR THE OPTIMAL
  181. C        (((MMAX+1)*MMAX)/2)         SUBSET. THIS ARRAY IS A  PARALLEL
  182. C                                    ARRAY FOR BVAL; I.E., IF BVAL(J)=2.
  183. C                                    AND IDEX(J)=7  THEN BETA(7)=2.7 IN
  184. C                                    THE ASSOCIATED OPTIMAL SUBSET.
  185. C
  186. C  ISTAT    INTEGER ARRAY    INPUT:  PARAMETER STATUS ARRAY.
  187. C            (MMAX)
  188. C
  189. C                                        1 IF BETA(J) IS REQUIRED
  190. C                                          IN EVERY MODEL
  191. C                     ISTAT(J)=
  192. C
  193. C                                        0 OTHERWISE
  194. C
  195. C  ZL       REAL ARRAY       OUTPUT: BEST OBJECTIVE VALUE FOR EACH SUBSE
  196. C            (MMAX)
  197. C                                    ZL(J) GIVES THE BEST OBJECTIVE VALU
  198. C                                    FOR THE SUBSET WITH M-J+1 PARAMETER
  199. C***********************************************************************
  200. C
  201. C    IMPLEMENTATION NOTES:
  202. C
  203. C        1. THE ROUTINE USES TWO MACHINE DEPENDENT VALUES ACU AND BIG.
  204. C           (THESE ARE SET AT THE BEGINNING OF THIS SUBROUTINE)
  205. C           A) ACU IS USED TO TEST FOR "ZERO". ACU SHOULD BE SET TO
  206. C              APPROXIMATELY 100 * THE RELATIVE MACHINE ACCURACY OF THE
  207. C              SYSTEM IN USE.
  208. C           B) BIG IS USED TO INITIALIZE THE ARRAY ZL (DESCRIBED ABOVE).
  209. C              BIG SHOULD BE SET TO THE LARGEST FLOATING POINT VALUE
  210. C              ASSIGNABLE.
  211. C
  212. C        2. BOTH SINGLE AND DOUBLE PRECISION VERSIONS ARE SUPPLIED. THIS
  213. C           VERSION IS IN DOUBLE PRECISION. TO GET A SINGLE PRECISION CO
  214. C           ALL STATEMENTS HAVING A "C++" IN COLUMNS 1-3 SHOULD BE MODIF
  215. C           BY DELETING THE "C++". THE CORRESPONDING DOUBLE PRECISION
  216. C           STATEMENTS SHOULD BE EITHER DELETED OR COMMENTED OUT.
  217. C
  218. C        3. ARRAY DIMENSIONS
  219. C           THE CODE IS CURRENTLY DIMENSIONED TO SOLVE PROBLEMS WITH UP
  220. C           20 PARAMETERS AND 300 OBSERVATIONS. THE DIMENSION SIZES ARE
  221. C           DETERMINED AS FOLLOWS
  222. C               20 MMAX
  223. C              300 NMAX
  224. C              210 (NMAX+1)*MMAX/2
  225. C             6000 NMAX*MMAX
  226. C               10 MAXIMUM VALUE OF IK .
  227. C
  228. C            IF THE DIMENSION SIZES ARE CHANGED THE COMMON AND DIMENSION
  229. C            STATEMENTS WILL NEED TO BE MODIFIED IN EACH SUBROUTINE.
  230. C
  231. C
  232. C***********************************************************************
  233. C
  234. C     SUBROUTINES:
  235. C
  236. C         CALBET   - BACK-SOLVES A SYSYTEM OF EQUATIONS
  237. C
  238. C         CALCPI   - FORWARD-SOLVES A SYSYTEM OF EQUATIONS
  239. C
  240. C         KBEST    - THE DRIVER
  241. C
  242. C         L1NORM   - SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
  243. C                    PARAMETERS INCLUDED IN THE MODEL
  244. C
  245. C         PHASE2   - SOLVES THE CURRENT REGRESSION PROBLEM WITH A PRIMAL
  246. C                    ALGORITHM
  247. C
  248. C         SETUP    - DETERMINES THE FORM OF THE CURRENT PROBLEM BY
  249. C                    CHOOSING THE PARAMETER TO LEAVE BASED ON A PENALTY
  250. C
  251. C         UPDATE   - UPDATES LU DECOMPOSITION MATRIX
  252. C
  253. C
  254. C***********************************************************************
  255. C
  256. C
  257. C     DESCRIPTION OF VARIABLES:
  258. C
  259. C     IK:LENGTH OF CANDIDATE LIST
  260. C     BETA:ESTIMATES TO THE CURRENT SUBPROBLEM
  261. C     SAD:THE MINIMUM TOTAL ABSOLUTE DEVIATION
  262. C     IBASE:THE INDEX ARRAY OF COLUMNS OF THE BASIS
  263. C
  264. C        N,M,NMAX,MMAX,X,Y, ARE NOT CHANGED IN THE SUBROUTINE;
  265. C          SAD IS UPDATED AT EACH ITERATION
  266. C
  267. C     LU:THE LU DECOMPOSITION OF THE CURRENT BASIS
  268. C     INDEX:THE INDEX ARRAY OF ROWS OF LU
  269. C     TOT:THE CURRENT RHS OF THE DUAL PROBLEM
  270. C     SIGMA:INDICATOR ARRAY TO SPECIFY WHETHER A NONBASIC DUAL VARIABLE
  271. C           IS AT UPPER OR LOWER BOUND(+1 IMPLIES UPPER;-1 IMPLIES
  272. C           LOWER)
  273. C     INEXT:A LOCAL ARRAY FOR THE SORT ROUTINE
  274. C     RHS: A LOCAL ARRAY FOR THE CALBET ROUTINE
  275. C
  276. C     IPAR(I)=-K IF THE K-TH PARAMETER (BETA) IS OUT OF MODEL AT LEVEL I
  277. C         = K IF THE K-TH PARAMETER IS IN
  278. C
  279. C     LABEL(I)  SAVES THE INDICES OF THE BASIC OBSERVATIONS IN THE
  280. C            PREDECESSOR PATH
  281. C     ZSAVE(I)  SAVES THE OBJECTIVE VALUES ON THE PREDECESSOR PATH
  282. C
  283. C**********************************************************************
  284. C
  285.       DIMENSION X(300,20), Y(300), ZL(20), BVAL(210), IDEX(210), ISTAT(2
  286.      10)
  287. C
  288.       DIMENSION BETA(20), TOT(20), PI(20), REPP(20)
  289. C
  290.       COMMON /BLOCK/
  291.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  292.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  293.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  294. C
  295.       INTEGER I,IBASE,IDEX,IFAULT,IK,IK1,INCOL,INDEX,INPROB,IPAR
  296.       INTEGER ISTAT,ITER,J,JINM,K2,K3,KKK1,L,LABEL,LEVEL,M,MININ
  297.       INTEGER MMAX,N,NFREE,NMAX,NUMIN
  298. C
  299.       DOUBLE PRECISION ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
  300. C++   REAL             ACU,BETA,BIG,BSAVE,BVAL,LU,PI,PIE,POPT,POPT1
  301.       DOUBLE PRECISION REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
  302. C++   REAL             REPP,SAD,SADT,SAVTOT,SIG,SIGMA,TOT,X,Y
  303.       DOUBLE PRECISION ZERO,ZL,ZLOW,ZSAVE
  304. C++   REAL             ZERO,ZL,ZLOW,ZSAVE
  305. C
  306.       LOGICAL DIRECT,INTL
  307. C
  308.       DATA ZERO/0.0D0/
  309. C++   DATA ZERO/0.0E0/
  310. C
  311. C     ASSIGN VALUES TO MACHINE DEPENDENT CONSTANTS
  312. C    SEE - FOX, P.A., A.D. HALL AND N.L. SCHRYER, "FRAMEWORK FOR
  313. C          A PORTABLE LIBRARY: ALGORITHM 528," ACM TRANSACTIONS
  314. C          ON MATHEMATICAL SOFTWARE, VOL 4, NO 2, JUNE 1978.
  315. C
  316.        ACU = D1MACH(1)
  317.        ACU = 1000.*ACU
  318.        BIG = D1MACH(2)
  319. C++    ACU = R1MACH(1)
  320. C++    ACU = 100.*ACU
  321. C++    BIG = R1MACH(2)
  322. C
  323. C
  324.       IFAULT=0
  325.       IK=5
  326.       IK1=IK-1
  327.       POPT1=(100.-POPT)/100.
  328. C
  329. C     TEST PROBLEM SIZE
  330. C
  331.       IF (N.LE.NMAX.AND.M.LE.MMAX) GO TO 10
  332.       IFAULT=2
  333.       RETURN
  334.    10 CONTINUE
  335.       DO 20 I=1,M
  336.          INDEX(I)=I
  337.          TOT(I)=ZERO
  338.    20 INCOL(I)=I
  339.       CALL L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
  340.       IF (IFAULT.GT.0) RETURN
  341. C
  342. C     INITIALIZATION
  343. C
  344. C     SAVE THE INITIAL SOLUTION
  345. C
  346.       JINM=0
  347.       DO 30 I=1,M
  348.          BSAVE(I)=BETA(I)
  349.          BVAL(I)=BETA(I)
  350.          PIE(I)=PI(I)
  351.          LABEL(I)=IBASE(I)
  352.          IDEX(I)=I
  353.          INPROB(I)=I
  354.          SAVTOT(I)=TOT(I)
  355.          ZL(I)=BIG
  356.          JINM=JINM+ISTAT(I)
  357.    30 CONTINUE
  358.       DO 40 I=1,N
  359.          SIG(I)=SIGMA(I)
  360.    40 CONTINUE
  361. C
  362.       LEVEL=0
  363.       DIRECT=.TRUE.
  364. C
  365.       IF (JINM.GT.MININ) MININ=JINM
  366. C
  367. C     INITIALIZE NUMBER OF FREE PARAMETERS IN MODEL
  368. C
  369.       NFREE=M-JINM
  370. C
  371.       ZL(M)=SAD
  372.       ZSAVE(M)=SAD
  373. C
  374.       IF (MININ.EQ.M) RETURN
  375. C
  376. C       NUMIN GIVES THE NUMBER OF PARAMETERS IN THE MODEL
  377. C
  378.       NUMIN=M
  379. C
  380. C     START THE MAIN LOOP
  381. C
  382.    50 NUMIN=NUMIN-1
  383.       LEVEL=LEVEL+1
  384.       IF (NUMIN.LT.MININ.OR.NFREE.EQ.0) GO TO 120
  385.       K2=(M*(M+1)-(NUMIN+1)*(NUMIN+2))/2
  386.       KKK1=NUMIN+1
  387.       ZLOW=ZL(NUMIN)*POPT1
  388.       IF (DIRECT) GO TO 80
  389.       DO 60 I=1,KKK1
  390.          J=K2+I
  391.          INDEX(I)=I
  392.          IBASE(I)=LABEL(J)
  393.          INCOL(I)=INPROB(J)
  394.          PI(I)=PIE(J)
  395.          BETA(I)=BSAVE(J)
  396.          TOT(I)=SAVTOT(J)
  397.    60 CONTINUE
  398. C
  399. C     LOAD IN THE BOUNDS FROM THE IMMEDIATE PREDECESSOR
  400. C
  401.       K3=(M-KKK1)*N
  402.       DO 70 I=1,N
  403.          K3=K3+1
  404.          SIGMA(I)=SIG(K3)
  405.    70 CONTINUE
  406.       SAD=ZSAVE(KKK1)
  407. C
  408. C     SET UP A NEW PROBLEM
  409. C
  410.    80 CALL SETUP (X,KKK1,IFAULT,ISTAT,BETA,TOT,PI,REPP)
  411.       IF (IFAULT.GT.0) RETURN
  412.       NFREE=NFREE-1
  413. C
  414. C     CALL PRIMAL L.P. CODE
  415. C
  416.       CALL PHASE2 (X,Y,NUMIN,N,ITER,IFAULT,BETA,TOT,PI,REPP)
  417. C
  418. C      SAVE THE SOLUTION DATA FOR LATER RECALL
  419. C
  420.       K2=K2+KKK1+1
  421.       J=K2
  422.       DO 90 I=1,NUMIN
  423.          L=K2+I
  424.          PIE(L)=PI(I)
  425.          LABEL(L)=IBASE(I)
  426.          BSAVE(L)=BETA(I)
  427.          INPROB(L)=INCOL(I)
  428.          SAVTOT(L)=TOT(I)
  429.    90 CONTINUE
  430. C
  431. C     SAVE THE OBJECTIVE VALUE
  432. C
  433.       ZSAVE(NUMIN)=SAD
  434. C
  435. C     SAVE THE NONBASIC BOUNDS
  436. C
  437.       K3=(M-NUMIN)*N
  438.       DO 100 I=1,N
  439.          K3=K3+1
  440.          SIG(K3)=SIGMA(I)
  441.   100 CONTINUE
  442. C
  443.       DIRECT=.TRUE.
  444. C
  445. C     CHECK FOR A BETTER SOLUTION
  446. C
  447.       IF (SADT.GE.ZLOW) GO TO 50
  448. C
  449. C     SAVE THE BETTER SOLUTION
  450. C
  451.       ZL(NUMIN)=SAD
  452.       K2=J
  453.       DO 110 I=1,NUMIN
  454.          L=K2+I
  455.          IDEX(L)=INCOL(I)
  456.          BVAL(L)=BETA(I)
  457.   110 CONTINUE
  458. C
  459. C     GO TO TOP OF LOOP
  460. C
  461.       GO TO 50
  462.   120 NUMIN=NUMIN+1
  463. C
  464.       DIRECT=.FALSE.
  465. C
  466. C     MUST WORK BACK UP THE TREE
  467. C
  468.       LEVEL=LEVEL-1
  469.   130 IF (IPAR(LEVEL).GT.0) GO TO 140
  470.       IPAR(LEVEL)=-IPAR(LEVEL)
  471.       J=IPAR(LEVEL)
  472.       ISTAT(J)=1
  473.       NUMIN=NUMIN+1
  474.       GO TO 50
  475.   140 J=IPAR(LEVEL)
  476.       ISTAT(J)=0
  477.       NFREE=NFREE+1
  478.       LEVEL=LEVEL-1
  479.       IF (LEVEL.GT.0) GO TO 130
  480.       RETURN
  481. C
  482.       END
  483.       SUBROUTINE SETUP (X,M,IFAULT,ISTAT,BETA,TOT,PI,REPP)
  484. C
  485. C     THIS SUBROUTINE DETERMINES THE FORM OF THE CURRENT SUBPROBLEM
  486. C
  487.       DIMENSION X(300,20), ISTAT(20), RHS(20), BETA(20), TOT(20), PI(20)
  488.      1, REPP(20)
  489. C
  490.       COMMON /BLOCK/
  491.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  492.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  493.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  494. C
  495.       INTEGER I,IBASE,IFAULT,IFIXI,II,IK,IK1,INCOL,INDEX,INPROB
  496.       INTEGER IOUT,IPAR,ISAVE,ISTAT,KKK,L,LABEL,LEAVE,LEVEL,LL,M
  497. C
  498.       DOUBLE PRECISION ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
  499. C++   REAL             ACU,BESIDE,BETA,BIG,BSAVE,LU,ONE,PENALT,PI,PIE
  500.       DOUBLE PRECISION RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
  501. C++   REAL             RATIO,REPP,RHO,RHS,SAD,SADT,SAVTOT,SIDE,SIG
  502.       DOUBLE PRECISION SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
  503. C++   REAL             SIGMA,SREPP,TEST,TOT,X,ZERO,ZLOW,ZSAVE
  504. C
  505.       LOGICAL DIRECT,INTL
  506. C
  507.       DATA ONE/1.0D0/,ZERO/0.0D0/
  508. C++   DATA ONE/1.0E0/,ZERO/0.0E0/
  509. C
  510. C     IF THE IMMEDIATE PREDECESSOR IS IN MEMORY GO ON
  511. C
  512.       IF (DIRECT) GO TO 10
  513. C
  514. C     RECONSTRUCT THE LU DECOMPOSITION
  515. C
  516.       KKK=1
  517.       CALL UPDATE (KKK,IFAULT,X,1,M)
  518.       IF (IFAULT.GT.0) RETURN
  519. C
  520. C     RECALCULATE THE BASIC PI VALUES
  521. C
  522.       KKK=0
  523. C
  524.       CALL CALCPI (KKK,PI,TOT,M)
  525. C
  526. C     DETERMINE PARAMETER TO LEAVE BASED ON PENALTY
  527. C
  528.    10 PENALT=BIG
  529.       KKK=0
  530.       DO 20 L=1,M
  531.          RHS(L)=ZERO
  532.    20 CONTINUE
  533. C
  534.       DO 40 I=1,M
  535.          II=INDEX(I)
  536.          LL=INCOL(II)
  537.          IF (ISTAT(LL).EQ.1) GO TO 40
  538. C++      RHO=SIGN(1.,BETA(II))
  539.          RHO=DSIGN(1.D0,BETA(II))
  540.          RHS(II)=ONE
  541.          KKK=I-1
  542.          CALL CALCPI (KKK,REPP,RHS,M)
  543.          RHS(II)=ZERO
  544.          TEST=BIG
  545.          DO 30 L=1,M
  546. C++         IF ( ABS(REPP(L)).LE.ACU) GO TO 30
  547.             IF (DABS(REPP(L)).LE.ACU) GO TO 30
  548. C++         SREPP=SIGN(1.,REPP(L))
  549.             SREPP=DSIGN(1.D0,REPP(L))
  550. C++         RATIO=ABS((1.-RHO*SREPP*PI(L))/REPP(L))
  551.             RATIO=DABS((1.-RHO*SREPP*PI(L))/REPP(L))
  552.             IF (RATIO.GE.TEST) GO TO 30
  553.             SIDE=RHO*SREPP
  554.             TEST=RATIO
  555.             IOUT=L
  556.    30    CONTINUE
  557. C
  558. C     SEE IF THE PENALTY IS LESS
  559. C
  560. C++      IF (TEST*ABS(BETA(II)).GE.PENALT) GO TO 40
  561.          IF (TEST*DABS(BETA(II)).GE.PENALT) GO TO 40
  562.          LEAVE=IOUT
  563.          IFIXI=II
  564.          BESIDE=SIDE
  565. C++      PENALT=TEST*ABS(BETA(II))
  566.          PENALT=TEST*DABS(BETA(II))
  567.    40 CONTINUE
  568. C
  569. C     UPDATE THE OBJECTIVE VALUE
  570. C
  571.       SAD=SAD+PENALT
  572. C
  573. C     PARAMETER INCOL(IFIXI) WILL LEAVE THE MODEL
  574. C     PI(IBASE(LEAVE)) WILL LEAVE THE BASIS AND GO TO BOUND
  575. C
  576. C     SWITCH UNWANTED PARAMETER AND PI TO THE END OF LIST
  577. C
  578.       ISAVE=INCOL(IFIXI)
  579.       INCOL(IFIXI)=INCOL(M)
  580.       INCOL(M)=ISAVE
  581.       TOT(IFIXI)=TOT(M)
  582. C
  583. C     LABEL THIS NODE AND FLAG THE OUTGOING PARAMETER
  584. C
  585.       IPAR(LEVEL)=-ISAVE
  586.       ISTAT(ISAVE)=-1
  587. C
  588.       ISAVE=IBASE(LEAVE)
  589.       IBASE(LEAVE)=IBASE(M)
  590.       IBASE(M)=ISAVE
  591. C
  592. C     PLACE THE OUTGOING PI AT THE CORRECT BOUND
  593. C
  594.       SIGMA(ISAVE)=BESIDE
  595. C
  596. C     FORM A NEW BASIS
  597. C
  598.       M=M-1
  599. C
  600.       DO 50 I=1,M
  601.    50 INDEX(I)=I
  602. C
  603.       KKK=1
  604.       CALL UPDATE (KKK,IFAULT,X,1,M)
  605.       IF (IFAULT.GT.0) RETURN
  606. C
  607. C     UPDATE THE TOT ARRAY
  608. C
  609.       DO 60 I=1,M
  610.          LL=INCOL(I)
  611.          TOT(I)=TOT(I)-BESIDE*X(ISAVE,LL)
  612.    60 CONTINUE
  613. C
  614.       RETURN
  615. C
  616.       END
  617.       SUBROUTINE PHASE2 (X,Y,M,N,ITER,IFAULT,BETA,TOT,PI,REPP)
  618. C
  619. C     SOLVES LP WITH PRIMAL ALGORITHM
  620. C
  621.       DIMENSION X(300,20), Y(300), BETA(20), TOT(20), PI(20), REPP(20),
  622.      1RHS(20),Z(10),KAN(10)
  623. C
  624.       COMMON /BLOCK/
  625.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  626.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  627.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  628. C
  629.       INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INDEX,INPROB,IOUT,IPAR
  630.       INTEGER ITER,J,K,KAN,KKK,KROW,L,LABEL,LEVEL,LL,M,N
  631. C
  632.       DOUBLE PRECISION ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
  633. C++   REAL             ACU,BEST,BETA,BIG,BSAVE,LU,ONE,PI,PIE,RATIO,REPP
  634.       DOUBLE PRECISION RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
  635. C++   REAL             RHO,RHS,SAD,SADT,SAVTOT,SIG,SIGMA,SREPP,TEST
  636.       DOUBLE PRECISION TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
  637. C++   REAL             TOT,X,XTWO,Y,Z,ZERO,ZC,ZLOW,ZSAVE,ZZ
  638. C
  639.       LOGICAL DIRECT,INTL
  640. C
  641.       DATA ZERO/0.0D0/,ONE/1.0D0/
  642. C++   DATA ZERO/0.0E0/,ONE/1.0E0/
  643. C
  644. C     Z(I) STORES THE REDUCED COST VALUES ON THE LIST
  645. C     KAN(I) STORES THE ROW NUMBER ON THE LIST
  646. C     IK IS THE LENGTH OF THE LIST
  647. C
  648. C     CALCULATE THE SIMPLEX MULTIPLIERS
  649. C
  650.       DO 10 I=1,M
  651.          K=IBASE(I)
  652.    10 PI(I)=Y(K)
  653.       KKK=0
  654.       CALL CALBET (KKK,BETA,PI,M)
  655. C
  656. C     RECALCULATE THE BASIC PI VALUES
  657. C
  658.       KKK=0
  659.       CALL CALCPI (KKK,PI,TOT,M)
  660. C
  661. C     STORE THE OBJECTIVE VALUE USED FOR TERMINATION CHECK
  662. C
  663.       SADT=SAD
  664.       IF (SAD.GE.ZLOW) RETURN
  665. C
  666. C     GENERATE THE CANDIDATE LIST
  667. C
  668.    20 DO 30 I=1,IK
  669.          Z(I)=ACU
  670.    30 CONTINUE
  671. C
  672. C     CALCULATE THE REDUCED COSTS
  673. C
  674.       DO 70 J=1,N
  675.          IF (SIGMA(J).EQ.ZERO) GO TO 70
  676.          ZC=ZERO
  677.          DO 40 I=1,M
  678.             LL=INCOL(I)
  679.    40    ZC=ZC+BETA(I)*X(J,LL)
  680.          ZZ=(ZC-Y(J))*SIGMA(J)
  681.          IF (ZZ.LE.Z(IK)) GO TO 70
  682. C
  683. C     RANK THE VALUES IN Z(I), IN DESCENDING ORDER
  684. C
  685.          DO 60 K=1,IK1
  686.             IF (ZZ.LE.Z(K)) GO TO 60
  687.             I=IK
  688.    50       Z(I)=Z(I-1)
  689.             KAN(I)=KAN(I-1)
  690.             I=I-1
  691.             IF (I.GT.K) GO TO 50
  692.             Z(K)=ZZ
  693.             KAN(K)=J
  694.             GO TO 70
  695.    60    CONTINUE
  696.          Z(IK)=ZZ
  697.          KAN(IK)=J
  698.    70 CONTINUE
  699.       IF (Z(1).LE.ACU) RETURN
  700.       KROW=1
  701.       ZZ=Z(1)
  702.       IIN=KAN(1)
  703.       ZC=ZZ*SIGMA(IIN)
  704.       GO TO 100
  705.    80 KROW=KROW+1
  706.       IF (KROW.GT.IK) GO TO 20
  707.       IF (Z(KROW).LE.ACU) GO TO 20
  708. C
  709. C     DETERMINE THE ENTERING VARIABLE FROM THE CANDIDATE LIST
  710. C
  711.       IIN=KAN(KROW)
  712.       ZC=ZERO
  713.       DO 90 I=1,M
  714.          LL=INCOL(I)
  715.    90 ZC=ZC+BETA(I)*X(IIN,LL)
  716.       ZC=ZC-Y(IIN)
  717.       ZZ=ZC*SIGMA(IIN)
  718.       IF (ZZ.LE.ACU) GO TO 80
  719.   100 BEST=ZZ
  720.       RHO=SIGMA(IIN)
  721. C
  722. C     RHO: SIGN OF THE INCOMING REDUCED COST
  723. C     IIN: INCOMING CANDIDATE
  724. C
  725. C     FIND THE REPRESENTATION OF THE ENTERING VARIABLE
  726. C
  727.       DO 110 I=1,M
  728.          LL=INCOL(I)
  729.          RHS(I)=X(IIN,LL)
  730.   110 CONTINUE
  731.       KKK=0
  732.       CALL CALCPI (KKK,REPP,RHS,M)
  733. C
  734. C     CALCULATE THE MIN RATIO TEST TO DETERMINE THE LEAVING VARIABLE
  735. C
  736.       TEST=2.0
  737.       DO 120 I=1,M
  738. C++      IF (ABS(REPP(I)).LE.ACU) GO TO 120
  739.          IF (DABS(REPP(I)).LE.ACU) GO TO 120
  740. C++      SREPP=SIGN(1.0,REPP(I))
  741.          SREPP=DSIGN(1.0D0,REPP(I))
  742. C++      RATIO=ABS((1.-RHO*SREPP*PI(I))/REPP(I))
  743.          RATIO=DABS((1.-RHO*SREPP*PI(I))/REPP(I))
  744.          IF (RATIO.GE.TEST) GO TO 120
  745.          TEST=RATIO
  746.          IOUT=I
  747.   120 CONTINUE
  748. C
  749. C     PERFORM FATHOMING TEST BEFORE PIVOT
  750. C
  751.       SADT=SAD+TEST*ZZ
  752.       IF (SADT.GE.ZLOW) RETURN
  753.       SAD=SADT
  754. C
  755.       DO 130 I=1,M
  756.   130 PI(I)=PI(I)+TEST*REPP(I)*RHO
  757.       IF (TEST.LT.2.0) GO TO 150
  758. C
  759.       SIGMA(IIN)=-SIGMA(IIN)
  760.       XTWO=SIGMA(IIN)+SIGMA(IIN)
  761.       DO 140 I=1,M
  762.          LL=INCOL(I)
  763.          TOT(I)=TOT(I)-XTWO*X(IIN,LL)
  764.   140 CONTINUE
  765.       GO TO 80
  766.   150 K=IBASE(IOUT)
  767. C++   SIGMA(K)=SIGN(1.0,PI(IOUT))
  768.       SIGMA(K)=DSIGN(1.0D0,PI(IOUT))
  769.       PI(IOUT)=SIGMA(IIN)*(1.0-TEST)
  770.       IBASE(IOUT)=IIN
  771. C
  772.       DO 160 I=1,M
  773.          LL=INCOL(I)
  774.          TOT(I)=TOT(I)+SIGMA(IIN)*X(IIN,LL)-SIGMA(K)*X(K,LL)
  775.   160 CONTINUE
  776. C
  777.       KKK=IOUT
  778.       CALL UPDATE (IOUT,IFAULT,X,1,M)
  779. C
  780.       DO 170 I=1,M
  781.          RHS(I)=ZERO
  782.   170 CONTINUE
  783.       RHS(KKK)=ONE
  784.       KKK=KKK-1
  785.       CALL CALBET (KKK,REPP,RHS,M)
  786.       DO 180 L=1,M
  787.   180 BETA(L)=BETA(L)-ZC*REPP(L)
  788.       SIGMA(IIN)=ZERO
  789.       ITER=ITER+1
  790.       GO TO 80
  791. C
  792.       END
  793.       SUBROUTINE L1NORM (X,Y,N,M,IFAULT,ITER,BETA,PI,TOT)
  794. C
  795. C     THIS ROUTINE  SOLVES THE INITIAL REGRESSION PROBLEM WITH ALL
  796. C     THE PARAMETERS INCLUDED IN THE MODEL
  797. C
  798.       DIMENSION X(300,20), Y(300), D(300), DELTA(300), INEXT(300), PRICE
  799.      1(300), RHS(20), BETA(20), TOT(20), PI(20)
  800. C
  801. C
  802.       COMMON /BLOCK/
  803.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  804.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  805.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  806. C
  807.       INTEGER I,IBASE,IFAULT,IIN,IK,IK1,INCOL,INEXT,INPROB,IOUT,IPAR
  808.       INTEGER INDEX,IPT,ITER,J,K,K1,KKK,KOUNT,L,LABEL,LEVEL,M,M1,N
  809. C
  810.       DOUBLE PRECISION ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
  811. C++   REAL             ACU,AHALF,AONE,BETA,BIG,BONE,BSAVE,D,DELTA,LU
  812.       DOUBLE PRECISION ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
  813. C++   REAL             ONE,PI,PIE,PRICE,RATIO,RHO,RHS,SAD,SADT,SAVTOT
  814.       DOUBLE PRECISION SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
  815. C++   REAL             SIG,SIGMA,SUBT,SUM,T,TEST,TOT,VAL,X,Y,ZERO,ZLOW
  816.       DOUBLE PRECISION ZSAVE,ZZZ
  817. C++   REAL             ZSAVE,ZZZ
  818. C
  819.       LOGICAL DIRECT,INTL
  820. C
  821.       DATA ONE/1.0D0/,ZERO/0.0D0/
  822. C++   DATA ONE/1.0E0/,ZERO/0.0E0/
  823. C
  824. C        INITIAL SETTINGS
  825. C
  826.       ITER=0
  827.       AONE=ONE+ACU
  828.       AHALF=.5+ACU
  829.       BONE=ONE-ACU
  830.       M1=M-1
  831. C
  832. C     SET UP INITIAL LU DECOMPOSITION
  833. C
  834.       INTL=.TRUE.
  835.       K=1
  836.       CALL UPDATE (K,IFAULT,X,N,M)
  837.       IF (IFAULT.NE.0) RETURN
  838.       INTL=.FALSE.
  839.       DO 10 I=1,M
  840.          K1=IBASE(I)
  841.          RHS(I)=Y(K1)
  842.    10 CONTINUE
  843.       KKK=0
  844.       CALL CALBET (KKK,BETA,RHS,M)
  845. C
  846. C        CALCULATE INITIAL D, TOT AND SIGMA VECTORS
  847. C
  848.       DO 30 J=1,N
  849.          VAL=ZERO
  850.          DO 20 I=1,M
  851.             VAL=VAL+BETA(I)*X(J,I)
  852.    20    CONTINUE
  853.          D(J)=Y(J)-VAL
  854. C++      SIGMA(J)=SIGN(1.,D(J))
  855.          SIGMA(J)=DSIGN(1.D0,D(J))
  856.    30 CONTINUE
  857.       DO 40 J=1,M
  858.          RHS(J)=ZERO
  859.          KKK=IBASE(J)
  860.          SIGMA(KKK)=ZERO
  861.    40 CONTINUE
  862.       DO 60 J=1,N
  863.          DO 50 I=1,M
  864.    50    TOT(I)=TOT(I)-SIGMA(J)*X(J,I)
  865.    60 CONTINUE
  866. C
  867. C        MAIN ITERATIVE LOOP BEGINS
  868. C
  869.    70 CONTINUE
  870. C
  871.       KKK=0
  872.       CALL CALCPI (KKK,PI,TOT,M)
  873. C
  874.       T=AONE
  875.       K=0
  876.       DO 80 J=1,M
  877. C++      IF (ABS(PI(J)).LT.T) GO TO 80
  878.          IF (DABS(PI(J)).LT.T) GO TO 80
  879.          K=J
  880. C++      T=ABS(PI(J))
  881.          T=DABS(PI(J))
  882. C++      RHO=-SIGN(1.,PI(J))
  883.          RHO=-DSIGN(1.D0,PI(J))
  884.    80 CONTINUE
  885.       IF (K.EQ.0) GO TO 220
  886.       KKK=K-1
  887.       RHS(K)=ONE
  888.       CALL CALBET (KKK,BETA,RHS,M)
  889.       RHS(K)=ZERO
  890.       DO 100 I=1,N
  891.          DELTA(I)=ZERO
  892.          IF (SIGMA(I).EQ.ZERO) GO TO 100
  893.          DO 90 J=1,M
  894.             DELTA(I)=DELTA(I)+BETA(J)*X(I,J)
  895.    90    CONTINUE
  896.          DELTA(I)=RHO*DELTA(I)
  897.   100 CONTINUE
  898. C
  899. C        PERFORM PARTIAL SORT OF RATIOS
  900. C
  901.       T=T*.5
  902.       KOUNT=0
  903.       RATIO=BIG
  904.       SUM=AHALF
  905.       SUBT=ZERO
  906.       DO 160 I=1,N
  907.          IF (DELTA(I)*SIGMA(I).LE.ACU) GO TO 160
  908.          TEST=D(I)/DELTA(I)
  909.          IF (TEST.GE.RATIO) GO TO 160
  910. C++      SUM=SUM+ABS(DELTA(I))
  911.          SUM=SUM+DABS(DELTA(I))
  912.          IF (SUM-SUBT.GE.T) GO TO 110
  913. C
  914. C     INSERT I IN LIST
  915. C
  916.          KOUNT=KOUNT+1
  917.          PRICE(KOUNT)=TEST
  918.          INEXT(KOUNT)=I
  919.          GO TO 160
  920. C
  921. C     UPDATE SUM AND KICK IIN OUT OF THE LIST
  922. C
  923.   110    SUM=SUM-SUBT
  924.          RATIO=TEST
  925.          IPT=0
  926.          KKK=0
  927. C
  928. C     IDENTIFY A NEW IIN
  929. C
  930.   120    KKK=KKK+1
  931.          IF (KKK.GT.KOUNT) GO TO 130
  932.          IF (PRICE(KKK).LE.RATIO) GO TO 120
  933.          RATIO=PRICE(KKK)
  934.          IPT=KKK
  935.          GO TO 120
  936.   130    IF (IPT.EQ.0) GO TO 150
  937. C
  938. C     SWITCH VALUES
  939. C
  940.          KKK=INEXT(IPT)
  941. C++      SUBT=ABS(DELTA(KKK))
  942.          SUBT=DABS(DELTA(KKK))
  943.          IF (SUM-SUBT.LT.T) GO TO 140
  944.          PRICE(IPT)=PRICE(KOUNT)
  945.          INEXT(IPT)=INEXT(KOUNT)
  946.          KOUNT=KOUNT-1
  947.          GO TO 110
  948.   140    IIN=INEXT(IPT)
  949.          INEXT(IPT)=I
  950.          PRICE(IPT)=TEST
  951.          GO TO 160
  952.   150    IIN=I
  953. C++      SUBT=ABS(DELTA(I))
  954.          SUBT=DABS(DELTA(I))
  955.   160 CONTINUE
  956. C
  957. C        UPDATE BASIC INDICATORS
  958. C
  959.       IF (KOUNT.EQ.0) GO TO 190
  960.       DO 180 J=1,KOUNT
  961.          KKK=INEXT(J)
  962.          ZZZ=SIGMA(KKK)
  963.          DO 170 L=1,M
  964.             SUBT=ZZZ*X(KKK,L)
  965.             TOT(L)=TOT(L)+SUBT+SUBT
  966.   170    CONTINUE
  967.          SIGMA(KKK)=-SIGMA(KKK)
  968.   180 CONTINUE
  969.   190 CONTINUE
  970.       IOUT=IBASE(K)
  971.       DELTA(IOUT)=RHO
  972.       DO 200 L=1,M
  973.   200 TOT(L)=TOT(L)+RHO*X(IOUT,L)+SIGMA(IIN)*X(IIN,L)
  974.       SIGMA(IOUT)=-RHO
  975.       IBASE(K)=IIN
  976.       CALL UPDATE (K,IFAULT,X,1,M)
  977.       DO 210 J=1,N
  978.   210 D(J)=D(J)-RATIO*DELTA(J)
  979.       SIGMA(IIN)=ZERO
  980.       ITER=ITER+1
  981.       GO TO 70
  982. C
  983. C        CALCULATE OPTIMAL BETA AND SUM OF ABSOLUTE DEVIATIONS
  984. C
  985.   220 DO 230 I=1,M
  986.          K1=IBASE(I)
  987.          RHS(I)=Y(K1)
  988.   230 CONTINUE
  989.       KKK=0
  990.       CALL CALBET (KKK,BETA,RHS,M)
  991.       SAD=ZERO
  992.       DO 240 J=1,N
  993. C++  240 SAD=SAD+ABS(D(J))
  994.  240  SAD=SAD+DABS(D(J))
  995.       RETURN
  996. C
  997.       END
  998.       SUBROUTINE CALBET (KKK,BETA,RHS,M)
  999. C
  1000. C     SUBROUTINE CALBET FOR BACK-SOLVING SYSTEM OF EQUATIONS
  1001. C
  1002.       DIMENSION BETA(20), RHS(20)
  1003. C
  1004.       COMMON /BLOCK/
  1005.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  1006.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  1007.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  1008. C
  1009.       INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2
  1010.       INTEGER KK,KKK,KKK1,LABEL,LEVEL,M,M1
  1011. C
  1012.       DOUBLE PRECISION ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
  1013. C++   REAL             ACU,BETA,BIG,BSAVE,LU,PIE,RHS,SAD,SADT,SAVTOT
  1014.       DOUBLE PRECISION SIG,SIGMA,ZERO,ZLOW,ZSAVE
  1015. C++   REAL             SIG,SIGMA,ZERO,ZLOW,ZSAVE
  1016. C
  1017.       LOGICAL DIRECT,INTL
  1018. C
  1019.       DATA ZERO/0.0D0/
  1020. C++   DATA ZERO/0.0E0/
  1021. C
  1022.       IF (M.GT.1) GO TO 10
  1023.       K=INDEX(1)
  1024.       BETA(K)=RHS(1)/LU(K,1)
  1025.       RETURN
  1026.    10 CONTINUE
  1027.       M1=M-1
  1028.       IF (KKK.EQ.0) GO TO 30
  1029.       DO 20 I=1,KKK
  1030.          K=INDEX(I)
  1031.          BETA(K)=ZERO
  1032.    20 CONTINUE
  1033.    30 KKK=KKK+1
  1034.       K=INDEX(KKK)
  1035.       BETA(K)=RHS(KKK)/LU(K,KKK)
  1036.       IF (KKK.EQ.M) GO TO 60
  1037.       KKK1=KKK+1
  1038.       DO 50 II=KKK1,M
  1039.          K=INDEX(II)
  1040.          BETA(K)=RHS(II)
  1041.          II1=II-1
  1042.          DO 40 I=KKK,II1
  1043.             KK=INDEX(I)
  1044.             BETA(K)=BETA(K)-LU(KK,II)*BETA(KK)
  1045.    40    CONTINUE
  1046.          BETA(K)=BETA(K)/LU(K,II)
  1047.    50 CONTINUE
  1048.    60 CONTINUE
  1049.       DO 70 II=1,M1
  1050.          K1=M-II
  1051.          K=INDEX(K1)
  1052.       DO 70 I=1,II
  1053.          KK=M-I+1
  1054.          K2=INDEX(KK)
  1055.          BETA(K)=BETA(K)-LU(K2,K1)*BETA(K2)
  1056.    70 CONTINUE
  1057.       RETURN
  1058. C
  1059.       END
  1060.       SUBROUTINE UPDATE (KKK,IFAULT,X,N,M)
  1061. C
  1062. C     SUBROUTINE UPDATE FOR UPDATING LU DECOMPOSITION MATRIX
  1063. C
  1064.       COMMON /BLOCK/
  1065.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  1066.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  1067.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  1068. C
  1069.       INTEGER I,IBASE,ICOL,IFAULT,II,II1,IK,IK1,INCOL,INDEX,INPROB
  1070.       INTEGER IPAR,IROW,ISAVE,J,K,KK,KKK,LABEL,LEVEL,LL,M,N
  1071. C
  1072.       DOUBLE PRECISION ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
  1073. C++   REAL             ACU,BIG,BSAVE,LU,PIE,PIVOT,SAD,SADT,SAVTOT
  1074.       DOUBLE PRECISION SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
  1075. C++   REAL             SIG,SIGMA,SUBT,X,ZLOW,ZSAVE
  1076. C
  1077.       LOGICAL DIRECT,INTL
  1078. C
  1079.       DIMENSION X(300,20)
  1080. C
  1081.       IROW=0
  1082.       DO 90 II=KKK,M
  1083.          IF (INTL) GO TO 10
  1084.          IROW=IBASE(II)
  1085.          GO TO 20
  1086.    10    IROW=IROW+1
  1087.          IBASE(II)=IROW
  1088.          IF (IROW.GT.N) IFAULT=1
  1089.          IF (IFAULT.NE.0) RETURN
  1090.    20    DO 30 I=1,M
  1091.             LL=INCOL(I)
  1092.             LU(I,II)=X(IROW,LL)
  1093.    30    CONTINUE
  1094. C
  1095. C     SET UP REPRESENTATION OF INCOMING ROW
  1096. C
  1097.          IF (II.EQ.1) GO TO 60
  1098.          II1=II-1
  1099.          DO 50 ICOL=1,II1
  1100.             K=INDEX(ICOL)
  1101.             SUBT=LU(K,II)
  1102.             J=ICOL+1
  1103.             DO 40 I=J,M
  1104.                K=INDEX(I)
  1105.                LU(K,II)=LU(K,II)-SUBT*LU(K,ICOL)
  1106.    40       CONTINUE
  1107.    50    CONTINUE
  1108. C
  1109. C     FIND MAXIMUM ENTRY
  1110. C
  1111.    60    PIVOT=ACU
  1112.          KK=0
  1113.          DO 70 I=II,M
  1114.             K=INDEX(I)
  1115.             IF (DABS(LU(K,II)).LE.PIVOT) GO TO 70
  1116. C++         IF (ABS(LU(K,II)).LE.PIVOT) GO TO 70
  1117. C++         PIVOT=ABS(LU(K,II))
  1118.             PIVOT=DABS(LU(K,II))
  1119.             KK=I
  1120.    70    CONTINUE
  1121.          IF (KK.EQ.0) GO TO 10
  1122. C
  1123. C     SWITCH ORDER
  1124. C
  1125.          ISAVE=INDEX(KK)
  1126.          INDEX(KK)=INDEX(II)
  1127.          INDEX(II)=ISAVE
  1128. C
  1129. C     PUT IN COLUMNS OF LU ONE AT A TIME
  1130. C
  1131.          IF (II.EQ.M) GO TO 90
  1132.          J=II+1
  1133.          DO 80 I=J,M
  1134.             K=INDEX(I)
  1135.             LU(K,II)=LU(K,II)/LU(ISAVE,II)
  1136.    80    CONTINUE
  1137.    90 CONTINUE
  1138.       KKK=IROW
  1139.       RETURN
  1140. C
  1141.       END
  1142.       SUBROUTINE CALCPI (KKK,PI,TOT,M)
  1143. C
  1144. C     THIS ROUTINE FORWARD SOLVES A LINEAR SYSTEM
  1145. C
  1146.       DIMENSION TOT(20), PI(20)
  1147. C
  1148.       COMMON /BLOCK/
  1149.      -    ACU,BIG,BSAVE(210),LU(20,20),PIE(210),SAD,SADT,SAVTOT(210),
  1150.      -SIG(6000),SIGMA(300),ZLOW,ZSAVE(20),IBASE(20),IK,IK1,INCOL(20),
  1151.      -INDEX(20),INPROB(210),IPAR(20),LABEL(210),LEVEL,DIRECT,INTL
  1152. C
  1153.       INTEGER I,IBASE,II,II1,IK,IK1,INCOL,INDEX,INPROB,IPAR,K,K1,K2,KKK
  1154.       INTEGER KKK1,LABEL,LEVEL,M,M1
  1155. C
  1156.       DOUBLE PRECISION ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
  1157. C++   REAL             ACU,BIG,BSAVE,LU,PI,PIE,SAD,SADT,SAVTOT,SIG
  1158.       DOUBLE PRECISION SIGMA,TOT,ZERO,ZLOW,ZSAVE
  1159. C++   REAL             SIGMA,TOT,ZERO,ZLOW,ZSAVE
  1160. C
  1161.       LOGICAL DIRECT,INTL
  1162. C
  1163. C++   DATA ZERO/0.0E0/
  1164.       DATA ZERO/0.0D0/
  1165.       M1=M-1
  1166. C
  1167.       IF (M.GT.1) GO TO 10
  1168.       K=INDEX(M)
  1169.       PI(M)=TOT(K)/LU(K,M)
  1170.       RETURN
  1171.    10 IF (KKK.EQ.0) GO TO 30
  1172.       DO 20 I=1,KKK
  1173.          PI(I)=ZERO
  1174.    20 CONTINUE
  1175.    30 CONTINUE
  1176.       KKK=KKK+1
  1177.       K=INDEX(KKK)
  1178.       PI(KKK)=TOT(K)
  1179.       IF (KKK.EQ.M) GO TO 60
  1180.       KKK1=KKK+1
  1181.       DO 50 II=KKK1,M
  1182.          K=INDEX(II)
  1183.          PI(II)=TOT(K)
  1184.          II1=II-1
  1185.          DO 40 I=KKK,II1
  1186.    40    PI(II)=PI(II)-LU(K,I)*PI(I)
  1187.    50 CONTINUE
  1188.    60 CONTINUE
  1189.       K=INDEX(M)
  1190.       PI(M)=PI(M)/LU(K,M)
  1191.       DO 80 II=1,M1
  1192.          K1=M-II
  1193.          K=INDEX(K1)
  1194.          DO 70 I=1,II
  1195.             K2=M-I+1
  1196.             PI(K1)=PI(K1)-LU(K,K2)*PI(K2)
  1197.    70    CONTINUE
  1198.          PI(K1)=PI(K1)/LU(K,K1)
  1199.    80 CONTINUE
  1200.       RETURN
  1201.       END
  1202.       DOUBLE PRECISION FUNCTION D1MACH(I)
  1203. C***BEGIN PROLOGUE  D1MACH
  1204. C***DATE WRITTEN   750101   (YYMMDD)
  1205. C***REVISION DATE  820801   (YYMMDD)
  1206. C***CATEGORY NO.  R1
  1207. C***KEYWORDS  MACHINE CONSTANTS
  1208. C***AUTHOR  FOX, P. A., (BELL LABS)
  1209. C           HALL, A. D., (BELL LABS)
  1210. C           SCHRYER, N. L., (BELL LABS)
  1211. C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
  1212. C***DESCRIPTION
  1213. C     D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
  1214. C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
  1215. C     SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
  1216. C     AS FOLLOWS, FOR EXAMPLE
  1217. C
  1218. C          D = D1MACH(I)
  1219. C
  1220. C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
  1221. C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
  1222. C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
  1223. C
  1224. C  DOUBLE-PRECISION MACHINE CONSTANTS
  1225. C  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
  1226. C  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
  1227. C  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
  1228. C  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
  1229. C  D1MACH( 5) = LOG10(B)
  1230. C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
  1231. C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
  1232. C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
  1233. C***ROUTINES CALLED  XERROR
  1234. C***END PROLOGUE  D1MACH
  1235. C
  1236.       INTEGER SMALL(4)
  1237.       INTEGER LARGE(4)
  1238.       INTEGER RIGHT(4)
  1239.       INTEGER DIVER(4)
  1240.       INTEGER LOG10(4)
  1241. C
  1242.       DOUBLE PRECISION DMACH(5)
  1243. C
  1244.       EQUIVALENCE (DMACH(1),SMALL(1))
  1245.       EQUIVALENCE (DMACH(2),LARGE(1))
  1246.       EQUIVALENCE (DMACH(3),RIGHT(1))
  1247.       EQUIVALENCE (DMACH(4),DIVER(1))
  1248.       EQUIVALENCE (DMACH(5),LOG10(1))
  1249. C
  1250. C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
  1251. C
  1252. C     DATA SMALL(1) / ZC00800000 /
  1253. C     DATA SMALL(2) / Z000000000 /
  1254. C
  1255. C     DATA LARGE(1) / ZDFFFFFFFF /
  1256. C     DATA LARGE(2) / ZFFFFFFFFF /
  1257. C
  1258. C     DATA RIGHT(1) / ZCC5800000 /
  1259. C     DATA RIGHT(2) / Z000000000 /
  1260. C
  1261. C     DATA DIVER(1) / ZCC6800000 /
  1262. C     DATA DIVER(2) / Z000000000 /
  1263. C
  1264. C     DATA LOG10(1) / ZD00E730E7 /
  1265. C     DATA LOG10(2) / ZC77800DC0 /
  1266. C
  1267. C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM.
  1268. C
  1269. C     DATA SMALL(1) / O1771000000000000 /
  1270. C     DATA SMALL(2) / O0000000000000000 /
  1271. C
  1272. C     DATA LARGE(1) / O0777777777777777 /
  1273. C     DATA LARGE(2) / O0007777777777777 /
  1274. C
  1275. C     DATA RIGHT(1) / O1461000000000000 /
  1276. C     DATA RIGHT(2) / O0000000000000000 /
  1277. C
  1278. C     DATA DIVER(1) / O1451000000000000 /
  1279. C     DATA DIVER(2) / O0000000000000000 /
  1280. C
  1281. C     DATA LOG10(1) / O1157163034761674 /
  1282. C     DATA LOG10(2) / O0006677466732724 /
  1283. C
  1284. C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS.
  1285. C
  1286. C     DATA SMALL(1) / O1771000000000000 /
  1287. C     DATA SMALL(2) / O7770000000000000 /
  1288. C
  1289. C     DATA LARGE(1) / O0777777777777777 /
  1290. C     DATA LARGE(2) / O7777777777777777 /
  1291. C
  1292. C     DATA RIGHT(1) / O1461000000000000 /
  1293. C     DATA RIGHT(2) / O0000000000000000 /
  1294. C
  1295. C     DATA DIVER(1) / O1451000000000000 /
  1296. C     DATA DIVER(2) / O0000000000000000 /
  1297. C
  1298. C     DATA LOG10(1) / O1157163034761674 /
  1299. C     DATA LOG10(2) / O0006677466732724 /
  1300. C
  1301. C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
  1302. C
  1303. C     DATA SMALL(1) / 00604000000000000000B /
  1304. C     DATA SMALL(2) / 00000000000000000000B /
  1305. C
  1306. C     DATA LARGE(1) / 37767777777777777777B /
  1307. C     DATA LARGE(2) / 37167777777777777777B /
  1308. C
  1309. C     DATA RIGHT(1) / 15604000000000000000B /
  1310. C     DATA RIGHT(2) / 15000000000000000000B /
  1311. C
  1312. C     DATA DIVER(1) / 15614000000000000000B /
  1313. C     DATA DIVER(2) / 15010000000000000000B /
  1314. C
  1315. C     DATA LOG10(1) / 17164642023241175717B /
  1316. C     DATA LOG10(2) / 16367571421742254654B /
  1317. C
  1318. C     MACHINE CONSTANTS FOR THE CRAY 1
  1319. C
  1320. C     DATA SMALL(1) / 200004000000000000000B /
  1321. C     DATA SMALL(2) / 00000000000000000000B /
  1322. C
  1323. C     DATA LARGE(1) / 577777777777777777777B /
  1324. C     DATA LARGE(2) / 000007777777777777777B /
  1325. C
  1326. C     DATA RIGHT(1) / 377214000000000000000B /
  1327. C     DATA RIGHT(2) / 000000000000000000000B /
  1328. C
  1329. C     DATA DIVER(1) / 377224000000000000000B /
  1330. C     DATA DIVER(2) / 000000000000000000000B /
  1331. C
  1332. C     DATA LOG10(1) / 377774642023241175717B /
  1333. C     DATA LOG10(2) / 000007571421742254654B /
  1334. C
  1335. C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  1336. C
  1337. C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
  1338. C     STATIC DMACH(5)
  1339. C
  1340. C     DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
  1341. C     DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
  1342. C     DATA LOG10/40423K,42023K,50237K,74776K/
  1343. C
  1344. C     MACHINE CONSTANTS FOR THE HARRIS 220
  1345. C
  1346. C     DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
  1347. C     DATA LARGE(1),LARGE(2) / [37777777, [37777577 /
  1348. C     DATA RIGHT(1),RIGHT(2) / [20000000, [00000333 /
  1349. C     DATA DIVER(1),DIVER(2) / [20000000, [00000334 /
  1350. C     DATA LOG10(1),LOG10(2) / [23210115, [10237777 /
  1351. C
  1352. C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
  1353. C
  1354. C     DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 /
  1355. C     DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 /
  1356. C     DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 /
  1357. C     DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 /
  1358. C     DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /
  1359. C
  1360. C      MACHINE CONSTANTS FOR THE HP 2100
  1361. C      THREE WORD DOUBLE PRECISION OPTION WITH FTN4
  1362. C
  1363. C      DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
  1364. C      DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
  1365. C      DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
  1366. C      DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
  1367. C      DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
  1368. C
  1369. C
  1370. C      MACHINE CONSTANTS FOR THE HP 2100
  1371. C      FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
  1372. C
  1373. C      DATA SMALL(1), SMALL(2) /  40000B,       0 /
  1374. C      DATA SMALL(3), SMALL(4) /       0,       1 /
  1375. C      DATA LARGE(1), LARGE(2) /  77777B, 177777B /
  1376. C      DATA LARGE(3), LARGE(4) / 177777B, 177776B /
  1377. C      DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
  1378. C      DATA RIGHT(3), RIGHT(4) /       0,    225B /
  1379. C      DATA DIVER(1), DIVER(2) /  40000B,       0 /
  1380. C      DATA DIVER(3), DIVER(4) /       0,    227B /
  1381. C      DATA LOG10(1), LOG10(2) /  46420B,  46502B /
  1382. C      DATA LOG10(3), LOG10(4) /  76747B, 176377B /
  1383. C
  1384. C
  1385. C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  1386. C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
  1387. C     THE PERKIN ELMER (INTERDATA) 7/32.
  1388. C
  1389. C     DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 /
  1390. C     DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
  1391. C     DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 /
  1392. C     DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 /
  1393. C     DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /
  1394. C
  1395. C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
  1396. C
  1397. C     DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 /
  1398. C     DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 /
  1399. C     DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 /
  1400. C     DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 /
  1401. C     DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /
  1402. C
  1403. C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR).
  1404. C
  1405. C     DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 /
  1406. C     DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 /
  1407. C     DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 /
  1408. C     DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 /
  1409. C     DATA LOG10(1),LOG10(2) / "177464202324, "476747767461 /
  1410. C
  1411. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  1412. C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  1413. C
  1414. C     DATA SMALL(1),SMALL(2) /    8388608,           0 /
  1415. C     DATA LARGE(1),LARGE(2) / 2147483647,          -1 /
  1416. C     DATA RIGHT(1),RIGHT(2) /  612368384,           0 /
  1417. C     DATA DIVER(1),DIVER(2) /  620756992,           0 /
  1418. C     DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /
  1419. C
  1420. C     DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 /
  1421. C     DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 /
  1422. C     DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 /
  1423. C     DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 /
  1424. C     DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /
  1425. C
  1426. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  1427. C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  1428. C
  1429. C     DATA SMALL(1),SMALL(2) /    128,      0 /
  1430. C     DATA SMALL(3),SMALL(4) /      0,      0 /
  1431. C
  1432. C     DATA LARGE(1),LARGE(2) /  32767,     -1 /
  1433. C     DATA LARGE(3),LARGE(4) /     -1,     -1 /
  1434. C
  1435. C     DATA RIGHT(1),RIGHT(2) /   9344,      0 /
  1436. C     DATA RIGHT(3),RIGHT(4) /      0,      0 /
  1437. C
  1438. C     DATA DIVER(1),DIVER(2) /   9472,      0 /
  1439. C     DATA DIVER(3),DIVER(4) /      0,      0 /
  1440. C
  1441. C     DATA LOG10(1),LOG10(2) /  16282,   8346 /
  1442. C     DATA LOG10(3),LOG10(4) / -31493, -12296 /
  1443. C
  1444. C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
  1445. C     DATA SMALL(3),SMALL(4) / O000000, O000000 /
  1446. C
  1447. C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
  1448. C     DATA LARGE(3),LARGE(4) / O177777, O177777 /
  1449. C
  1450. C     DATA RIGHT(1),RIGHT(2) / O022200, O000000 /
  1451. C     DATA RIGHT(3),RIGHT(4) / O000000, O000000 /
  1452. C
  1453. C     DATA DIVER(1),DIVER(2) / O022400, O000000 /
  1454. C     DATA DIVER(3),DIVER(4) / O000000, O000000 /
  1455. C
  1456. C     DATA LOG10(1),LOG10(2) / O037632, O020232 /
  1457. C     DATA LOG10(3),LOG10(4) / O102373, O147770 /
  1458. C
  1459. C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
  1460. C
  1461. C     DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 /
  1462. C     DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 /
  1463. C     DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 /
  1464. C     DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 /
  1465. C     DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /
  1466. C
  1467. C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FOR COMPILER
  1468. C
  1469. C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
  1470. C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
  1471. C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
  1472. C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
  1473. C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572/
  1474. C
  1475. C
  1476. C     MACHINE CONSTANTS FOR VAX 11/780
  1477. C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  1478. C    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
  1479. C    *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
  1480. C
  1481. C     DATA SMALL(1), SMALL(2) /        128,           0 /
  1482. C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
  1483. C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
  1484. C     DATA DIVER(1), DIVER(2) /       9472,           0 /
  1485. C     DATA LOG10(1), LOG10(2) /  546979738,  -805665541 /
  1486. C
  1487. C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
  1488. C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
  1489. C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
  1490. C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
  1491. C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFFA84FB /
  1492. C
  1493. C***FIRST EXECUTABLE STATEMENT  D1MACH
  1494. C
  1495.       D1MACH = DMACH(I)
  1496.       RETURN
  1497. C
  1498.       END
  1499.       REAL FUNCTION R1MACH(I)
  1500. C***BEGIN PROLOGUE  R1MACH
  1501. C***DATE WRITTEN   790101   (YYMMDD)
  1502. C***REVISION DATE  820801   (YYMMDD)
  1503. C***CATEGORY NO.  R1
  1504. C***KEYWORDS  MACHINE CONSTANTS
  1505. C***AUTHOR  FOX, P. A., (BELL LABS)
  1506. C           HALL, A. D., (BELL LABS)
  1507. C           SCHRYER, N. L., (BELL LABS)
  1508. C***PURPOSE  RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS
  1509. C***DESCRIPTION
  1510. C     R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
  1511. C     FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
  1512. C     SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
  1513. C     AS FOLLOWS, FOR EXAMPLE
  1514. C
  1515. C          A = R1MACH(I)
  1516. C
  1517. C     WHERE I=1,...,5.  THE (OUTPUT) VALUE OF A ABOVE IS
  1518. C     DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
  1519. C     VARIOUS VALUES OF I ARE DISCUSSED BELOW.
  1520. C
  1521. C  SINGLE-PRECISION MACHINE CONSTANTS
  1522. C  R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
  1523. C  R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
  1524. C  R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
  1525. C  R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
  1526. C  R1MACH(5) = LOG10(B)
  1527. C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
  1528. C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
  1529. C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
  1530. C                 PP. 177-188.
  1531. C***ROUTINES CALLED  XERROR
  1532. C***END PROLOGUE  R1MACH
  1533. C
  1534.       INTEGER SMALL(2)
  1535.       INTEGER LARGE(2)
  1536.       INTEGER RIGHT(2)
  1537.       INTEGER DIVER(2)
  1538.       INTEGER LOG10(2)
  1539. C
  1540.       REAL RMACH(5)
  1541. C
  1542.       EQUIVALENCE (RMACH(1),SMALL(1))
  1543.       EQUIVALENCE (RMACH(2),LARGE(1))
  1544.       EQUIVALENCE (RMACH(3),RIGHT(1))
  1545.       EQUIVALENCE (RMACH(4),DIVER(1))
  1546.       EQUIVALENCE (RMACH(5),LOG10(1))
  1547. C
  1548. C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
  1549. C
  1550. C     DATA RMACH(1) / Z400800000 /
  1551. C     DATA RMACH(2) / Z5FFFFFFFF /
  1552. C     DATA RMACH(3) / Z4E9800000 /
  1553. C     DATA RMACH(4) / Z4EA800000 /
  1554. C     DATA RMACH(5) / Z500E730E8 /
  1555. C
  1556. C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS.
  1557. C
  1558. C     DATA RMACH(1) / O1771000000000000 /
  1559. C     DATA RMACH(2) / O0777777777777777 /
  1560. C     DATA RMACH(3) / O1311000000000000 /
  1561. C     DATA RMACH(4) / O1301000000000000 /
  1562. C     DATA RMACH(5) / O1157163034761675 /
  1563. C
  1564. C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
  1565. C
  1566. C     DATA RMACH(1) / 00014000000000000000B /
  1567. C     DATA RMACH(2) / 37767777777777777777B /
  1568. C     DATA RMACH(3) / 16404000000000000000B /
  1569. C     DATA RMACH(4) / 16414000000000000000B /
  1570. C     DATA RMACH(5) / 17164642023241175720B /
  1571. C
  1572. C     MACHINE CONSTANTS FOR THE CRAY 1
  1573. C
  1574. C     DATA RMACH(1) / 200004000000000000000B /
  1575. C     DATA RMACH(2) / 577777777777777777777B /
  1576. C     DATA RMACH(3) / 377214000000000000000B /
  1577. C     DATA RMACH(4) / 377224000000000000000B /
  1578. C     DATA RMACH(5) / 377774642023241175720B /
  1579. C
  1580. C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  1581. C
  1582. C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
  1583. C     STATIC RMACH(5)
  1584. C
  1585. C     DATA SMALL/20K,0/,LARGE/77777K,177777K/
  1586. C     DATA RIGHT/35420K,0/,DIVER/36020K,0/
  1587. C     DATA LOG10/40423K,42023K/
  1588. C
  1589. C     MACHINE CONSTANTS FOR THE HARRIS 220
  1590. C
  1591. C     DATA SMALL(1),SMALL(2) / [20000000, [00000201 /
  1592. C     DATA LARGE(1),LARGE(2) / [37777777, [00000177 /
  1593. C     DATA RIGHT(1),RIGHT(2) / [20000000, [00000352 /
  1594. C     DATA DIVER(1),DIVER(2) / [20000000, [00000353 /
  1595. C     DATA LOG10(1),LOG10(2) / [23210115, [00000377 /
  1596. C
  1597. C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES.
  1598. C
  1599. C     DATA RMACH(1) / O402400000000 /
  1600. C     DATA RMACH(2) / O376777777777 /
  1601. C     DATA RMACH(3) / O714400000000 /
  1602. C     DATA RMACH(4) / O716400000000 /
  1603. C     DATA RMACH(5) / O776464202324 /
  1604. C
  1605. C     MACHINE CONSTANTS FOR THE HP 2100
  1606. C
  1607. C     3 WORD DOUBLE PRECISION WITH FTN4
  1608. C
  1609. C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
  1610. C     DATA LARGE(1), LARGE(2) / 77777B, 177776B /
  1611. C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
  1612. C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
  1613. C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
  1614. C
  1615. C     MACHINE CONSTANTS FOR THE HP 2100
  1616. C     4 WORD DOUBLE PRECISION WITH FTN4
  1617. C
  1618. C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
  1619. C     DATA LARGE91), LARGE(2) / 77777B, 177776B /
  1620. C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
  1621. C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
  1622. C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
  1623. C
  1624. C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  1625. C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86  AND
  1626. C     THE PERKIN ELMER (INTERDATA) 7/32.
  1627. C
  1628. C     DATA RMACH(1) / Z00100000 /
  1629. C     DATA RMACH(2) / Z7FFFFFFF /
  1630. C     DATA RMACH(3) / Z3B100000 /
  1631. C     DATA RMACH(4) / Z3C100000 /
  1632. C     DATA RMACH(5) / Z41134413 /
  1633. C
  1634. C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR).
  1635. C
  1636. C     DATA RMACH(1) / "000400000000 /
  1637. C     DATA RMACH(2) / "377777777777 /
  1638. C     DATA RMACH(3) / "146400000000 /
  1639. C     DATA RMACH(4) / "147400000000 /
  1640. C     DATA RMACH(5) / "177464202324 /
  1641. C
  1642. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  1643. C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  1644. C
  1645. C     DATA SMALL(1) /    8388608 /
  1646. C     DATA LARGE(1) / 2147483647 /
  1647. C     DATA RIGHT(1) /  880803840 /
  1648. C     DATA DIVER(1) /  889192448 /
  1649. C     DATA LOG10(1) / 1067065499 /
  1650. C
  1651. C     DATA RMACH(1) / O00040000000 /
  1652. C     DATA RMACH(2) / O17777777777 /
  1653. C     DATA RMACH(3) / O06440000000 /
  1654. C     DATA RMACH(4) / O06500000000 /
  1655. C     DATA RMACH(5) / O07746420233 /
  1656. C
  1657. C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  1658. C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
  1659. C
  1660. C     DATA SMALL(1),SMALL(2) /   128,     0 /
  1661. C     DATA LARGE(1),LARGE(2) / 32767,    -1 /
  1662. C     DATA RIGHT(1),RIGHT(2) / 13440,     0 /
  1663. C     DATA DIVER(1),DIVER(2) / 13568,     0 /
  1664. C     DATA LOG10(1),LOG10(2) / 16282,  8347 /
  1665. C
  1666. C     DATA SMALL(1),SMALL(2) / O000200, O000000 /
  1667. C     DATA LARGE(1),LARGE(2) / O077777, O177777 /
  1668. C     DATA RIGHT(1),RIGHT(2) / O032200, O000000 /
  1669. C     DATA DIVER(1),DIVER(2) / O032400, O000000 /
  1670. C     DATA LOG10(1),LOG10(2) / O037632, O020233 /
  1671. C
  1672. C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES.
  1673. C
  1674. C     DATA RMACH(1) / O000400000000 /
  1675. C     DATA RMACH(2) / O377777777777 /
  1676. C     DATA RMACH(3) / O146400000000 /
  1677. C     DATA RMACH(4) / O147400000000 /
  1678. C     DATA RMACH(5) / O177464202324 /
  1679. C
  1680. C     MACHINE CONSTANTS FOR THE VAX 11/780
  1681. C    (EXPRESSED IN INTEGER AND HEXADECIMAL)
  1682. C  ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS***
  1683. C  *** THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
  1684. C
  1685. C     DATA SMALL(1) /       128 /
  1686. C     DATA LARGE(1) /    -32769 /
  1687. C     DATA RIGHT(1) /     13440 /
  1688. C     DATA DIVER(1) /     13568 /
  1689. C     DATA LOG10(1) / 547045274 /
  1690. C
  1691. C     DATA SMALL(1) / Z00000080 /
  1692. C     DATA LARGE(1) / ZFFFF7FFF /
  1693. C     DATA RIGHT(1) / Z00003480 /
  1694. C     DATA DIVER(1) / Z00003500 /
  1695. C     DATA LOG10(1) / Z209B3F9A /
  1696. C
  1697. C***FIRST EXECUTABLE STATEMENT  R1MACH
  1698. C
  1699.       R1MACH = RMACH(I)
  1700.       RETURN
  1701. C
  1702.       END
  1703.     3   DL1BL3= BLISS VOL II P.294... ALL DATA
  1704.    60 1  5.0
  1705.   (1X,F6.3,F3.0,2F6.3)
  1706.     125  1   300   402
  1707.     164  1   156   258
  1708.    -065  1   401   503
  1709.     120  1   137   239
  1710.     034  1   214   316
  1711.     010  1   084   186
  1712.    -079  1   137   239
  1713.    -533  1   031   133  X
  1714.    -611  1   227   329  X
  1715.     076  1   092   194
  1716.     319  1   427   228
  1717.     465  1   487   288
  1718.     313  1   487   288
  1719.     193  1   487   288
  1720.     129  1   589   390
  1721.     209  1   571   372
  1722.     284  1   446   247
  1723.     173  1   570   371
  1724.     264  1   472   273
  1725.     132  1   476   277
  1726.     324  1   696   321
  1727.     367  1   729   354
  1728.     321  1   509   134
  1729.     375  1   559   184
  1730.     345  1   679   304
  1731.     341  1   583   208
  1732.     327  1   742   367
  1733.     256  1   781   406
  1734.     256  1   739   364
  1735.     214  1   865   490
  1736.     501  1   723   223
  1737.     318  1   940   440
  1738.     317  1   903   403
  1739.     349  1   910   410
  1740.     402  1   684   184
  1741.     374  1   904   404
  1742.     340  1   887   387
  1743.     598  1   593   093
  1744.     444  1   640   140
  1745.     543  1   512   012
  1746.     559  1   832   156
  1747.     705  1   817   141
  1748.     585  1   881   205
  1749.     588  1   848   172
  1750.     532  1   857   181
  1751.     593  1   808   132
  1752.     593  1   829   153
  1753.     559  1   872   196
  1754.     611  1   805   129
  1755.     561  1   904   228
  1756.     561  1  1047   246
  1757.     733  1   942   141
  1758.     553  1  1194   393
  1759.     593  1  1090   289
  1760.     525  1  1000   199
  1761.     580  1  1054   253
  1762.     711  1   956   155
  1763.     546  1  1101   300
  1764.     632  1  1060   259
  1765.     733  1   893   092
  1766.  
  1767.  
  1768.  
  1769.  
  1770.  
  1771.  
  1772.  
  1773.  
  1774.  
  1775.  
  1776.  
  1777.  
  1778.  
  1779.  
  1780.  
  1781.  
  1782.  
  1783.     3   DL1BL3= BLISS VOL II P.294... ALL DATA
  1784.         PROBLEM TITLE   DL1BL3= BLISS VOL II P.294... ALL DATA
  1785.     NUMBER OF OBSERVATION=   60
  1786.     NUMBER OF PARAMETERS=    3
  1787.     MINIMUM NUMBER OF PARAMETERS CONSIDERED=    1
  1788.     PERCENTAGE DEVIATION FROM OPTIMALITY ALLOWED=  5.00
  1789.         **BEST SUBSET LAV PROGRAM
  1790.  IFAULT=  0
  1791.  
  1792.  
  1793.    BEST RESULTS FOR LAV SUBSET OF SIZE=  1
  1794.     **************SUM OF ABSOLUTE VALUES=          7.659
  1795.  
  1796.  BETA(  2)           .550
  1797.  
  1798.  
  1799.    BEST RESULTS FOR LAV SUBSET OF SIZE=  2
  1800.     **************SUM OF ABSOLUTE VALUES=          5.407
  1801.  
  1802.  BETA(  2)           .818
  1803.  
  1804.  BETA(  3)          -.782
  1805.  
  1806.  
  1807.    BEST RESULTS FOR LAV SUBSET OF SIZE=  3
  1808.     **************SUM OF ABSOLUTE VALUES=          3.877
  1809.  
  1810.  BETA(  3)         -1.116
  1811.  
  1812.  BETA(  2)           .628
  1813.  
  1814.  BETA(  1)           .235
  1815.  
  1816.  
  1817.      ITERATION COUNT=    106
  1818.     0
  1819.