C     THIS PROGRAM IS A DRIVER FOR THE SUBROUTINE  L1
C WHICH SOLVES AN OVERDETERMINED SYSTEM OF LINEAR EQUATIONS
C IN THE L1 NORM, USING A DUAL SIMPLEX METHOD.
C     THE OVERDETERMINED SYSTEM HAS THE FORM  CA=F
C   C IS A GIVEN REAL N BY M MATRIX OF RANK K, K.LE.M.LE.N
C   F IS A GIVEN REAL N-VECTOR.
C   A IS THE SOLUTION M-VECTOR.
C
C
C
      DIMENSION AA(252,27)
      DIMENSION FAY(8,250),EF(250)
      DIMENSION CT(25,250),F(250),R(250),A(25),BINV(25,25)
      DIMENSION P(350),IC(250),IB(250),Y(250),TH(250),IR(25)
      DIMENSION U(25),V(25),W(25),GINV(25,25),VB(25)
      DIMENSION BV(25),IBOUND(250),ICBAS(25),IRBAS(25),ZC(250)
      DIMENSION F1(50),F2(50),BA(6,15)
      INTEGER S(250)
C
  101 FORMAT(10H1EXAMPLE  ,I3,8H SIZE   ,I3,8H  BY    ,I3)
  102 FORMAT(40H0   R.H.S.   OF THE SYSTEM              )
  103 FORMAT(40H0TRANSPOSE OF COEFFICIENT MATRIX        )
  104 FORMAT(17H0EXECUTION TIME =,F12.5,13HSECONDS      )
  105 FORMAT(9H0L1 NORM=,F10.5,6H,RANK=,I4,6H,ITER=,I4,5H,IND=,I4)
  106 FORMAT(40H0THE ANSWER :  A(I)   OR   X(I)         )
  107 FORMAT(40H0RESIDUALS  R(J)    OR E(J)             )
  108 FORMAT(40H0IC(I)                                  )
  109 FORMAT(40H0IR(I)                                  )
  110 FORMAT(1H ,10F12.5)
  111 FORMAT(1H ,20I4)
  115 FORMAT(9H0L1 NORM=,F10.5,6H,RANK=,F3.0,6H,ITER=,F3.0,5H,IND=,F2.0)
  120 FORMAT(1H ,8F15.5)
  125 FORMAT(5F10.6)
  126 FORMAT(6F10.6)
C
      CALL SUNDER
      PREC=1.E-6
      EPS=1.E-4
      TOLER=1.0E-4
      MM=25
      MM2=MM+2
      MMM=(MM*(MM+3))/2
      NN=250
      NN2=NN+2
      IEXMPL=0
    1 IEXMPL=IEXMPL+1
      GO TO (2,4,5,6,7,8,9,13,17,19,20,21,22,26,27,28,29,32,35,100),
     *IEXMPL
C EXAMPLE 1.
    2 N=201
      DX=0.02
      DO 3 I=1,N
      X1=DX*FLOAT(I-1)
      X2=X1*X1
      X3=X2*X1
      X4=X2*X2
      X5=X3*X2
      X6=X3*X3
      X7=X3*X4
      FAY(1,I)=1.
      FAY(2,I)=X1
      FAY(3,I)=X2
      FAY(4,I)=X3
      FAY(5,I)=X4
      FAY(6,I)=X5
      FAY(7,I)=X6
    3 EF(I)= (EXP(-X1))*SIN(X1)
      M=1
      GO TO 10
    4 M=2
      GO TO 10
    5 M=3
      GO TO 10
    6 M=4
      GO TO 10
    7 M=5
      GO TO 10
    8 M=6
      GO TO 10
    9 M=7
   10 DO 12 J=1,N
      F(J)=EF(J)
        DO 11 I=1,M
   11   CT(I,J)=FAY(I,J)
   12 CONTINUE
      GO TO 90
C EXAMPLE 2.
   13 N=100
      N2=50
      M=6
      READ(1,125)(F1(I), I=1,N2)
      READ(1,125)(F2(I), I=1,N2)
      DO 14 J=1,N2
      F(J)=F1(J)
      JN2=J+N2
   14 F(JN2)=F2(J)
      DO 16 J=1,10
      D=0.1*FLOAT(J)
      DD=D*D
        DO 15 I=1,10
        E=0.1*FLOAT(I)
        EE=E*E
        K=10*(J-1)+I
        CT(1,K)=1.0
        CT(2,K)=E
        CT(3,K)=D
        CT(4,K)=E*D
        CT(5,K)=EE
        CT(6,K)=DD
   15   CONTINUE
   16 CONTINUE
      GO TO 90
C EXAMPLE 3.
   17 N=51
      M=8
      DO 18 I=1,N
      D=0.02*FLOAT(I-1)
      DD=D*D
      DDD=D*DD
      E1=D-0.1
      E2=D-0.2
      E3=D-0.4
      E4=D-0.7
      E13=E1*E1*E1
      E23=E2*E2*E2
      E33=E3*E3*E3
      E43=E4*E4*E4
      CT(1,I)=1.0
      CT(2,I)=D
      CT(3,I)=DD
      CT(4,I)=DDD
      IF(E1.GT.0.0) CT(5,I)=E13
      IF(E1.LE.0.0) CT(5,I)=0.
      IF(E2.GT.0.0) CT(6,I)=E23
      IF(E2.LE.0.0) CT(6,I)=0.
      IF(E3.GT.0.0) CT(7,I)=E33
      IF(E3.LE.0.0) CT(7,I)=0.
      IF(E4.GT.0.0) CT(8,I)=E43
      IF(E4.LE.0.0) CT(8,I)=0.
   18 F(I)= SQRT(D)
      GO TO 90
C EXAMPLE 4.
   19 M=6
      PI=4.0*ATAN(1.0)
      X01=PI/2.0
      X02=X01*X01
      X03=X02*X01
      X04=X02*X02
      N=23
      GO TO 23
   20 N=53
      GO TO 23
   21 N=103
      GO TO 23
   22 N=203
   23 N2=N-2
      N3=N-3
      DX=PI/FLOAT(N3)
      DO 24 I=1,N2
      X1=-X01+DX*FLOAT(I-1)
      X2=X1*X1
      X3=X2*X1
      X4=X2*X2
      X5=X3*X2
      CT(1,I)=1.0
      CT(2,I)=X1
      CT(3,I)=X2
      CT(4,I)=X3
      CT(5,I)=X4
      CT(6,I)=X5
   24 F(I)= SIN(X1)
      G=-1.0
      DO 25 I=1,2
      IF(I.EQ.2) G=1.0
      J=N2+I
      CT(1,J)=0.0
      CT(2,J)=1000.0
      CT(3,J)=2000.0*G*X01
      CT(4,J)=3000.0*X02
      CT(5,J)=4000.0*G*X03
      CT(6,J)=5000.0*X04
   25 F(J)=0.0
      GO TO 90
C EXAMPLE 5.
   26 M=11
      D=20.0
      DX=1.0/D
      N=21
      GO TO 30
   27 D=50.0
      DX=1.0/D
      N=51
      GO TO 30
   28 D=100.0
      DX=1.0/D
      N=101
      GO TO 30
   29 D=200.0
      DX=1.0/D
      N=201
   30 DO 31 I=1,N
      X1=DX*FLOAT(I-1)
      X2=X1+X1
      X3=X1+X2
      X4=X2+X2
      X5=X2+X3
      CT(1,I)=1.0
      CT(2,I)= SIN(X1)
      CT(3,I)= COS(X1)
      CT(4,I)= SIN(X2)
      CT(5,I)= COS(X2)
      CT(6,I)= SIN(X3)
      CT(7,I)= COS(X3)
      CT(8,I)= SIN(X4)
      CT(9,I)= COS(X4)
      CT(10,I)= SIN(X5)
      CT(11,I)= COS(X5)
      G= EXP(X1)
      G1= EXP(0.5)
      F(I)=G1
      IF(G.LT.G1) F(I)=G
   31 CONTINUE
      GO TO 90
C EXAMPLE 6.
   32 M=5
      M1=M+1
      N=15
      DO 34 J=1,N
      READ(1,126)(BA(I,J),I=1,M1)
        DO 33 I=1,M
   33   CT(I,J)=BA(I,J)
      F(J)=BA(M1,J)
   34 CONTINUE
      GO TO 90
C EXAMPLE 7.
   35 M=5
      N=51
      DX=0.02
      DO 36 I=1,N
      X1=DX*FLOAT(I-1)
      X2=X1*X1
      X3=X2*X1
      X4=X2*X2
      CT(1,I)=1.0
      CT(2,I)=X1
      CT(3,I)=X2
      CT(4,I)=X3
      CT(5,I)=X4
      GX=0.0
      IF(X1.GT.0.93) GX=5.0
      F(I)=X1*(1.0+X1*(1.0+X1*(1.0+X1)))+GX
   36 CONTINUE
C
   90 WRITE (3,101)IEXMPL,N,M
C     WRITE(3,102)
C     WRITE(3,110)(F(J),J=1,N)
C     WRITE(3,103)
C     DO 93 I=1,M
C  93 WRITE(3,110)(CT(I,J),J=1,N)
      CALL SETCLK
      CALL L1(M,N,MM,MMM,NN,CT,F,PREC,EPS,IC,IR,IB,
     *U,V,W,Y,TH,P,GINV,VB,IRANK,ITER,R,A,Z,IND)
      CALL RDCLK(TX)
      WRITE(3,108)
      WRITE(3,111)(IC(I),I=1,M)
      WRITE(3,109)
      WRITE(3,111)(IR(I),I=1,M)
      WRITE(3,107)
      WRITE(3,110)(R(J),J=1,N)
      WRITE(3,106)
      WRITE(3,120)(A(I),I=1,M)
      WRITE(3,105)Z,IRANK,ITER,IND
      WRITE(3,104)TX
      GO TO 1
  100 STOP
      END
      SUBROUTINE PSLV(ID,K,MM,MMM,P,B,X)                                PSLV0010
C     THIS SUBROUTINE SOLVES THE SQUARE NON-SINGULAR SYSTEM
C OF LINEAR EQUATIONS
C                            P*X=B,
C OR THE SQUARE NON-SINGULAR SYSTEM OF LINEAR EQUATIONS
C                         P(TRANSPOSE)*X=B,
C WHERE P IS AN UPPER TRIANGULAR MATRIX, B IS THE RIGHT HAND
C SIDE VECTOR AND X IS THE SOLUTION VECTOR.
C     THE INPUT DATA TO THE SUBROUTINE.
C ID     AN INTEGER INDICATOR SPECIFIED BY THE USER.
C        IF ID=1 THE EQUATION P*X=B IS SOLVED.
C        IF ID= ANY INTEGER OTHER THAN 1, THE EQUATION
C        P(TRANSPOSE)*X=B IS SOLVED.
C K      AN INTEGER = THE NUMBER OF EQUATIONS OF THE GIVEN
C        SYSTEM.
C MM     AN INTEGER GREATER THAN OR EQUAL TO K.
C MMM    AN INTEGER = (MM*(MM+3))/2
C P      AN MMM-VECTOR. ITS FIRST (K+1) ELEMENTS CONTAIN THE
C        FIRST K ELEMENTS OF ROW 1 OF THE UPPER TRIANGULAR
C        MATRIX PLUS AN EXTRA LOCATION TO THE RIGHT. ITS
C        NEXT K ELEMENTS CONTAIN THE (K-1) ELEMENTS OF ROW 2
C        OF THE UPPER TRIANGULAR MATRIX PLUS AN EXTRA
C        LOCATION TO THE RIGHT, ..., ETC.
C B      AN MM-VECTOR. ITS FIRST K ELEMENTS CONTAIN THE
C        R.H.S. VECTOR OF THE GIVEN SYSTEM.
C     THE OUTPUT OF THE SUBROUTINE.
C X      AN MM-VECTOR. ON EXIT, ITS FIRST K ELEMENTS CONTAIN
C        THE SOLUTION TO THE GIVEN SYSTEM.
      DOUBLE PRECISION S,SA,SB
      DIMENSION P(MMM),B(MM),X(MM)
      IF(ID.NE.1) GO TO 3
C SOLUTION OF THE UPPER TRIANGULAR SYSTEM.
      L=(K-1)+(K*(K+1))/2
      X(K)=B(K)/P(L)
      IF(K.EQ.1) RETURN
      KD=3
      KM1=K-1
      DO 2 I=1,KM1
        J=K-I
        L=L-KD
        KD=KD+1
        S=B(J)
        LL=L
        JJ=J
        DO 1 KK=1,I
          LL=LL+1
          SA=P(LL)
          JJ=JJ+1
          SB=X(JJ)
    1     S=S-SA*SB
        X(J)=S/P(L)
    2   CONTINUE
      RETURN
C SOLUTION OF THE LOWER TRIANGULAR SYSTEM.
    3 X(1)=B(1)/P(1)
      IF(K.EQ.1) RETURN
      L=1
      KD=K+1
      DO 5 I=2,K
        L=L+KD
        KD=KD-1
        S=B(I)
        KK=I
        KKM1=I-1
        KKD=K
        DO 4 J=1,KKM1
          SA=P(KK)
          SB=X(J)
          S=S-SA*SB
          KK=KK+KKD
          KKD=KKD-1
    4     CONTINUE
    5   X(I)=S/P(L)
      RETURN
      END
      SUBROUTINE L1(M,N,MM,MMM,NN,CT,F,PREC,EPS,IC,IR,IB,               L1  0010
     *UF,BP,XP,T,ALFA,P,GINV,VB,IRANK,ITER,R,A,Z,IND)
C     THIS SUBROUTINE SOLVES AN OVERDETERMINED SYSTEM OF
C LINEAR EQUATIONS IN THE L1 NORM BY USING A DUAL SIMPLEX
C ALGORITHM TO THE LINEAR PROGRAMMING FORMULATION OF THE
C GIVEN PROBLEM. IN THIS ALGORITHM CERTAIN INTERMEDIATE
C SIMPLEX ITERATIONS ARE SKIPPED.
C     FOR PURPOSE OF NUMERICAL STABILITY, THIS SUBROUTINE
C USES A TRIANGULAR DECOMPOSITION TO THE BASIS MATRIX.
C     THE SYSTEM OF LINEAR EQUATIONS HAS THE FORM
C                          C*A=F,
C WHERE C IS A GIVEN REAL N BY M MATRIX OF RANK K.LE.M.LE.N
C AND F IS A GIVEN REAL N-VECTOR.
C     THE PROBLEM TO BE SOLVED IS TO CALCULATE THE ELEMENTS
C OF THE M-VECTOR  A*  WHICH GIVES THE MINIMUM NORM Z.
C               Z= ABS(R(1)) + ABS(R(2)) + ... + ABS(R(N)) ,
C WHERE R(I) IS THE ITH RESIDUAL AND IS GIVEN BY
C   R(I)=C(I,1)*A(1)+C(I,2)*A(2)+ ... +C(I,M)*A(M)-F(I).
C     SUBROUTINE L1 IS COMPLETELY SELF CONTAINED (CONSISTS
C OF TWO SUBROUTINES L1   AND PSLV).
C     THE INPUT DATA TO THE SUBROUTINE.
C M      THE NUMBER OF COLUMNS OF MATRIX C.
C N      THE NUMBER OF ROWS OF MATRIX C.
C MM     AN INTEGER GREATER THAN OR EQUAL TO M.
C MMM    AN INTEGER = (MM*(MM+3))/2 .
C NN     AN INTEGER GREATER THAN OR EQUAL TO N.
C CT     A MATRIX OF DIMENSIONS MM BY NN. ON ENTRY, ITS FIRST
C        M ROWS AND FIRST N COLUMNS CONTAIN THE TRANSPOSE OF
C        MATRIX C IN THE GIVEN SYSTEM C*A=F.
C F      AN NN-VECTOR. ON ENTRY, ITS FIRST N ELEMENTS CONTAIN
C        THE R.H.S. OF THE GIVEN SYSTEM C*A=F.
C PREC   THE ROUND-OFF LEVEL OF THE COMPUTER. FOR THE IBM
C        360/67 COMPUTER, PREC IS ABOUT 1.E-6 AND 1.E-16 FOR
C        SINGLE AND DOUBLE PRECISION RESPECTIVELY.
C EPS    A SPECIFIED TOLERANCE SUCH THAT A CALCULATED NUMBER
C        X IS CONSIDERED ZERO IF  ABS(X) <EPS. FOR THE IBM 360/67
C        COMPUTER, EPS IS USUALLY TAKEN 1.E-4 AND 1.E-11
C        RESPECTIVELY.
C     THE RESULTS OF THE PROBLEM.
C IRANK  THE CALCULATED RANK OF MATRIX C.
C ITER   THE NUMBER OF ITERATIONS WHICH THE SOLUTION NEEDED.
C A      AN MM-VECTOR. ITS FIRST M ELEMENTS ARE THE SOLUTION
C        VECTOR A*.
C R      AN NN-VECTOR. ITS FIRST N ELEMENTS CONTAIN THE
C        THE RESIDUAL VECTOR R=C*A-F.
C Z      THE MINIMUM L1 NORM OF THE RESIDUAL VECTOR R.
C IND    A RETURN INDICATOR. IND=0 INDICATES THAT THE
C        SOLUTION VECTOR A* IS UNIQUE. IND=1, INDICATES THAT
C        A* IS MOST PROBABLY NOT UNIQUE. IND=-1 INDICATES
C        PREMATURE TERMINATION OF THE CALCULATION DUE TO
C        VERY ILL-CONDITIONING OF MATRIX C.
C     THE MEANING OF THE OTHER PARAMETERS.
C GINV   AN MM-SQUARE MATRIX. ITS FIRST IRANK COLUMNS AND
C        FIRST IRANK ROWS CONTAIN THE INVERSE OF THE INITIAL
C        BASIS MATRIX AND ITS UPDATE.
C VB     AN MM-VECTOR. ITS FIRST IRANK ELEMENTS CONTAIN THE
C        INITIAL BASIC SOLUTION AND ITS UPDATE.
C T      AN NN-VECTOR. ITS FIRST N ELEMENTS CONTAIN THE
C        ELEMENTS OF THE ROW IN THE SIMPLEX TABLEAU, THAT
C        CORRESPONDS TO THE COLUMN WHICH LEAVES THE BASIS.
C ALFA   AN NN-VECTOR. ITS FIRST N ELEMENTS CONTAIN THE
C        RATIOS:  ALFA(J) = R(J)/T(J).
C IC     AN NN-VECTOR. ITS FIRST N ELEMENTS CONTAIN THE
C        COLUMN INDICES OF MATRIX CT.
C IR     AN MM VECTOR. ITS FIRST IRANK ELEMENTS CONTAIN THE
C        ROW INDICES OF THE LINEARLY INDEPENDENT ROWS OF CT.
C IB     A SIGN NN-VECTOR. ITS FIRST N ELEMENTS HAVE THE
C        VALUES +1 OR -1. IB(J)=+1 INDICATES THAT COLUMN J
C        IN THE TABLEAU IS AT ITS LOWER BOUND 0. IB(J)=-1
C        INDICATES THAT COLUMN J IS AT ITS UPPER BOUND 2.
C P      AN MMM-VECTOR. ITS FIRST ((IRANK*(IRANK+3))/2)-1
C        ELEMENTS CONTAIN THE (IRANK*(IRANK+1))/2 ELEMENTS OF
C        THE UPPER TRIANGULAR MATRIX + EXTRA (IRANK-1) WORKING
C        LOCATIONS. SEE THE COMMENTS IN SUBROUTINE PSLV.
C BP     AN MM-VECTOR. ITS FIRST IRANK ELEMENTS ARE THE
C        R.H.S. OF THE TRIANGULAR EQUATIONS AS P*XP=BP.
C XP     AN MM-VECTOR . ITS FIRST IRANK ELEMENTS ARE THE
C        SOLUTION OF THE TRIANGULAR EQUATIONS AS P*XP=BP.
C UF     AN MM-WORKING VECTOR.
      DOUBLE PRECISION S,SA,SB
      DIMENSION CT(MM,NN),F(NN),A(MM),GINV(MM,MM),P(MMM)
      DIMENSION IC(NN),IB(NN),R(NN),T(NN),ALFA(NN),IR(MM)
      DIMENSION UF(MM),BP(MM),XP(MM),VB(MM)
      IND=0
      TPEPS=2.+EPS
      NMM=(M*(M+3))/2
      IRANK=M
      ITER=0
      DO 1 J=1,N
        IB(J)=1
    1   IC(J)=J
      DO 3 J=1,M
        IR(J)=J
        A(J)=0.
        DO 2 I=1,M
    2     GINV(I,J)=0.
    3   GINV(J,J)=1.
      IOUT=0
C PART 1 OF THE ALGORITHM.
    4 IOUT=IOUT+1
      IF(IOUT.GT.IRANK) GO TO 16
      PIV=0.
      DO 6 J=IOUT,N
        ICJ=IC(J)
        DO 5 I=IOUT,IRANK
          D=CT(I,ICJ)
          IF(D.LT.0.0) D=-D
          IF(D.LE.PIV) GO TO 5
          LI=I
          JIN=ICJ
          LJ=J
          PIV=D
    5     CONTINUE
    6   CONTINUE
C DETECTION OF RANK DEFICIENCY.
      IF(PIV.GT.EPS) GO TO 7
      IRANK=IOUT-1
      IND=1
      GO TO 16
    7 IF(LI.EQ.IOUT) GO TO 10
      DO 8 J=1,N
        G=CT(LI,J)
        CT(LI,J)=CT(IOUT,J)
    8   CT(IOUT,J)=G
      K=IR(LI)
      IR(LI)=IR(IOUT)
      IR(IOUT)=K
      IF(IOUT.EQ.1) GO TO 10
      K=IOUT-1
      DO 9 J=1,K
        D=GINV(LI,J)
        GINV(LI,J)=GINV(IOUT,J)
    9   GINV(IOUT,J)=D
C A GAUSS-JORDAN ELIMINATION STEP.
   10 PIVOT=CT(IOUT,JIN)
      DO 11 J=1,N
   11   CT(IOUT,J)=CT(IOUT,J)/PIVOT
      DO 12 J=1,IOUT
   12   GINV(IOUT,J)=GINV(IOUT,J)/PIVOT
      DO 15 I=1,IRANK
        IF(I.EQ.IOUT) GO TO 15
        D=CT(I,JIN)
        DO 13 J=1,N
   13     CT(I,J)=CT(I,J)-D*CT(IOUT,J)
        DO 14 J=1,IOUT
   14     GINV(I,J)=GINV(I,J)-D*GINV(IOUT,J)
   15   CONTINUE
      ITER=ITER+1
      K=IC(LJ)
      IC(LJ)=IC(IOUT)
      IC(IOUT)=K
      GO TO 4
C PART 2 OF THE ALGORITHM.
   16 IRANK1=IRANK+1
      IRNKM1=IRANK-1
C INITIAL RESIDUALS AND INITIAL BASIC SOLUTION.
      DO 17 J=1,IRANK
        ICJ=IC(J)
        R(ICJ)=0.
   17   UF(J)=F(ICJ)
      DO 19 J=IRANK1,N
        ICJ=IC(J)
        S=-F(ICJ)
        DO 18 I=1,IRANK
          SA=UF(I)
          SB=CT(I,ICJ)
   18     S=S+SA*SB
        R(ICJ)=S
        IF(S.GE.0.0) GO TO 19
        IB(ICJ)=-1
   19   CONTINUE
      DO 21 I=1,IRANK
        S=0.
        DO 20 J=1,N
          SA=CT(I,J)
          IF(IB(J).EQ.(-1)) SA=-SA
   20     S=S+SA
   21   VB(I)=S
C INITIALIZING THE TRIANGULAR MATRIX.
      DO 22 I=1,NMM
   22   P(I)=0.
      K=1
      KD=IRANK1
      DO 23 I=1,IRANK
        P(K)=1.
        K=K+KD
   23   KD=KD-1
C DETERMINE THE VECTOR WHICH LEAVES THE BASIS.
   24 IVO=0
      CALL PSLV(1,IRANK,MM,MMM,P,VB,XP)
      G=1.
      DO 27 I=1,IRANK
        E=XP(I)
        IF(E.LT.(-EPS)) GO TO 25
        IF(E.LE.TPEPS) GO TO 27
        D=2.0-E
        IF(D.GE.G) GO TO 27
        IVO=1
        GO TO 26
   25   D=E
        IF(D.GE.G) GO TO 27
        IVO=-1
   26   G=D
        IOUT=I
        XB=E
   27   CONTINUE
      IF(IVO.EQ.0) GO TO 66
C CALCULATION OF ROW (IOUT) IN THE TABLEAU.
      ICIOUT=IC(IOUT)
      T(ICIOUT)=1.
      BXB=XB
      DO 28 I=1,IRANK
   28   BP(I)=0.
      BP(IOUT)=1.
      CALL PSLV(2,IRANK,MM,MMM,P,BP,XP)
      ALFMX=0.0
      DO 31 J=IRANK1,N
        ICJ=IC(J)
        ALFA(ICJ)=0.
        S=0.
        DO 29 I=IOUT,IRANK
          SA=XP(I)
          SB=CT(I,ICJ)
   29     S=S+SA*SB
        E=S
        T(ICJ)=E
        IF(E.LT.EPS.AND.E.GT.(-EPS)) GO TO 31
        D=R(ICJ)
        IF(D.NE.0.0) GO TO 30
        D=PREC*PREC*FLOAT(IB(ICJ))
   30   ALFA(ICJ)=D/E
        GG=ALFA(ICJ)
        IF(GG.LT.0.0) GG=-GG
        IF(GG.LE.ALFMX) GO TO 31
        ALFMX=GG
   31   CONTINUE
      PIVOTO=1.
      ITEST=0
C DETERMINE THE VECTOR WHICH ENTERS THE BASIS.
      GG=ALFMX+ALFMX
   32 ALFMX=GG
      ALFMN=-GG
      DO 35 J=IRANK1,N
        ICJ=IC(J)
        E=ALFA(ICJ)
        D=E*FLOAT(IVO)
        IF(D.LE.0.0) GO TO 35
        IF(IVO.EQ.1) GO TO 33
        IF(E.LE.ALFMN) GO TO 35
        ALFMN=E
        GO TO 34
   33   IF(E.GE.ALFMX) GO TO 35
        ALFMX=E
   34   JIN=J
        ITEST=1
   35   CONTINUE
      IF(ITEST.EQ.1) GO TO 36
C NO FEASIBLE SOLUTION HAS BEEN FOUND.
      IND=-1
      GO TO 66
   36 ICJIN=IC(JIN)
      PIVOT=T(ICJIN)
      ALPHA=ALFA(ICJIN)
      PIVOTN=PIVOT/PIVOTO
      IF(XB.GT.TPEPS) GO TO 37
      IF(PIVOTN.GT.0.) GO TO 39
      GO TO 41
   37 DO 38 I=1,IRANK
        E=CT(I,ICIOUT)
   38   VB(I)=VB(I)-E-E
      E=T(ICIOUT)
      BXB=BXB-E-E
      IB(ICIOUT)=-1
      IF(PIVOTN.GT.0.) GO TO 41
   39 DO 40 I=1,IRANK
        E=CT(I,ICJIN)
   40   VB(I)=VB(I)+E+E
      E=T(ICJIN)
      BXB=BXB+E+E
      IB(ICJIN)=1
   41 XB=BXB/PIVOT
      IF(XB.GE.(-EPS).AND.XB.LE.TPEPS) GO TO 42
      ITEST=0
   42 ALFA(ICJIN)=0.
      IF(ITEST.EQ.1) GO TO 43
      PIVOTO=PIVOT
      ICIOUT=ICJIN
      GO TO 32
   43 R(ICJIN)=0.
      ITER=ITER+1
      IF(IOUT.EQ.IRANK) GO TO 46
C UPDATING MATRIX (P,GINV,VB,CT).
      DO 45 J=IOUT,IRNKM1
        K=J
        K1=K+1
        KD=IRANK
        L=IC(K)
        IC(K)=IC(K1)
        IC(K1)=L
        DO 44 I=1,K1
          P(K)=P(K+1)
          K=K+KD
   44     KD=KD-1
   45   CONTINUE
   46 L=IC(IRANK)
      IC(IRANK)=IC(JIN)
      IC(JIN)=L
      K=IRANK
      KD=IRANK
      DO 47 I=1,IRANK
        P(K)=CT(I,ICJIN)
        K=K+KD
   47   KD=KD-1
      IF(IOUT.EQ.IRANK) GO TO 58
      DO 57 I=IOUT,IRNKM1
        II=I
        I1=I+1
        K=0
        KD=IRANK1
        DO 48 J=1,II
          K=K+KD
   48     KD=KD-1
        KK=K
        KL=K-KD
        L=KL
        G=P(K)
        D=P(L)
        IF(G.LT.0.0) G=-G
        IF(D.LT.0.0) D=-D
        IF(G.LE.D) GO TO 53
        DO 49 J=II,IRANK
          E=P(K)
          P(K)=P(L)
          P(L)=E
          K=K+1
   49     L=L+1
        J=IR(I)
        IR(I)=IR(I1)
        IR(I1)=J
        DO 50 J=1,N
          E=CT(I,J)
          CT(I,J)=CT(I1,J)
   50     CT(I1,J)=E
        DO 51 J=1,IRANK
          E=GINV(I,J)
          GINV(I,J)=GINV(I1,J)
   51     GINV(I1,J)=E
        DO 52 J=1,IRANK
          E=GINV(J,I)
          GINV(J,I)=GINV(J,I1)
   52     GINV(J,I1)=E
        E=VB(I)
        VB(I)=VB(I1)
        VB(I1)=E
   53   E=P(KK)/P(KL)
        IF(E.LT.PREC.AND.E.GT.(-PREC)) GO TO 57
        K=KK
        L=KL
        DO 54 J=II,IRANK
          P(K)=P(K)-E*P(L)
          K=K+1
   54     L=L+1
        DO 55 J=1,N
   55     CT(I1,J)=CT(I1,J)-E*CT(I,J)
        DO 56 J=1,IRANK
   56     GINV(I1,J)=GINV(I1,J)-E*GINV(I,J)
        VB(I1)=VB(I1)-E*VB(I)
   57   CONTINUE
   58 DO 59 J=1,IRANK
        ICJ=IC(J)
   59   UF(J)=F(ICJ)
      IF(ALPHA.GT.(1.0).OR.ALPHA.LT.(-1.0)) GO TO 61
      DO 60 J=IRANK1,N
        ICJ=IC(J)
   60   R(ICJ)=R(ICJ)-ALPHA*T(ICJ)
      GO TO 64
   61 CALL PSLV(2,IRANK,MM,MMM,P,UF,XP)
      DO 63 J=IRANK1,N
        ICJ=IC(J)
        S=-F(ICJ)
        DO 62 I=1,IRANK
          SA=XP(I)
          SB=CT(I,ICJ)
   62     S=S+SA*SB
   63   R(ICJ)=S
   64 DO 65 J=IRANK1,N
        ICJ=IC(J)
        D=R(ICJ)*FLOAT(IB(ICJ))
        IF(D.GE.0.0) GO TO 65
        R(ICJ)=0.0
   65   CONTINUE
      GO TO 24
C CALCULATING THE ANSWER OF THE PROBLEM.
   66 CALL PSLV(2,IRANK,MM,MMM,P,UF,VB)
      DO 68 I=1,IRANK
        S=0.
        DO 67 K=1,IRANK
          SA=VB(K)
          SB=GINV(K,I)
   67     S=S+SA*SB
        K=IR(I)
   68   A(K)=S
      S=0.
      DO 69 J=1,N
        SA=R(J)
        IF(SA.LT.0.0) SA=-SA
   69   S=S+SA
      Z=S
      IF(IND.NE.0) RETURN
      E=2.-EPS
      DO 70 I=1,IRANK
        D=XP(I)
        IF(D.LT.EPS.OR.D.GT.E) IND=1
   70   CONTINUE
      RETURN
      END
  1.478689  1.113268  0.855048  0.673991  0.541967
  0.442372  0.365257  0.304338  0.255453  0.215730
  1.027547  0.884824  0.730009  0.598490  0.492769
  0.408461  0.340883  0.286251  0.241692  0.205051
  0.692608  0.652965  0.577505  0.496487  0.422211
  0.357993  0.303718  0.258208  0.220102  0.188148
  0.443503  0.452751  0.427474  0.386389  0.341112
  0.297441  0.257770  0.222784  0.192393  0.166192
  0.251763  0.286357  0.291899  0.279505  0.257866
  0.232626  0.207020  0.182721  0.160481  0.140547
  0.100448  0.149070  0.173775  0.181324  0.177820
  0.167924  0.154817  0.140517  0.126221  0.112594
 -0.020740  0.035871  0.072667  0.093931  0.103914
  0.106235  0.103674  0.098224  0.091240  0.083607
 -0.118431 -0.057358 -0.012925  0.017703  0.037524
  0.049293  0.055315  0.057387  0.056853  0.054674
 -0.197141 -0.133831 -0.084691 -0.047765 -0.020895
 -0.001986  0.010826  0.019096  0.024061  0.026676
 -0.260102 -0.196048 -0.144211 -0.103189 -0.071395
 -0.047222 -0.029174 -0.015940 -0.006419  0.000282
        5.        3.        4.       12.        4.        7.
        9.        7.        3.       19.       13.        4.
        6.        6.        0.       12.       12.        2.
        9.        9.        7.       25.       11.        7.
        3.        0.        1.        4.        2.        7.
        8.        1.        8.       17.        1.        7.
        1.        9.        8.       18.        2.        3.
        0.        9.        3.       12.        6.        3.
        3.        1.        1.        5.        3.        5.
        6.        7.        6.       19.        7.        1.
        6.        1.        9.       16.       -2.        4.
        0.        4.        8.       12.       -4.        1.
        0.        5.        7.       12.       -2.        6.
        7.        3.        2.       12.        8.        6.
        5.        4.        9.       18.        0.        0.
      SUBROUTINE SUNDER                                                 SUN 0010
C     THIS IS A DUMMY SUBROUTINE. ITS FUNCTION IS TO
C SUPPRESS THE UNDERFLOWS OCCURING IN THE CALCULATION.
      RETURN
      END
      SUBROUTINE SETCLK                                                 SET 0010
C     THIS IS A DUMMY SUBROUTINE. ITS FUNCTION IS TO START
C READING THE CPU TIME.
      RETURN
      END
      SUBROUTINE RDCLK(TX)                                              RDC 0010
C     THIS IS A DUMMY SUBROUTINE. ITS FUNCTION IS TO RECORD
C THE CPU TIME  TX  IN SECONDS, WHERE TX IS A REAL VARIABLE.
      TX=0.00
      RETURN
      END