[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / h / dsmplx.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:49  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
#if defined(CERNLIB_DOUBLE)
      SUBROUTINE DSMPLX(A,B,C,Z0,IDA,M,M1,N,N1,LW,IDW,W,X,Z,ITYPE)
#endif
#if !defined(CERNLIB_DOUBLE)
      SUBROUTINE RSMPLX(A,B,C,Z0,IDA,M,M1,N,N1,LW,IDW,W,X,Z,ITYPE)
#endif
#include "gen/imp64.inc"

      LOGICAL L1,L2,L3
      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
#if defined(CERNLIB_DOUBLE)
      PARAMETER (NAME = 'DSMPLX')
#endif
#if !defined(CERNLIB_DOUBLE)
      PARAMETER (NAME = 'RSMPLX')
#endif
      PARAMETER (R10 = 10, IE0 = 15)

      DIMENSION A(IDA,*),B(*),C(*),W(*),X(*),LW(IDW,*)

      IF(M .LE. 0 .OR. N .LE. 0) THEN
       WRITE(ERRTXT,101) M,N
       CALL MTLPRT(NAME,'H101.1',ERRTXT)
       RETURN
      ENDIF
      IF(.NOT.(0 .LE. M1 .AND. M1 .LE. M .AND.
     1         0 .LE. N1 .AND. N1 .LE. N)) THEN
       WRITE(ERRTXT,102) M,N,M1,N1
       CALL MTLPRT(NAME,'H101.2',ERRTXT)
       RETURN
      ENDIF

      EPS=0
      DO 5 I = 1,M
      DO 5 K = 1,N
    5 EPS=EPS+ABS(A(I,K))
      EPSL=LOG10(2*EPS/(M*N))
      IEXP=INT(EPSL)-IE0
      IF(EPSL .LT. 0) IEXP=IEXP-1
      EPS=R10**IEXP

      Z=Z0
      DO 10 I = 1,M
      LW(I,1)=0
      IF(I .GT. M1) LW(I,1)=1
   10 LW(I,4)=I
      DO 11 K = 1,N
      LW(K,2)=0
      IF(K .GT. N1) LW(K,2)=-1
      LW(K,5)=M+K
   11 LW(K,3)=0

C     LW(I,1), LW(K,2) =  0   x >= 0
C     LW(I,1), LW(K,2) = -1   x = 0
C     LW(I,1), LW(K,2) =  1   x unrestricted
C     LW(I,1), LW(K,2) =  2   x linearly independent, no influence
C     LW(I,1)          = -2   x unrestricted, cannot be eliminated

C     Elimination of the equality constraints from the basis

      IF(N1 .NE. N) THEN
   55  KP=0
       DO 60 K = 1,N
       IF(LW(K,2) .EQ. -1) KP=K
   60  CONTINUE

C     If KP = 0, no more equality constraints in the basis

       IF(KP .NE. 0) THEN
        DMAX=0
        DO 61  I = 1, M
        IF(LW(I,1) .NE. -1) THEN
         DMAX=MAX(ABS(A(I,KP)),DMAX)
         IF(ABS(A(I,KP)) .EQ. DMAX) IP=I
        ENDIF
   61   CONTINUE
        IF(DMAX .LT. EPS) DMAX=0
        IF(DMAX .EQ. 0) THEN
         IF(ABS(C(KP)) .LT. EPS) C(KP)=0
         IF(C(KP) .EQ. 0) THEN
          LW(KP,2)=1
          GO TO 55
         ENDIF
         ITYPE=4
         GO TO 9
        ELSE

C     Homogeneous part of at least 2 equ. is linearly dependent
C     Inhomogeneous part is contradictory.

         CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
         LW(KP,2)=LW(IP,1)
         LW(IP,1)=-1
         GO TO 55
        ENDIF
       ENDIF

C     Non-basic variables only equation variables?

       IND=0
       DO 81 I = 1,M
       IF(LW(I,1) .NE. -1)  IND=IND+1
   81  CONTINUE

C     If non-basic variables only equation variables,
C     is constraint violated?

       IF(IND .EQ. 0) THEN
        DMIN=0
        DO 85  K = 1, N
        IF(LW(K,2) .NE. 1) THEN
         IF(ABS(C(K)) .LT. EPS) C(K)=0
         DMIN=MIN(C(K),DMIN)
        ENDIF
   85   CONTINUE
        ITYPE=1
        IF(DMIN .LT. 0) ITYPE=4
        GO TO 9
       ENDIF
      ENDIF

C     Eliminate unrestricted variables

      IF(M1 .EQ. M) GO TO 200
  111 IP=0
      DO 105 I = 1,M
      IF(LW(I,1) .EQ. 1) IP=I
  105 CONTINUE
      IF(IP .EQ. 0) GO TO 200

C     If there are no free variables, or if they have been
C     exchanged with equations, go to 200.
C     Present tableau represents feasable initial solution?

      DMAX=0
      DMIN=0
      DO 110  K = 1, N
      IF(LW(K,2) .NE. 1) THEN
       IF(ABS(C(K)) .LT. EPS) C(K)=0
       DMIN=MIN(C(K),DMIN)
       DMAX=MAX(C(K),DMAX)
      ENDIF
  110 CONTINUE
      IND=0
      IF(DMIN .GE. 0) IND=1
      Q=DMAX

C     If tableau represents an initial solution (IND = 1),
C     elimination of remaining unrestricted variables.

      DMAX=0
      DMIN=0
      DO 112  K = 1, N
      IF(LW(K,2) .NE. 1) THEN
       IF(ABS(A(IP,K)) .LT. EPS) A(IP,K)=0
       DMAX=MAX(A(IP,K),DMAX)
       IF(DMAX .EQ. A(IP,K)) KPMAX=K
       DMIN=MIN(A(IP,K),DMIN)
       IF(DMIN .EQ. A(IP,K)) KPMIN=K
      ENDIF
  112 CONTINUE
      IF(DMAX .NE. 0 .OR. DMIN .NE. 0) THEN
       IF(IND .NE. 0) THEN
        IF(ABS(B(IP)) .LT. EPS)  B(IP)=0
        IF(B(IP) .LT. 0 .AND. DMAX .GT. 0) THEN
         Q=Q/DMAX
         DO 150 K = 1,N
         IF(LW(K,2) .NE. 1 .AND. A(IP,K) .GT. 0) THEN
          QQ=C(K)/A(IP,K)
          Q=MIN(QQ,Q)
          IF(Q .EQ. QQ) KP=K
         ENDIF
  150    CONTINUE
        ELSEIF(B(IP) .GT. 0 .AND. DMIN .LT. 0) THEN
         Q=Q/DMIN
         DO 170 K = 1,N
         IF(LW(K,2) .NE. 1 .AND. A(IP,K) .LT. 0) THEN
          QQ=C(K)/A(IP,K)
          Q=MAX(QQ,Q)
          IF(Q .EQ. QQ) KP=K
         ENDIF
  170    CONTINUE
        ELSEIF(B(IP) .EQ. 0) THEN
         Q=Q/MAX(DMAX,ABS(DMIN))
         DO 190 K = 1,N
         IF(LW(K,2) .NE. 1 .AND. A(IP,K) .NE. 0) THEN
          QQ=ABS(C(K)/A(IP,K))
          Q=MIN(QQ,Q)
          IF(Q .EQ. QQ) KP=K
         ENDIF
  190    CONTINUE
        ENDIF
       ELSE
        KP=KPMAX
        IF(DMAX .LT. ABS(DMIN)) KP=KPMIN
       ENDIF
       CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
       LW(KP,2)=1
       LW(IP,1)=0
      ELSE

C     If all pivot elements for free variable = 0, it has to
C     remain out of the basis (LW(.,1) = -2).

       LW(IP,1)=-2
      ENDIF
      GO TO 111

C     Feasable initial solution by mutiphase method.

  200 IF(M1 .EQ. M .AND. N1 .EQ. N) GO TO 250
      KP=0
      DO 201 K = 1,N
      IF(LW(K,2) .EQ. 0) KP=KP+1
  201 CONTINUE
      IF(KP .EQ. 0) THEN
       IP=0
       DO 202 I = 1,M
       IF(LW(I,1) .NE. -1) IP=IP+1
  202  CONTINUE
       IF(IP .EQ. 0) THEN
        ITYPE=1
        GO TO 9
       ENDIF

       L1=.FALSE.
       L2=.FALSE.
       L3=.FALSE.
       DO 206 I = 1,M
       IF(LW(I,1) .NE. -1) THEN
        IF(ABS(B(I)) .LT. EPS) B(I)=0
        L1=LW(I,1) .EQ.  0 .AND. B(I) .GT. 0
        L2=LW(I,1) .EQ.  0 .AND. B(I) .EQ. 0
        L3=LW(I,1) .EQ.  0 .AND. B(I) .LT. 0
        L2=LW(I,1) .EQ. -2 .AND. B(I) .EQ. 0
        L3=LW(I,1) .EQ. -2 .AND. B(I) .NE. 0
       ENDIF
  206  CONTINUE
       IF(L1) ITYPE=1
       IF(L2) ITYPE=2
       IF(L3) ITYPE=3
       GO TO 9
      ENDIF

C     Only unrestricted variables in the basis,
C     hence final tableau obtained.

      IP=0
      DO 211 I = 1,M
      IF(LW(I,1) .EQ. 0)  IP=IP+1
  211 CONTINUE
      IF(IP .NE. 0) GO TO 250
      DMIN=0
      DO 212 K = 1,N
      IF(LW(K,2) .NE. 1) THEN
       IF(ABS(C(K)) .LT. EPS) C(K)=0
       DMIN=MIN(C(K),DMIN)
      ENDIF
  212 CONTINUE
      IF(DMIN .LT. 0) THEN
       ITYPE=4
       GO TO 9
      ENDIF

      IP=0
      DO 216 I = 1,M
      IF(LW(I,1) .EQ. -2) IP=IP+1
  216 CONTINUE
      IF(IP .EQ. 0) THEN
       ITYPE=1
       GO TO 9
      ENDIF

      L2=.FALSE.
      L3=.FALSE.
      DO 225 I = 1,M
      IF(LW(I,1) .NE. -1) THEN
       IF(ABS(B(I)) .LT. EPS) B(I)=0
       L2=B(I) .EQ. 0
       L3=.NOT.L2
      ENDIF
  225 CONTINUE
      IF(L2) ITYPE=2
      IF(L3) ITYPE=3
      GO TO 9

C     Variables out of the basis either unrestricted or 0,
C     they cannot be exchanged, hence final tableau obtained.
C
C     Tableau representing initial solution?

  250 DMIN=0
      DO 255 K = 1,N
      IF(LW(K,2) .NE. 1) THEN
       IF(ABS(C(K)) .LT. EPS) C(K)=0
       DMIN=MIN(C(K),DMIN)
       IF(C(K) .EQ. DMIN) KK=K
      ENDIF
  255 CONTINUE

C     All C(K) for constraints X >= 0, hence initial solution found.
C     Otherwise column with index KK new objective function.

      IF(DMIN .LT. 0) THEN
  270  DMIN=0
       DO 300 I = 1,M
       IF(LW(I,1) .EQ. 0) THEN
        IF(ABS(A(I,KK)) .LT. EPS) A(I,KK)=0
        DMIN=MIN(A(I,KK),DMIN)
        IF(A(I,KK) .EQ. DMIN) IP=I
       ENDIF
  300  CONTINUE
       IF(DMIN .GE. 0) THEN
        ITYPE=4
        GO TO 9
       ENDIF

       CALL H101S2(A,B,C,M,N,IDA,IP,KP,KK,LW,IDW,W,X,EPS,0)
       IF(KP .EQ. 0) KP=KK
       CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)

       IF(ABS(C(KK)) .LT. EPS) C(KK)=0
       IF(C(KK) .LT. 0) GO TO 270
       GO TO 250
      ENDIF

C     Initial solution found, hence algorithm may stop.
C     Maximum obtained?

  500 KASE=0
      DO 510 I = 1,M
      IF(LW(I,1) .EQ. -2) THEN
       IF(ABS(B(I)) .LT. EPS) B(I)=0
       IF(B(I) .NE. 0) THEN
        ITYPE=3
        GO TO 9
       ENDIF
       KASE=1
      ENDIF
  510 CONTINUE

C     KASE = 1: There is a free variable out of the basis
C     which does not influence the maximum (ITYPE = 2).
C     Take out negative B(I).

      IP=0
      DO 555 I = 1,M
      IF(LW(I,1) .EQ. 0) THEN
       IF(ABS(B(I)) .LT. EPS) B(I)=0
       IF(B(I) .LT. 0)  IP=I
      ENDIF
  555 CONTINUE

C     All B(I) non-negative, hence final tableau obtained
C     Solution unique?

      IF(IP .NE. 0) THEN
       CALL H101S2(A,B,C,M,N,IDA,IP,KP,0,LW,IDW,W,X,EPS,0)
       IF(KP .EQ. 0) THEN
        ITYPE=3
        GO TO 9
       ENDIF
       CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
       GO TO 500
      ENDIF

      IF(KASE .NE. 0) THEN
       ITYPE=2
       GO TO 9
      ENDIF
      IND=0
      DO 562 K = 1,N
      IF(LW(K,2) .EQ. 0) THEN
       IF(ABS(C(K)) .LT. EPS) C(K)=0
       IF(C(K) .EQ. 0) THEN
        IND=IND+1
        LW(K,2)=2
       ENDIF
      ENDIF
  562 CONTINUE
      KASE=IND
      IND=0
      DO 565 I = 1,M
      IF(LW(I,1) .EQ. 0) THEN
      IF(ABS(B(I)) .LT. EPS)  B(I)=0
       IF(B(I) .EQ. 0) THEN
        IND=IND+1
        LW(I,1)=2
       ENDIF
      ENDIF
  565 CONTINUE
      IF(IND .EQ. 0) THEN
       ITYPE=1
       GO TO 9
      ELSEIF(KASE .EQ. 0) THEN
       ITYPE=2
       GO TO 9
      ENDIF

C     Some C(K) = 0 and some B(I) = 0 in the final tableau.
C     Solution unique?

  575 DO 577 I = 1,M
      IF(LW(I,1) .EQ. 2) THEN
       DMAX=0
       DO 576 K = 1, N
       IF(LW(K,2) .EQ. 2) THEN
        IF(ABS(A(I,K)) .LT. EPS) A(I,K)=0
        DMAX=MAX(A(I,K),DMAX)
       ENDIF
  576  CONTINUE
       IF(DMAX .LE. 0) THEN
        ITYPE=2
        GO TO 9
       ENDIF
      ENDIF
  577 CONTINUE
      DO 581 K = 1,N
      IF(LW(K,2) .EQ. 2) THEN
       DMIN=1
       DO 578 I = 1,M
       IF(LW(I,1) .EQ. 2) DMIN=MIN(A(I,K),DMIN)
  578  CONTINUE
       IF(DMIN .GT. 0) THEN
        ITYPE=1
        GO TO 9
       ELSEIF(DMIN .EQ. 0) THEN
        DO 580  I = 1, M
        IF(LW(I,1) .EQ. 2 .AND. A(I,K) .GT. 0) LW(I,1)=-2
  580   CONTINUE
        LW(K,2)=1
       ENDIF
      ENDIF
  581 CONTINUE
      NC=0
      DO 582 K = 1,N
      IF(LW(K,2) .EQ. 2) NC=NC+1
  582 CONTINUE
      NR=0
      DO 583 I = 1,M
      IF(LW(I,1) .EQ. 2) THEN
       NR=NR+1
       IP=I
      ENDIF
  583 CONTINUE
      IF(NR .EQ. 0 .OR. NC .EQ. 0) THEN
       ITYPE=1
       GO TO 9
      ENDIF

      CALL H101S2(A,B,C,M,N,IDA,IP,KP,0,LW,IDW,W,X,EPS,2)

      IF(KP .EQ. 0) THEN
       ITYPE=2
       GO TO 9
      ENDIF

      CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
      GO TO 575

    9 DO 1 I = 1,M
    1 X(LW(I,4))=0
      DO 2 K = 1,N
    2 X(LW(K,5))=C(K)
      RETURN
  101 FORMAT('M = ',I4,'  OR  N = ',I4,'  ILLEGAL')
  102 FORMAT('M = ',I4,',  N = ',I4,':  M1 = ',I4,'  OR  N1 = ',I4,
     1       '  ILLEGAL')
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:02:49 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	#if defined(CERNLIB_DOUBLE)
	11	SUBROUTINE DSMPLX(A,B,C,Z0,IDA,M,M1,N,N1,LW,IDW,W,X,Z,ITYPE)
	12	#endif
	13	#if !defined(CERNLIB_DOUBLE)
	14	SUBROUTINE RSMPLX(A,B,C,Z0,IDA,M,M1,N,N1,LW,IDW,W,X,Z,ITYPE)
	15	#endif
	16	#include "gen/imp64.inc"
	17
	18	LOGICAL L1,L2,L3
	19	CHARACTER NAME()
	20	CHARACTER*80 ERRTXT
	21	#if defined(CERNLIB_DOUBLE)
	22	PARAMETER (NAME = 'DSMPLX')
	23	#endif
	24	#if !defined(CERNLIB_DOUBLE)
	25	PARAMETER (NAME = 'RSMPLX')
	26	#endif
	27	PARAMETER (R10 = 10, IE0 = 15)
	28
	29	DIMENSION A(IDA,),B(),C(),W(),X(),LW(IDW,)
	30
	31	IF(M .LE. 0 .OR. N .LE. 0) THEN
	32	WRITE(ERRTXT,101) M,N
	33	CALL MTLPRT(NAME,'H101.1',ERRTXT)
	34	RETURN
	35	ENDIF
	36	IF(.NOT.(0 .LE. M1 .AND. M1 .LE. M .AND.
	37	1 0 .LE. N1 .AND. N1 .LE. N)) THEN
	38	WRITE(ERRTXT,102) M,N,M1,N1
	39	CALL MTLPRT(NAME,'H101.2',ERRTXT)
	40	RETURN
	41	ENDIF
	42
	43	EPS=0
	44	DO 5 I = 1,M
	45	DO 5 K = 1,N
	46	5 EPS=EPS+ABS(A(I,K))
	47	EPSL=LOG10(2EPS/(MN))
	48	IEXP=INT(EPSL)-IE0
	49	IF(EPSL .LT. 0) IEXP=IEXP-1
	50	EPS=R10**IEXP
	51
	52	Z=Z0
	53	DO 10 I = 1,M
	54	LW(I,1)=0
	55	IF(I .GT. M1) LW(I,1)=1
	56	10 LW(I,4)=I
	57	DO 11 K = 1,N
	58	LW(K,2)=0
	59	IF(K .GT. N1) LW(K,2)=-1
	60	LW(K,5)=M+K
	61	11 LW(K,3)=0
	62
	63	C LW(I,1), LW(K,2) = 0 x >= 0
	64	C LW(I,1), LW(K,2) = -1 x = 0
65	C LW(I,1), LW(K,2) = 1 x unrestricted
66	C LW(I,1), LW(K,2) = 2 x linearly independent, no influence
67	C LW(I,1) = -2 x unrestricted, cannot be eliminated
68
69	C Elimination of the equality constraints from the basis
70
71	IF(N1 .NE. N) THEN
72	55 KP=0
73	DO 60 K = 1,N
74	IF(LW(K,2) .EQ. -1) KP=K
75	60 CONTINUE
76
77	C If KP = 0, no more equality constraints in the basis
78
79	IF(KP .NE. 0) THEN
80	DMAX=0
81	DO 61 I = 1, M
82	IF(LW(I,1) .NE. -1) THEN
83	DMAX=MAX(ABS(A(I,KP)),DMAX)
84	IF(ABS(A(I,KP)) .EQ. DMAX) IP=I
85	ENDIF
86	61 CONTINUE
87	IF(DMAX .LT. EPS) DMAX=0
88	IF(DMAX .EQ. 0) THEN
89	IF(ABS(C(KP)) .LT. EPS) C(KP)=0
90	IF(C(KP) .EQ. 0) THEN
91	LW(KP,2)=1
92	GO TO 55
93	ENDIF
94	ITYPE=4
95	GO TO 9
96	ELSE
97
98	C Homogeneous part of at least 2 equ. is linearly dependent
99	C Inhomogeneous part is contradictory.
100
101	CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
102	LW(KP,2)=LW(IP,1)
103	LW(IP,1)=-1
104	GO TO 55
105	ENDIF
106	ENDIF
107
108	C Non-basic variables only equation variables?
109
110	IND=0
111	DO 81 I = 1,M
112	IF(LW(I,1) .NE. -1) IND=IND+1
113	81 CONTINUE
114
115	C If non-basic variables only equation variables,
116	C is constraint violated?
117
118	IF(IND .EQ. 0) THEN
119	DMIN=0
120	DO 85 K = 1, N
121	IF(LW(K,2) .NE. 1) THEN
122	IF(ABS(C(K)) .LT. EPS) C(K)=0
123	DMIN=MIN(C(K),DMIN)
124	ENDIF
125	85 CONTINUE
126	ITYPE=1
127	IF(DMIN .LT. 0) ITYPE=4
128	GO TO 9
129	ENDIF
130	ENDIF
131
132	C Eliminate unrestricted variables
133
134	IF(M1 .EQ. M) GO TO 200
135	111 IP=0
136	DO 105 I = 1,M
137	IF(LW(I,1) .EQ. 1) IP=I
138	105 CONTINUE
139	IF(IP .EQ. 0) GO TO 200
140
141	C If there are no free variables, or if they have been
142	C exchanged with equations, go to 200.
143	C Present tableau represents feasable initial solution?
144
145	DMAX=0
146	DMIN=0
147	DO 110 K = 1, N
148	IF(LW(K,2) .NE. 1) THEN
149	IF(ABS(C(K)) .LT. EPS) C(K)=0
150	DMIN=MIN(C(K),DMIN)
151	DMAX=MAX(C(K),DMAX)
152	ENDIF
153	110 CONTINUE
154	IND=0
155	IF(DMIN .GE. 0) IND=1
156	Q=DMAX
157
158	C If tableau represents an initial solution (IND = 1),
159	C elimination of remaining unrestricted variables.
160
161	DMAX=0
162	DMIN=0
163	DO 112 K = 1, N
164	IF(LW(K,2) .NE. 1) THEN
165	IF(ABS(A(IP,K)) .LT. EPS) A(IP,K)=0
166	DMAX=MAX(A(IP,K),DMAX)
167	IF(DMAX .EQ. A(IP,K)) KPMAX=K
168	DMIN=MIN(A(IP,K),DMIN)
169	IF(DMIN .EQ. A(IP,K)) KPMIN=K
170	ENDIF
171	112 CONTINUE
172	IF(DMAX .NE. 0 .OR. DMIN .NE. 0) THEN
173	IF(IND .NE. 0) THEN
174	IF(ABS(B(IP)) .LT. EPS) B(IP)=0
175	IF(B(IP) .LT. 0 .AND. DMAX .GT. 0) THEN
176	Q=Q/DMAX
177	DO 150 K = 1,N
178	IF(LW(K,2) .NE. 1 .AND. A(IP,K) .GT. 0) THEN
179	QQ=C(K)/A(IP,K)
180	Q=MIN(QQ,Q)
181	IF(Q .EQ. QQ) KP=K
182	ENDIF
183	150 CONTINUE
184	ELSEIF(B(IP) .GT. 0 .AND. DMIN .LT. 0) THEN
185	Q=Q/DMIN
186	DO 170 K = 1,N
187	IF(LW(K,2) .NE. 1 .AND. A(IP,K) .LT. 0) THEN
188	QQ=C(K)/A(IP,K)
189	Q=MAX(QQ,Q)
190	IF(Q .EQ. QQ) KP=K
191	ENDIF
192	170 CONTINUE
193	ELSEIF(B(IP) .EQ. 0) THEN
194	Q=Q/MAX(DMAX,ABS(DMIN))
195	DO 190 K = 1,N
196	IF(LW(K,2) .NE. 1 .AND. A(IP,K) .NE. 0) THEN
197	QQ=ABS(C(K)/A(IP,K))
198	Q=MIN(QQ,Q)
199	IF(Q .EQ. QQ) KP=K
200	ENDIF
201	190 CONTINUE
202	ENDIF
203	ELSE
204	KP=KPMAX
205	IF(DMAX .LT. ABS(DMIN)) KP=KPMIN
206	ENDIF
207	CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
208	LW(KP,2)=1
209	LW(IP,1)=0
210	ELSE
211
212	C If all pivot elements for free variable = 0, it has to
213	C remain out of the basis (LW(.,1) = -2).
214
215	LW(IP,1)=-2
216	ENDIF
217	GO TO 111
218
219	C Feasable initial solution by mutiphase method.
220
221	200 IF(M1 .EQ. M .AND. N1 .EQ. N) GO TO 250
222	KP=0
223	DO 201 K = 1,N
224	IF(LW(K,2) .EQ. 0) KP=KP+1
225	201 CONTINUE
226	IF(KP .EQ. 0) THEN
227	IP=0
228	DO 202 I = 1,M
229	IF(LW(I,1) .NE. -1) IP=IP+1
230	202 CONTINUE
231	IF(IP .EQ. 0) THEN
232	ITYPE=1
233	GO TO 9
234	ENDIF
235
236	L1=.FALSE.
237	L2=.FALSE.
238	L3=.FALSE.
239	DO 206 I = 1,M
240	IF(LW(I,1) .NE. -1) THEN
241	IF(ABS(B(I)) .LT. EPS) B(I)=0
242	L1=LW(I,1) .EQ. 0 .AND. B(I) .GT. 0
243	L2=LW(I,1) .EQ. 0 .AND. B(I) .EQ. 0
244	L3=LW(I,1) .EQ. 0 .AND. B(I) .LT. 0
245	L2=LW(I,1) .EQ. -2 .AND. B(I) .EQ. 0
246	L3=LW(I,1) .EQ. -2 .AND. B(I) .NE. 0
247	ENDIF
248	206 CONTINUE
249	IF(L1) ITYPE=1
250	IF(L2) ITYPE=2
251	IF(L3) ITYPE=3
252	GO TO 9
253	ENDIF
254
255	C Only unrestricted variables in the basis,
256	C hence final tableau obtained.
257
258	IP=0
259	DO 211 I = 1,M
260	IF(LW(I,1) .EQ. 0) IP=IP+1
261	211 CONTINUE
262	IF(IP .NE. 0) GO TO 250
263	DMIN=0
264	DO 212 K = 1,N
265	IF(LW(K,2) .NE. 1) THEN
266	IF(ABS(C(K)) .LT. EPS) C(K)=0
267	DMIN=MIN(C(K),DMIN)
268	ENDIF
269	212 CONTINUE
270	IF(DMIN .LT. 0) THEN
271	ITYPE=4
272	GO TO 9
273	ENDIF
274
275	IP=0
276	DO 216 I = 1,M
277	IF(LW(I,1) .EQ. -2) IP=IP+1
278	216 CONTINUE
279	IF(IP .EQ. 0) THEN
280	ITYPE=1
281	GO TO 9
282	ENDIF
283
284	L2=.FALSE.
285	L3=.FALSE.
286	DO 225 I = 1,M
287	IF(LW(I,1) .NE. -1) THEN
288	IF(ABS(B(I)) .LT. EPS) B(I)=0
289	L2=B(I) .EQ. 0
290	L3=.NOT.L2
291	ENDIF
292	225 CONTINUE
293	IF(L2) ITYPE=2
294	IF(L3) ITYPE=3
295	GO TO 9
296
297	C Variables out of the basis either unrestricted or 0,
298	C they cannot be exchanged, hence final tableau obtained.
299	C
300	C Tableau representing initial solution?
301
302	250 DMIN=0
303	DO 255 K = 1,N
304	IF(LW(K,2) .NE. 1) THEN
305	IF(ABS(C(K)) .LT. EPS) C(K)=0
306	DMIN=MIN(C(K),DMIN)
307	IF(C(K) .EQ. DMIN) KK=K
308	ENDIF
309	255 CONTINUE
310
311	C All C(K) for constraints X >= 0, hence initial solution found.
312	C Otherwise column with index KK new objective function.
313
314	IF(DMIN .LT. 0) THEN
315	270 DMIN=0
316	DO 300 I = 1,M
317	IF(LW(I,1) .EQ. 0) THEN
318	IF(ABS(A(I,KK)) .LT. EPS) A(I,KK)=0
319	DMIN=MIN(A(I,KK),DMIN)
320	IF(A(I,KK) .EQ. DMIN) IP=I
321	ENDIF
322	300 CONTINUE
323	IF(DMIN .GE. 0) THEN
324	ITYPE=4
325	GO TO 9
326	ENDIF
327
328	CALL H101S2(A,B,C,M,N,IDA,IP,KP,KK,LW,IDW,W,X,EPS,0)
329	IF(KP .EQ. 0) KP=KK
330	CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
331
332	IF(ABS(C(KK)) .LT. EPS) C(KK)=0
333	IF(C(KK) .LT. 0) GO TO 270
334	GO TO 250
335	ENDIF
336
337	C Initial solution found, hence algorithm may stop.
338	C Maximum obtained?
339
340	500 KASE=0
341	DO 510 I = 1,M
342	IF(LW(I,1) .EQ. -2) THEN
343	IF(ABS(B(I)) .LT. EPS) B(I)=0
344	IF(B(I) .NE. 0) THEN
345	ITYPE=3
346	GO TO 9
347	ENDIF
348	KASE=1
349	ENDIF
350	510 CONTINUE
351
352	C KASE = 1: There is a free variable out of the basis
353	C which does not influence the maximum (ITYPE = 2).
354	C Take out negative B(I).
355
356	IP=0
357	DO 555 I = 1,M
358	IF(LW(I,1) .EQ. 0) THEN
359	IF(ABS(B(I)) .LT. EPS) B(I)=0
360	IF(B(I) .LT. 0) IP=I
361	ENDIF
362	555 CONTINUE
363
364	C All B(I) non-negative, hence final tableau obtained
365	C Solution unique?
366
367	IF(IP .NE. 0) THEN
368	CALL H101S2(A,B,C,M,N,IDA,IP,KP,0,LW,IDW,W,X,EPS,0)
369	IF(KP .EQ. 0) THEN
370	ITYPE=3
371	GO TO 9
372	ENDIF
373	CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
374	GO TO 500
375	ENDIF
376
377	IF(KASE .NE. 0) THEN
378	ITYPE=2
379	GO TO 9
380	ENDIF
381	IND=0
382	DO 562 K = 1,N
383	IF(LW(K,2) .EQ. 0) THEN
384	IF(ABS(C(K)) .LT. EPS) C(K)=0
385	IF(C(K) .EQ. 0) THEN
386	IND=IND+1
387	LW(K,2)=2
388	ENDIF
389	ENDIF
390	562 CONTINUE
391	KASE=IND
392	IND=0
393	DO 565 I = 1,M
394	IF(LW(I,1) .EQ. 0) THEN
395	IF(ABS(B(I)) .LT. EPS) B(I)=0
396	IF(B(I) .EQ. 0) THEN
397	IND=IND+1
398	LW(I,1)=2
399	ENDIF
400	ENDIF
401	565 CONTINUE
402	IF(IND .EQ. 0) THEN
403	ITYPE=1
404	GO TO 9
405	ELSEIF(KASE .EQ. 0) THEN
406	ITYPE=2
407	GO TO 9
408	ENDIF
409
410	C Some C(K) = 0 and some B(I) = 0 in the final tableau.
411	C Solution unique?
412
413	575 DO 577 I = 1,M
414	IF(LW(I,1) .EQ. 2) THEN
415	DMAX=0
416	DO 576 K = 1, N
417	IF(LW(K,2) .EQ. 2) THEN
418	IF(ABS(A(I,K)) .LT. EPS) A(I,K)=0
419	DMAX=MAX(A(I,K),DMAX)
420	ENDIF
421	576 CONTINUE
422	IF(DMAX .LE. 0) THEN
423	ITYPE=2
424	GO TO 9
425	ENDIF
426	ENDIF
427	577 CONTINUE
428	DO 581 K = 1,N
429	IF(LW(K,2) .EQ. 2) THEN
430	DMIN=1
431	DO 578 I = 1,M
432	IF(LW(I,1) .EQ. 2) DMIN=MIN(A(I,K),DMIN)
433	578 CONTINUE
434	IF(DMIN .GT. 0) THEN
435	ITYPE=1
436	GO TO 9
437	ELSEIF(DMIN .EQ. 0) THEN
438	DO 580 I = 1, M
439	IF(LW(I,1) .EQ. 2 .AND. A(I,K) .GT. 0) LW(I,1)=-2
440	580 CONTINUE
441	LW(K,2)=1
442	ENDIF
443	ENDIF
444	581 CONTINUE
445	NC=0
446	DO 582 K = 1,N
447	IF(LW(K,2) .EQ. 2) NC=NC+1
448	582 CONTINUE
449	NR=0
450	DO 583 I = 1,M
451	IF(LW(I,1) .EQ. 2) THEN
452	NR=NR+1
453	IP=I
454	ENDIF
455	583 CONTINUE
456	IF(NR .EQ. 0 .OR. NC .EQ. 0) THEN
457	ITYPE=1
458	GO TO 9
459	ENDIF
460
461	CALL H101S2(A,B,C,M,N,IDA,IP,KP,0,LW,IDW,W,X,EPS,2)
462
463	IF(KP .EQ. 0) THEN
464	ITYPE=2
465	GO TO 9
466	ENDIF
467
468	CALL H101S1(A,B,C,Z,M,N,IDA,IP,KP,LW,IDW,EPS)
469	GO TO 575
470
471	9 DO 1 I = 1,M
472	1 X(LW(I,4))=0
473	DO 2 K = 1,N
474	2 X(LW(K,5))=C(K)
475	RETURN
476	101 FORMAT('M = ',I4,' OR N = ',I4,' ILLEGAL')
477	102 FORMAT('M = ',I4,', N = ',I4,': M1 = ',I4,' OR N1 = ',I4,
478	1 ' ILLEGAL')
479	END