MINICERN/mathlib/gen/h/dsmplx.F

   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(2*EPS/(M*N))
  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