← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/cwrsk.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/cwrsk.f
 
1
*DECK CWRSK
 
2
      SUBROUTINE CWRSK (ZR, FNU, KODE, N, Y, NZ, CW, TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  CWRSK
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to CBESI and CBESK
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      ALL (CWRSK-A, ZWRSK-A)
 
8
C***AUTHOR  Amos, D. E., (SNL)
 
9
C***DESCRIPTION
 
10
C
 
11
C     CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
 
12
C     NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
 
13
C
 
14
C***SEE ALSO  CBESI, CBESK
 
15
C***ROUTINES CALLED  CBKNU, CRATI, R1MACH
 
16
C***REVISION HISTORY  (YYMMDD)
 
17
C   830501  DATE WRITTEN
 
18
C   910415  Prologue converted to Version 4.0 format.  (BAB)
 
19
C***END PROLOGUE  CWRSK
 
20
      COMPLEX CINU, CSCL, CT, CW, C1, C2, RCT, ST, Y, ZR
 
21
      REAL ACT, ACW, ALIM, ASCLE, ELIM, FNU, S1, S2, TOL, YY, R1MACH
 
22
      INTEGER I, KODE, N, NW, NZ
 
23
      DIMENSION Y(N), CW(2)
 
24
C***FIRST EXECUTABLE STATEMENT  CWRSK
 
25
C-----------------------------------------------------------------------
 
26
C     I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
 
27
C     Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
 
28
C     WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
 
29
C-----------------------------------------------------------------------
 
30
      NZ = 0
 
31
      CALL CBKNU(ZR, FNU, KODE, 2, CW, NW, TOL, ELIM, ALIM)
 
32
      IF (NW.NE.0) GO TO 50
 
33
      CALL CRATI(ZR, FNU, N, Y, TOL)
 
34
C-----------------------------------------------------------------------
 
35
C     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
 
36
C     R(FNU+J-1,Z)=Y(J),  J=1,...,N
 
37
C-----------------------------------------------------------------------
 
38
      CINU = CMPLX(1.0E0,0.0E0)
 
39
      IF (KODE.EQ.1) GO TO 10
 
40
      YY = AIMAG(ZR)
 
41
      S1 = COS(YY)
 
42
      S2 = SIN(YY)
 
43
      CINU = CMPLX(S1,S2)
 
44
   10 CONTINUE
 
45
C-----------------------------------------------------------------------
 
46
C     ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
 
47
C     THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
 
48
C     SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
 
49
C     THE RESULT IS ON SCALE.
 
50
C-----------------------------------------------------------------------
 
51
      ACW = ABS(CW(2))
 
52
      ASCLE = 1.0E+3*R1MACH(1)/TOL
 
53
      CSCL = CMPLX(1.0E0,0.0E0)
 
54
      IF (ACW.GT.ASCLE) GO TO 20
 
55
      CSCL = CMPLX(1.0E0/TOL,0.0E0)
 
56
      GO TO 30
 
57
   20 CONTINUE
 
58
      ASCLE = 1.0E0/ASCLE
 
59
      IF (ACW.LT.ASCLE) GO TO 30
 
60
      CSCL = CMPLX(TOL,0.0E0)
 
61
   30 CONTINUE
 
62
      C1 = CW(1)*CSCL
 
63
      C2 = CW(2)*CSCL
 
64
      ST = Y(1)
 
65
C-----------------------------------------------------------------------
 
66
C     CINU=CINU*(CONJG(CT)/ABS(CT))*(1.0E0/ABS(CT) PREVENTS
 
67
C     UNDER- OR OVERFLOW PREMATURELY BY SQUARING ABS(CT)
 
68
C-----------------------------------------------------------------------
 
69
      CT = ZR*(C2+ST*C1)
 
70
      ACT = ABS(CT)
 
71
      RCT = CMPLX(1.0E0/ACT,0.0E0)
 
72
      CT = CONJG(CT)*RCT
 
73
      CINU = CINU*RCT*CT
 
74
      Y(1) = CINU*CSCL
 
75
      IF (N.EQ.1) RETURN
 
76
      DO 40 I=2,N
 
77
        CINU = ST*CINU
 
78
        ST = Y(I)
 
79
        Y(I) = CINU*CSCL
 
80
   40 CONTINUE
 
81
      RETURN
 
82
   50 CONTINUE
 
83
      NZ = -1
 
84
      IF(NW.EQ.(-2)) NZ=-2
 
85
      RETURN
 
86
      END