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