Changes needed by ICC/IFC compiler (Intel)
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / h / dsmplx.F
CommitLineData
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(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
63C LW(I,1), LW(K,2) = 0 x >= 0
64C LW(I,1), LW(K,2) = -1 x = 0
65C LW(I,1), LW(K,2) = 1 x unrestricted
66C LW(I,1), LW(K,2) = 2 x linearly independent, no influence
67C LW(I,1) = -2 x unrestricted, cannot be eliminated
68
69C 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
77C 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
98C Homogeneous part of at least 2 equ. is linearly dependent
99C 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
108C 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
115C If non-basic variables only equation variables,
116C 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
132C 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
141C If there are no free variables, or if they have been
142C exchanged with equations, go to 200.
143C 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
158C If tableau represents an initial solution (IND = 1),
159C 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
212C If all pivot elements for free variable = 0, it has to
213C remain out of the basis (LW(.,1) = -2).
214
215 LW(IP,1)=-2
216 ENDIF
217 GO TO 111
218
219C 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
255C Only unrestricted variables in the basis,
256C 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
297C Variables out of the basis either unrestricted or 0,
298C they cannot be exchanged, hence final tableau obtained.
299C
300C 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
311C All C(K) for constraints X >= 0, hence initial solution found.
312C 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
337C Initial solution found, hence algorithm may stop.
338C 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
352C KASE = 1: There is a free variable out of the basis
353C which does not influence the maximum (ITYPE = 2).
354C 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
364C All B(I) non-negative, hence final tableau obtained
365C 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
410C Some C(K) = 0 and some B(I) = 0 in the final tableau.
411C 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