MINICERN/mathlib/gen/f/comlr.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:33  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE COMLR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
  11       INTEGER I,J,L,M,N,EN,LL,MM,NM,IGH,IM1,ITS,LOW,MP1,ENM1,IERR
  12       REAL HR(NM,N),HI(NM,N),WR(N),WI(N)
  13       REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,MACHEP
  14       COMPLEX X,Y,Z
  15       REAL T1(2),T2(2),T3(2)
  16       EQUIVALENCE (X,T1(1),XR),(T1(2),XI),(Y,T2(1),YR),(T2(2),YI),
  17      X            (Z,T3(1),ZZR),(T3(2),ZZI)
  18 #if defined(CERNLIB_CDC)
  19       MACHEP=2.**(-47)
  20 #endif
  21 #if !defined(CERNLIB_CDC)
  22       MACHEP=2.**(-23)
  23 #endif
  24       IERR = 0
  25 C*UL  180 DO 200 I = 1, N
  26       DO 200 I = 1, N
  27          IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
  28          WR(I) = HR(I,I)
  29          WI(I) = HI(I,I)
  30   200 CONTINUE
  31       EN = IGH
  32       TR = 0.0
  33       TI = 0.0
  34   220 IF (EN .LT. LOW) GO TO 1001
  35       ITS = 0
  36       ENM1 = EN - 1
  37   240 DO 260 LL = LOW, EN
  38          L = EN + LOW - LL
  39          IF (L .EQ. LOW) GO TO 300
  40          IF (ABS(HR(L,L-1)) + ABS(HI(L,L-1)) .LE.
  41      X      MACHEP * (ABS(HR(L-1,L-1)) + ABS(HI(L-1,L-1))
  42      X             + ABS(HR(L,L)) +ABS(HI(L,L)))) GO TO 300
  43   260 CONTINUE
  44   300 IF (L .EQ. EN) GO TO 660
  45       IF (ITS .EQ. 30) GO TO 1000
  46       IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
  47       SR = HR(EN,EN)
  48       SI = HI(EN,EN)
  49       XR = HR(ENM1,EN) * HR(EN,ENM1) - HI(ENM1,EN) * HI(EN,ENM1)
  50       XI = HR(ENM1,EN) * HI(EN,ENM1) + HI(ENM1,EN) * HR(EN,ENM1)
  51       IF (XR .EQ. 0.0 .AND. XI .EQ. 0.0) GO TO 340
  52       YR = (HR(ENM1,ENM1) - SR) / 2.0
  53       YI = (HI(ENM1,ENM1) - SI) / 2.0
  54       Z = SQRT(CMPLX(YR**2-YI**2+XR,2.0*YR*YI+XI))
  55       IF (YR * ZZR + YI * ZZI .LT. 0.0) Z = -Z
  56       X = X / (Y + Z)
  57       SR = SR - XR
  58       SI = SI - XI
  59       GO TO 340
  60   320 SR = ABS(HR(EN,ENM1)) + ABS(HR(ENM1,EN-2))
  61       SI = ABS(HI(EN,ENM1)) + ABS(HI(ENM1,EN-2))
  62   340 DO 360 I = LOW, EN
  63          HR(I,I) = HR(I,I) - SR
  64          HI(I,I) = HI(I,I) - SI
  65   360 CONTINUE
  66       TR = TR + SR
  67       TI = TI + SI
  68       ITS = ITS + 1
  69       XR = ABS(HR(ENM1,ENM1)) + ABS(HI(ENM1,ENM1))
  70       YR = ABS(HR(EN,ENM1)) + ABS(HI(EN,ENM1))
  71       ZZR = ABS(HR(EN,EN)) + ABS(HI(EN,EN))
  72       DO 380 MM = L, ENM1
  73          M = ENM1 + L - MM
  74          IF (M .EQ. L) GO TO 420
  75          YI = YR
  76          YR = ABS(HR(M,M-1)) + ABS(HI(M,M-1))
  77          XI = ZZR
  78          ZZR = XR
  79          XR = ABS(HR(M-1,M-1)) + ABS(HI(M-1,M-1))
  80          IF (YR .LE. MACHEP * ZZR / YI * (ZZR + XR + XI)) GO TO 420
  81   380 CONTINUE
  82   420 MP1 = M + 1
  83       DO 520 I = MP1, EN
  84          IM1 = I - 1
  85          XR = HR(IM1,IM1)
  86          XI = HI(IM1,IM1)
  87          YR = HR(I,IM1)
  88          YI = HI(I,IM1)
  89          IF (ABS(XR) + ABS(XI) .GE. ABS(YR) + ABS(YI)) GO TO 460
  90          DO 440 J = IM1, N
  91             ZZR = HR(IM1,J)
  92             HR(IM1,J) = HR(I,J)
  93             HR(I,J) = ZZR
  94             ZZI = HI(IM1,J)
  95             HI(IM1,J) = HI(I,J)
  96             HI(I,J) = ZZI
  97   440    CONTINUE
  98          Z = X / Y
  99          WR(I) = 1.0
 100          GO TO 480
 101   460    Z = Y / X
 102          WR(I) = -1.0
 103   480    HR(I,IM1) = ZZR
 104          HI(I,IM1) = ZZI
 105          DO 500 J = I, EN
 106             HR(I,J) = HR(I,J) - ZZR * HR(IM1,J) + ZZI * HI(IM1,J)
 107             HI(I,J) = HI(I,J) - ZZR * HI(IM1,J) - ZZI * HR(IM1,J)
 108   500   CONTINUE
 109   520 CONTINUE
 110       DO 640 J = MP1, EN
 111          XR = HR(J,J-1)
 112          XI = HI(J,J-1)
 113          HR(J,J-1) = 0.0
 114          HI(J,J-1) = 0.0
 115          IF (WR(J) .LE. 0.0) GO TO 580
 116          DO 540 I = L, J
 117             ZZR = HR(I,J-1)
 118             HR(I,J-1) = HR(I,J)
 119             HR(I,J) = ZZR
 120             ZZI = HI(I,J-1)
 121             HI(I,J-1) = HI(I,J)
 122             HI(I,J) = ZZI
 123   540    CONTINUE
 124   580    DO 600 I = L, J
 125             HR(I,J-1) = HR(I,J-1) + XR * HR(I,J) - XI * HI(I,J)
 126             HI(I,J-1) = HI(I,J-1) + XR * HI(I,J) + XI * HR(I,J)
 127   600    CONTINUE
 128   640 CONTINUE
 129       GO TO 240
 130   660 WR(EN) = HR(EN,EN) + TR
 131       WI(EN) = HI(EN,EN) + TI
 132       EN = ENM1
 133       GO TO 220
 134  1000 IERR = EN
 135  1001 RETURN
 136       END