C ALGORITHM 458 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 16, NO. 10, C P. 629. PROGRAM MAIN(INPUT,OUTPUT,PUNCH,TAPE1=INPUT, MAIN0001 1 TAPE2=PUNCH,TAPE3=OUTPUT) C C C PROGRAM MAIN IS THE CALLING PROGRAM FOR SUBROUTINE C APPROX. THE PROBLEM TO BE SOLVED IS C C MINIMIZE Z = ABS(E(1)) + ... + ABS(E(N)) C SUBJECT TO C FX + E = T C X AND E ARE VECTORS OF VARIABLES HAVING DIMENSION M AND N C RESPECTIVELY. C C THREE USER OPTIONS ARE PROVIDED WHICH CAN BE USED C IN ANY COMBINATION. USE OF THESE OPTIONS IS CONTROLLED C BY THE PARAMETERS ON THE PROBLEM ID CARDS (SEE DESCRIPTION C BELOW). C C PUNCH OPTION - IF THERE IS A POSSIBILITY THAT THE PROBLEM C WILL BE RESOLVED IN A LATER RUN ( TO FIND SOME HIGHER C ORDER APPROXIMATION) THE PUNCH OPTION SHOULD BE USED. C WHEN USED, THE FINAL VALUES FOR INBASE, AINV, AND Y ARE C PUNCHED ON CARDS IN A FORMAT WHICH CAN BE READ AT THE C BEGINNING OF A LATER RUN. THE FIRST PART OF THE PROBLEM ID C CARD IS ALSO PUNCHED. C C ADVANCED START OPTION - THIS OPTION CAN BE USED ONLY IF C VALUES FOR INBASE, AINV, AND Y ARE AVAILABLE FROM A PRE- C VIOUS RUN. WHEN THIS OPTION IS USED, THE PROGRAM BEGINS C BY READING THESE VALUES IN THE SAME FORMAT AS THEY WERE C PUNCHED IN THE PREVIOUS RUN. C C STEP-WISE OPTION - THIS OPTION PROVIDES THE CAPABILITY C OF INCREASING THE ORDER OF THE APPROXIMATION ITERATIVELY C UNTIL EITHER THE PROGRAM RUNS OUT OF DATA OR A DESIRED C APPROXIMATION ACCURACY IS REACHED. IF THIS OPTION IS C USED, TWO PARAMETERS CONTROL THE OPERATION OF THE C ALGORITHM. INCRM IS THE NUMBER OF COLUMNS WHICH WILL BE C ADDED TO F AT EACH PASS. TEST IS THE DESIRED APPROX- C IMATION ACCURACY. THE PROGRAM WILL KEEP STEPPING UP THE C ORDER OF THE APPROXIMATION BY INCRM UNTIL EITHER THE DATA C IS EXHAUSTED OR ZOPT .LE. TEST. C C C THE FIRST CARD READ BY MAIN SPECIFIES HOW MANY C PROBLEMS ARE TO BE SOLVED DURING THIS RUN (FORMAT I3). C C FOR EACH PROBLEM TO BE SOLVED, THE FOLLOWING DATA IS C REQUIRED- C C PROBLEM ID CARD - CONTAINS PROBLEM IDENTIFICATION AND C PARAMETERS AS DESCRIBED BELOW (FORMAT A5,5I5,F10.0). C INBASE - (FOR ADVANCED START OPTION ONLY),(FORMAT 16I5). C Y - (FOR ADVANCED START OPTION ONLY),(FORMAT 5E15.7). C T - (FORMAT 8F10.0). C ROWS OF FT - (FORMAT 8F10.0). C ROWS TO BE ADDED TO FT - (FOR STEP-WISE OPTION ONLY), C (FORMAT 8F10.0). C END OF PROBLEM DATA CARD - (999999999. IN COLUMNS 1-10). C C C C PROGRAM MAIN MUST CONTAIN THE FOLLOWING DIMENSION C STATEMENT- C DIMENSION T(N),FT(M,N),INBASE(N),AINV(N,N),Y(N),XOPT(N), C TEMP(8) C WHERE M AND N ARE LARGE ENOUGH TO SOLVE THE LARGEST C PROBLEM. DIMENSION T(50),FT(15,50),INBASE(50),AINV(50,50),Y(50) 1,XOPT(50),TEMP(8) INTEGER PUNCH END = 999999999. READ (1,10) NP C NP = THE NUMBER OF PROBLEMS TO BE SOLVED. 10 FORMAT (I3) DO 350 INP = 1,NP C READ PROBLEM IDENTIFICATION CARD. READ (1,20) ID,MD,M,N,PUNCH,INCRM,TEST 20 FORMAT (A5,5I5,F10.0) C ID = 5 CHARACTER ALPHANUMERIC PROBLEM IDENTIFICATION. C MD = MODE OF OPERATION INDICATOR. C IF MD = 1, PROBLEM IS TO BE SOLVED FROM THE START. C IF MD = 2, PROBLEM IS TO BE SOLVED FROM AN ADVANCED C START. (DATA FOR INBASE, AINV, AND Y AT BEST KNOWN C SOLUTION ARE REQUIRED BY PROGRAM WHEN MD = 2.) C M = INITIAL ORDER OF APPROXIMATION. C N = NUMBER OF DATA POINTS TO BE USED. C PUNCH = PUNCHING INDICATOR. C IF PUNCH = 0, INFORMATION AT OPTIMAL SOLUTION IS C NOT SAVED FOR FUTURE ADVANCED START. OTHERWISE, C INFORMATION IS PUNCHED. C INCRM = PARAMETER WHICH INDICATES HOW MUCH TO STEP UP C ORDER OF APPROXIMATION EACH RUN IF SUBROUTINE IS TO C BE USED STEP-WISE. IF INCRM = 0, PROGRAM WILL CALL C APPROX ONLY ONCE AND ZOPT HAS NO MEANING. C TEST = DESIRED ACCURACY OF APPROXIMATION. IF INCRM .NE. 0 C PROGRAM WILL KEEP INCREASING ORDER OF APPROXIMATION C UNTIL EITHER DATA IS EXHAUSTED OR ZOPT .LE. TEST. TEMP(1)=0. IF (MD.EQ.1) GO TO 50 C PROBLEM TO USE ADVANCED START OPTION. C READ INBASE. READ (1,30) (INBASE(I),I=1,N) 30 FORMAT (16I5) C READ AINV. READ (1,40) ((AINV(I,J),J=1,N),I=1,N) 40 FORMAT (5E15.7) C READ Y. READ (1,40) (Y(I),I=1,N) C PROBLEM TO START FROM BEGINNING. C READ T. 50 READ (1,60) (T(I),I=1,N) 60 FORMAT (8F10.0) C READ M ROWS OF FT. DO 65 I=1,M READ (1,60) (FT(I,J),J=1,N) 65 CONTINUE C WRITE PROBLEM DATA ON OUTPUT LISTING. WRITE (3,70) ID 70 FORMAT (1H1,//10X,9HPROBLEM ,A5) WRITE (3,80) 80 FORMAT (///20X, 8HT VECTOR) WRITE (3,90) (T(I),I=1,N) 90 FORMAT (1H0,8(2X,F10.4)) WRITE (3,100) 100 FORMAT (///20X,18HMATRIX F TRANSPOSE) DO 110 I=1,M WRITE (3,90) (FT(I,J),J=1,N) 110 CONTINUE 120 IER=0 C INBASE, AINV, AND Y ARE REQUIRED AT CALL ONLY IF MD = 2. CALL APPROX(MD,M,N,T,FT,INBASE,AINV,Y,XOPT,ZOPT,IER) C TEST FOR ABNORMAL RETURN. IF (IER-1) 170,150,130 130 WRITE (3,140) 140 FORMAT(1H0,10X,33H**ERROR**ITERATION LIMIT EXCEEDED) GO TO 300 150 WRITE(3,160) 160 FORMAT (1H0,10X,24H**ERROR**SINGULAR MATRIX, 1 27H GENERATED-CHECK INPUT DATA) GO TO 300 C WRITE OPTIMAL SOLUTION. 170 WRITE (3,180) 180 FORMAT (1H0,10X,16HOPTIMAL SOLUTION) DO 190 I=1,M WRITE (3,192) I,XOPT(I) 190 CONTINUE 192 FORMAT (1H0,2HX(,I2,4H) = ,E15.7) WRITE (3,194) ZOPT 194 FORMAT (1H0,10HMINIMUM = ,E15.7) C TEST FOR STEP-WISE OPERATION. IF (INCRM.EQ.0) GO TO 300 C STEP-WISE OPERATION INDICATED. IF (ZOPT .LE. TEST) GO TO 300 C INCREASE ORDER OF APPROXIMATION. WRITE (3,200) 200 FORMAT (1H0,10X,30HACCURACY BELOW SPECIFIED LIMIT) IF (M+INCRM .LE. 100) GO TO 220 WRITE (3,210) 210 FORMAT (1H0,10X,26HCANNOT ADD INCRM MORE ROWS) GO TO 300 220 MP1 = M+1 C READ INCRM ADDITIONAL ROWS TO BE ADDED TO FT. M=M+INCRM DO 260 I=MP1,M J=0 225 READ (1,60) (TEMP(L),L=1,8) C TEST FOR END OF DATA. IF (TEMP(1) .NE. END) GO TO 240 WRITE (3,230) 230 FORMAT (1H0,10X,11HOUT OF DATA) GO TO 300 240 DO 250 L=1,8 J=J+1 IF (J .GT. N) GO TO 260 250 FT(I,J) = TEMP(L) GO TO 225 260 CONTINUE WRITE (3,270) INCRM 270 FORMAT (1H0,10X,22HORDER OF APPROXIMATION, 1 13H INCREASED BY,2X,I5) WRITE (3,280) 280 FORMAT (1H0,10X,23HFOLLOWING ROWS ADDED TO, 1 19H MATRIX F TRANSPOSE) DO 290 I=MP1,M WRITE (3,90) (FT(I,J),J=1,N) 290 CONTINUE MD = 2 GO TO 120 C PROBLEM COMPLETED. 300 IF (PUNCH .EQ. 0) GO TO 310 C PUNCH OPTION INDICATED. C PUNCH DATA FOR FUTURE ADVANCED START. MD = 2 WRITE (2,20) ID,MD,M,N WRITE (2,30) (INBASE(I),I=1,N) WRITE (2,40) ((AINV(I,J),J=1,N),I=1,N) WRITE (2,40) (Y(I),I=1,N) 310 IF (TEMP(1) .EQ. END) GO TO 350 320 READ (1,60) (TEMP(1)) IF (TEMP(1) .NE. END) GO TO 320 350 CONTINUE STOP END SUBROUTINE APPROX (MD,M,N,T,FT,INBASE,AINV,Y,XOPT, APRX0001 1 ZOPT,IER) C C C THIS SUBROUTINE SOLVES THE DISCRETE LINEAR L1 C APPROXIMATION PROBLEM USING THE SUBOPTIMIZATION METHOD OF C INTERVAL LINEAR PROGRAMMING. THE PROBLEM TO BE SOLVED IS C C MINIMIZE Z = ABS(E(1)) + ... + ABS(E(N)) C SUBJECT TO C FX + E = T C C WHERE F IS A GIVEN N BY M MATRIX, T IS A GIVEN N VECTOR, C X AND E ARE VECTORS OF VARIABLES HAVING DIMENSION M AND N C RESPECTIVELY. C C SUBROUTINE APPROX IS DESIGNED TO BE USED IN TWO MODES- C (1) TO SOLVE A PROBLEM FROM SCRATCH, AND C (2) TO SOLVE A PROBLEM USING AN ADVANCED START FROM A C PREVIOUS RUN. C THE ADVANCED START MODE IS USEFUL IF THE OPTIMAL VALUE OF C Z IS TOO LARGE ON A GIVEN PROBLEM (I.E. THE APPROXIMATION C IS TOO POOR) AND THE PROBLEM IS TO BE RERUN AFTER ADDING C ADDITIONAL COLUMNS TO THE F MATRIX (I.E. INCREASING C THE ORDER OF THE APPROXIMATION). C C SUBROUTINE APPROX IS COMPLETELY SELF-CONTAINED AND C COMMUNICATION IS ACHIEVED SOLELY THROUGH THE ARGUMENT C LIST. THE MEANING OF THE PARAMETERS ARE AS FOLLOWS- C MD = THE MODE OF OPERATION INDICATOR. C IF MD = 1, THE PROBLEM IS SOLVED FROM THE BEGINNING. C IF MD = 2, THE PROBLEM IS TO BE SOLVED USING AN C ADVANCED START FROM A PREVIOUS RUN. C M = THE NUMBER OF COLUMNS IN THE F MATRIX (IF MD = 2, C M MUST BE THE MODIFIED VALUE.) C N = THE NUMBER OF ROWS IN MATRIX F. C T = THE RIGHT HAND SIDE VECTOR FOR THE PROBLEM. C FT = THE TRANSPOSE OF MATRIX F. C INBASE = A VECTOR WHICH CONTAINS INDICES OF BASIC COLUMNS C IN THE OPTIMAL SOLUTION TO THE LINEAR PROGRAM WHEN C APPROX RETURNS CONTROL. C AINV = A MATRIX WHICH CONTAINS THE INVERSE OF THE MATRIX C OF BASIC COLUMNS WHEN APPROX RETURNS CONTROL. C Y = A VECTOR CONTAINING THE OPTIMAL DUAL SOLUTION WHEN C APPROX RETURNS CONTROL. C XOPT = A VECTOR CONTAINING THE OPTIMAL X-VALUES WHEN C APPROX RETURNS CONTROL. C ZOPT = THE OPTIMAL VALUE OF THE OBJECTIVE FUNCTION C WHEN APPROX RETURNS CONTROL. C IER = ERROR INDICATOR WHEN APPROX RETURNS CONTROL. IER=0 C INDICATES NORMAL RETURN. C C NO INITIAL VALUES ARE REQUIRED FOR INBASE, AINV, AND C Y WHEN APPROX IS CALLED WITH MD = 1. WHEN MD = 2, AN C ADVANCED START IS INDICATED. THESE VARIABLES MUST THEN C CONTAIN THEIR FINAL VALUES FROM THE PREVIOUS RUN. THE C USER WILL THUS WANT TO MAKE PROVISIOUS FOR SAVING THESE C VALUES IN THE CALLING PROGRAM SO THAT THEY CAN BE REUSED C IF NEEDED. C C C THE CALLING PROGRAM AND APPROX SHOULD CONTAIN THE C FOLLOWING DIMENSION STATEMENT- C DIMENSION T(N),FT(M,N),INBASE(N),AINV(N,N),Y(N),XOPT(M) C C APPROX MUST ALSO CONTAIN THE FOLLOWING DIMENSION C STATEMENT- C DIMENSION BP(N),BM(N),AR(N),ARAINV(N),Q(N),GAMMA(N),DEL(N) C ,TEMP(N) DIMENSION T(50),FT(15,50),INBASE(50),AINV(50,50),Y(50) 1,XOPT(15) DIMENSION BP(50),BM(50),AR(50),ARAINV(50),Q(50), = GAMMA(50),DEL(50),TEMP(50) INTEGER ENT,QQ,ADBASE,Q,P EQUIVALENCE (GAMMA,DEL) C EPSI IS THE SINGULAR MATRIX ERROR MESSAGE CRITERION. C THE VALUE OF EPSI CAN BE REDUCED FOR ILL CONDITIONED C PROBLEMS. EPSI = .0000001 IF (MD.EQ.2) GO TO 70 C PROBLEM TO BE SOLVED FROM THE BEGINNING. C DEFINE INITIAL SUBPROBLEM. IT = 1 ADBASE = N+1 DO 20 I=1,N BP(I) = 1.0 BM(I) = -1.0 AR(I) = FT(1,I) INBASE(I) = I C INITIALIZE AINV AS THE IDENTITY MATRIX. DO 10 J=1,N 10 AINV(I,J) = 0.0 20 AINV(I,I) = 1.0 C FIND THE INITIAL Y VECTOR. DO 60 I=1,N IF(T(I)) 30,40,50 30 Y(I) = -1.0 GO TO 60 40 Y(I) = 0. GO TO 60 50 Y(I) = 1.0 60 CONTINUE GO TO 100 C PROBLEM TO BE SOLVED FROM AN ADVANCED START. 70 ADBASE = M+N DO 90 I=1,N IF (INBASE(I).LE.N) GO TO 80 BP(I) = 0.0 BM(I) = 0.0 GO TO 90 80 BP(I) = 1.0 BM(I) =-1.0 90 AR(I) = FT(M,I) IT = IT+1 100 BRM = 0.0 BRP = 0.0 C BEGIN GENERAL ITERATION. C DETERMINE DEL (THE AMOUNT OF INFEASIBILITY IN THE BOTTOM C CONSTRAINT OF THE CURRENT SUBPROBLEM). 110 CONTINUE S = 0.0 DO 120 I=1,N 120 S = S + AR(I)*Y(I) D = S-BRP IF (D.GT.0.) GO TO 130 D = S - BRM IF (D.GE.0.) GO TO 430 130 CONTINUE DO 140 I=1,N ARAINV(I) = 0.0 DO 140 J=1,N 140 ARAINV(I) = ARAINV(I)+AR(J)*AINV(J,I) C CALCULATE GAMMA VECTOR (THE VECTOR OF MARGINAL COSTS FOR C MOVING TOWARD FEASIBILITY). DO 170 I=1,N TEMP(I) = 0.0 DO 150 J=1,N 150 TEMP(I) = TEMP(I)+T(J)*AINV(J,I) IF (ARAINV(I).NE.0.0)GO TO 160 GAMMA(I) = -1.0 GO TO 170 160 GAMMA(I) = TEMP(I)/ARAINV(I) IF (D.LT.0.0) GAMMA(I)=-GAMMA(I) 170 CONTINUE C FIND Q VECTOR (THE VECTOR OF INDICES WHICH INDICATE THE C VARIABLES WHICH CAN BE CHANGED TO MOVE TOWARD FEASIBILITY). QQ=0 DO 210 L=1,N DO 180 I=1,N IF (GAMMA(I).LT.0.0) GO TO 180 S = GAMMA(I) J=I GO TO 190 180 CONTINUE GO TO 215 190 DO 200 I=1,N IF (GAMMA(I).LT.0. .OR. GAMMA(I).GE.S)GO TO 200 S = GAMMA(I) J=I 200 CONTINUE QQ=QQ+1 Q(L) = J 210 GAMMA(J) = -1.0 C CALCULATE DELTA VECTOR (THE VECTOR INDICATING THE MAXIMUM C PERMISSABLE CHANGES IN THE VARIABLES). 215 DO 260 I=1,QQ K=Q(I) S = 0.0 II = INBASE(K) IF (II.LE.N) GO TO 230 L = II-N DO 220 J=1,N 220 S = S + FT(L,J)*Y(J) GO TO 240 230 S = Y(II) 240 IF(D*ARAINV(K).LE.0.) GO TO 250 DEL(K) = BM(K)-S GO TO 260 250 DEL(K) = BP(K)-S 260 CONTINUE C DETERMINE P (THE NUMBER OF VARIABLES CHANGED C THIS ITERATION). DO 280 I=1,QQ P = I S = 0.0 DO 270 J=1,I K = Q(J) 270 S = S+DEL(K)*ARAINV(K) IF (ABS(S).GE.ABS(D)) GO TO 290 280 CONTINUE C CALCULATE THETA (THE AMOUNT WHICH THE PTH VECTOR IS C CHANGED). 290 L = P-1 S=0.0 IF (L .LT. 1) GO TO 310 DO 300 J=1,L K = Q(J) 300 S = S+DEL(K)*ARAINV(K) 310 K=Q(P) THETA = -(D+S)/ARAINV(K) C UPDATE Y VECTOR (THE OPTIMAL SOLUTION TO THE CURRENT C SUBPROBLEM). DO 320 I=1,N 320 TEMP(I) = 0.0 IF (L .LT. 1) GO TO 340 DO 330 I=1,L K = Q(I) 330 TEMP(K) = DEL(K) 340 K = Q(P) TEMP(K) = THETA DO 360 I=1,N S = 0.0 DO 350 J=1,N 350 S = S+AINV(I,J)*TEMP(J) 360 Y(I) = Y(I)+S K = Q(P) BP(K) = BRP BM(K) = BRM INBASE(K) = ADBASE C CALCULATE NEW AINV MATRIX. DO 370 I=1,N TEMP(I) = 0.0 DO 370 L=1,N 370 TEMP(I) = TEMP(I)+AR(L)*AINV(L,I) IF (ABS(TEMP(K)).GT. EPSI ) GO TO 380 C SINGULAR MATRIX INDICATED. SET ERROR TAG AND TERMINATE. IER = 1 RETURN 380 DO 390 I=1,N 390 AINV(I,K) = AINV(I,K)/TEMP(K) DO 420 J=1,N IF(J.EQ.K) GO TO 410 DO 400 I=1,N AINV(I,J) = AINV(I,J)-AINV(I,K)*TEMP(J) 400 CONTINUE 410 CONTINUE 420 CONTINUE C FIND S (THE LARGEST INFEASIBILITY), AND ENT(THE INDEX OF C THE CORRESPONDING CONSTRAINT). 430 TEMP(1) = 0. L = M+N DO 510 I=1,L DO 440 J=1,N IF(INBASE(J).EQ.I) GO TO 510 440 CONTINUE IF(I.LE.N) GO TO 470 S = 0. II = I-N DO 450 J=1,N 450 S = S + FT(II,J)*Y(J) IF(S.EQ.0.) GO TO 510 460 IF(ABS(S) .LE. TEMP(1)) GO TO 510 TEMP(1) = ABS(S) ENT = I GO TO 510 470 S = Y(I) IF (S-1.) 490,490,480 480 S = S-1. GO TO 460 490 IF (S+1.) 500,510,510 500 S = S+1. GO TO 460 510 CONTINUE S = TEMP(1) IF (S.EQ.0.) GO TO 560 C PRESENT SOLUTION INFEASIBLE. START THE NEXT ITERATION. IT = IT+1 C DEFINE THE NEXT SUBPROBLEM. IF (ENT.LE.N) GO TO 530 BRP = 0.0 BRM = 0.0 L = ENT-N DO 520 J=1,N 520 AR(J) = FT(L,J) GO TO 550 530 BRP = 1.0 BRM = -1.0 DO 540 J=1,N 540 AR(J) = 0.0 AR(ENT) = 1.0 550 ADBASE = ENT IF (IT.LE.10*(M+1)) GO TO 110 C ITERATION LIMIT EXCEEDED. SET ERROR TAG AND TERMINATE. IER = 2 RETURN C OPTIMAL DUAL SOLUTION FOUND. CALCULATE PRIMAL SOLUTION. 560 DO 600 J=1,M L = J+N DO 570 I=1,N IF(INBASE(I).EQ.L) GO TO 580 570 CONTINUE XOPT(J) = 0. GO TO 600 580 TEMP(I) = 0. DO 590 L=1,N 590 TEMP(I) = TEMP(I)+T(L)*AINV(L,I) XOPT(J) = TEMP(I) 600 CONTINUE ZOPT = 0. DO 610 I=1,N 610 ZOPT = ZOPT + Y(I)*T(I) RETURN END