This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / f / lihoin64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:40  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_DOUBLE)
11       SUBROUTINE DLHOIN(A,MA,M,N,MAXV,V,NV,NVEC,EPS,IOUT,W,IW)
12
13 C     Solving Systems of Homogeneous Linear Inequalities.
14 C     Based on K.S. Koelbig and F. Schwarz, A Programm for Solving
15 C     Systems of Homogeneous Linear Inequalities,
16 C     Computer Phys. Comm. 17 (1979) 375-382
17
18 #include "gen/imp64.inc"
19       CHARACTER TEXT*(*)
20       CHARACTER NAME*(*)
21       CHARACTER*80 ERRTXT
22       PARAMETER (NAME = 'DLHOIN')
23       PARAMETER (TEXT = '+++++ CERN F500 DLHOIN ... ')
24
25       DIMENSION A(MA,*),V(NV,*),IW(MA,*),W(MAXV,*)
26
27       IF(MAXV .LT. N) THEN
28        CALL MTLPRT(NAME,'F500.1','MAXV TOO SMALL')
29        RETURN
30       END IF
31 C
32 C*****SETS INITIAL VALUES FOR BOOKKEEPING
33 C
34       DO 3 I = 1,M
35       IW(I,1)=0
36     3 IW(I,2)=I
37       DO 8 I = 1,N
38       IW(I,2)=0
39     8 IW(I,5)=I
40       NINC=N
41       NVEC=N
42 C
43 C*****DETERMINES N BASIS VECTORS OF THE INITIAL POLYHEDRAL CONE
44 C
45       CALL DMCPY(N,N,A(1,1),A(1,2),A(2,1),V(1,1),V(1,2),V(2,1))
46       CALL DINV(N,V,NV,W,IFAIL)
47       IF(IFAIL .EQ. -1) THEN
48        CALL MTLPRT(NAME,'F500.2','MATRIX A(N,N) SINGULAR')
49        RETURN
50       END IF
51       DO 1 I = 1,NVEC
52       S=DVMPY(N,V(1,I),V(2,I),V(1,I),V(2,I))
53     1 CALL DVSCL(N,1/SQRT(S),V(1,I),V(2,I),V(1,I),V(2,I))
54
55       IF(IOUT .EQ. 1) THEN
56        WRITE(6,111) TEXT
57        DO 80 I = 1,NVEC,7
58        DO 81 J = 1,N
59    81  WRITE(6,'(1X,I9,7E15.6)') J,(V(J,I1),I1=I,MIN(NVEC,I+6))
60    80  WRITE(6,'(1X)')
61       END IF
62 C
63 C*****COMPUTES MATRIX OF SCALAR PRODUCTS
64 C
65    17 DO 20 I = 1,NINC
66       DO 20 K = 1,NVEC
67    20 W(K,I)=DVMPY(N,A(IW(I,5),1),A(IW(I,5),2),V(1,K),V(2,K))
68
69       IF(IOUT .EQ. 1) THEN
70        WRITE(6,112) TEXT
71        DO 82 J = 1,NVEC,7
72        DO 83 I = 1,NINC
73    83  WRITE(6,'(1X,I9,7E15.6)') I,(W(J2,I),J2=J,MIN(NVEC,J+6))
74    82  WRITE(6,'(1X)')
75       END IF
76
77 C*****DETERMINES REDUNDANT INEQUALITIES AND CHOOSES NEW ONE
78
79       DO 40 K = 1,M
80       IW(K,3)=0
81       IF(IW(K,2) .EQ. 0) GO TO 40
82       DO 45 I = 1,NVEC
83       IF(DVMPY(N,A(K,1),A(K,2),V(1,I),V(2,I)) .GT. 0) IW(K,3)=IW(K,3)+1
84    45 CONTINUE
85    40 CONTINUE
86       NNEG=NVEC
87       DO 48 K = 1,M
88       IF(IW(K,2) .EQ. 0) GO TO 48
89       IF(IW(K,3) .EQ. 0) THEN
90        WRITE(ERRTXT,103) K
91        CALL MTLPRT(NAME,'F500.3',ERRTXT)
92        RETURN
93       END IF
94       IF(IW(K,3) .EQ. NVEC) THEN
95        IW(K,1)=K
96        IW(K,2)=0
97       END IF
98       IF(IW(K,3) .LT. NNEG) THEN
99        NNEG=IW(K,3)
100        KNEW=K
101       END IF
102    48 CONTINUE
103       IF(NNEG .EQ. NVEC) THEN
104        DO 74 I = 1,M
105        IW(I,1)=I
106        DO 75 J = 1,NVEC
107        IF(DVMPY(N,A(I,1),A(I,2),V(1,J),V(2,J)) .GE. EPS) GO TO 75
108        IW(I,1)=0
109        GO TO 74
110    75  CONTINUE
111    74  CONTINUE
112        RETURN
113       END IF
114
115       IF(IOUT .EQ. 1) WRITE(6,113) TEXT,KNEW
116       IW(KNEW,2)=0
117 C
118 C*****COMPUTES VECTOR OF SCALAR PRODUCTS
119 C
120       DO 50 I = 1,NVEC
121    50 W(I,M+1)=DVMPY(N,A(KNEW,1),A(KNEW,2),V(1,I),V(2,I))
122
123       IF(IOUT .EQ. 1) THEN
124        WRITE(6,114) TEXT
125        DO 84 J = 1,NVEC,7
126    84  WRITE(6,'(1X,I9,7E15.6)') J,(W(J2,M+1),J2=J,MIN(NVEC,J+6))
127        WRITE(6,'(1X)')
128       END IF
129 C
130 C*****DETERMINES BASIS VECTORS FOR NEW CONE
131 C
132       NTVE=NVEC
133       DO 60 I = 1,NVEC-1
134       DO 60 J = I+1,NVEC
135       IF(W(I,M+1)*W(J,M+1) .GT. 0) GO TO 60
136       NT=0
137       DO 62 L = 1,NINC
138       IW(L,4)=1
139       IF(ABS(W(I,L)) .GT. EPS .OR. ABS(W(J,L)) .GT. EPS) GO TO 62
140       NT=NT+1
141       IW(L,4)=0
142    62 CONTINUE
143       IF(NT .LT. N-2) GO TO 60
144       DO 63 K = 1,NVEC
145       IF(K .EQ. I .OR. K .EQ. J) GO TO 63
146       MT=0
147       DO 64 L = 1,NINC
148       IF(IW(L,4) .EQ. 0 .AND. ABS(W(K,L)) .LT. EPS) MT=MT+1
149    64 CONTINUE
150       IF(MT .EQ. N-2) GO TO 60
151    63 CONTINUE
152       NTVE=NTVE+1
153       IF(NTVE .GT. MAXV) THEN
154        CALL MTLPRT(NAME,'F500.1','MAXV TOO SMALL')
155        RETURN
156       END IF
157       DO 65 L = 1,N
158    65 V(L,NTVE)=ABS(W(J,M+1))*V(L,I)+ABS(W(I,M+1))*V(L,J)
159    60 CONTINUE
160       DO 66 I = 1,NTVE
161       S=DVMPY(N,V(1,I),V(2,I),V(1,I),V(2,I))
162    66 CALL DVSCL(N,1/SQRT(S),V(1,I),V(2,I),V(1,I),V(2,I))
163 C
164 C*****ELIMINATES VECTORS WITH NEGATIVE SCALAR PRODUCT
165 C
166       NNEW=0
167       DO 70 I = 1,NVEC
168       IF(W(I,M+1) .LT. 0) GO TO 70
169       NNEW=NNEW+1
170       CALL DVCPY(N,V(1,I),V(2,I),V(1,NNEW),V(2,NNEW))
171    70 CONTINUE
172       DO 71 I = NVEC+1,NTVE
173       NNEW=NNEW+1
174    71 CALL DVCPY(N,V(1,I),V(2,I),V(1,NNEW),V(2,NNEW))
175       CALL DMSET(N,NTVE-NNEW,0D0,V(1,NNEW+1),V(1,NNEW+2),V(2,NNEW+1))
176       NVEC=NNEW
177       NINC=NINC+1
178       IW(NINC,5)=KNEW
179
180       IF(IOUT .EQ. 1) THEN
181        WRITE(6,115) TEXT
182        DO 86 I = 1,NVEC,7
183        DO 87 J = 1,N
184    87  WRITE(6,'(1X,I9,7E15.6)') J,(V(J,I2),I2=I,MIN(NVEC,I+6))
185    86  WRITE(6,'(1X)')
186       END IF
187       GO TO 17
188   103 FORMAT('INEQUALITY ',I5,' IS INCONSISTENT')
189   111 FORMAT(7X,A27,'THE N BASIS VECTORS OF THE FIRST CONE'/)
190   112 FORMAT(7X,A27,'THE MATRIX OF SCALAR PRODUCTS'/)
191   113 FORMAT(7X,A27,'THE NEW INEQUALITY HAS INDEX',I5/)
192   114 FORMAT(7X,A27,'THE VECTOR OF SCALAR PRODUCTS'/)
193   115 FORMAT(7X,A27,'THE MATRIX OF BASIS VECTORS'/)
194       END
195 #endif