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