This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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