← Back to branch summary

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

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/casyi.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/casyi.f
 
1
*DECK CASYI
 
2
      SUBROUTINE CASYI (Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  CASYI
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  Subsidiary to CBESI and CBESK
 
6
C***LIBRARY   SLATEC
 
7
C***TYPE      ALL (CASYI-A, ZASYI-A)
 
8
C***AUTHOR  Amos, D. E., (SNL)
 
9
C***DESCRIPTION
 
10
C
 
11
C     CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
 
12
C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE
 
13
C     REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
 
14
C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
 
15
C
 
16
C***SEE ALSO  CBESI, CBESK
 
17
C***ROUTINES CALLED  R1MACH
 
18
C***REVISION HISTORY  (YYMMDD)
 
19
C   830501  DATE WRITTEN
 
20
C   910415  Prologue converted to Version 4.0 format.  (BAB)
 
21
C***END PROLOGUE  CASYI
 
22
      COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2,
 
23
     * Y, Z
 
24
      REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU,
 
25
     * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X,
 
26
     * YY, R1MACH
 
27
      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
 
28
      DIMENSION Y(N)
 
29
      DATA PI, RTPI  /3.14159265358979324E0 , 0.159154943091895336E0 /
 
30
      DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
 
31
C***FIRST EXECUTABLE STATEMENT  CASYI
 
32
      NZ = 0
 
33
      AZ = ABS(Z)
 
34
      X = REAL(Z)
 
35
      ARM = 1.0E+3*R1MACH(1)
 
36
      RTR1 = SQRT(ARM)
 
37
      IL = MIN(2,N)
 
38
      DFNU = FNU + (N-IL)
 
39
C-----------------------------------------------------------------------
 
40
C     OVERFLOW TEST
 
41
C-----------------------------------------------------------------------
 
42
      AK1 = CMPLX(RTPI,0.0E0)/Z
 
43
      AK1 = CSQRT(AK1)
 
44
      CZ = Z
 
45
      IF (KODE.EQ.2) CZ = Z - CMPLX(X,0.0E0)
 
46
      ACZ = REAL(CZ)
 
47
      IF (ABS(ACZ).GT.ELIM) GO TO 80
 
48
      DNU2 = DFNU + DFNU
 
49
      KODED = 1
 
50
      IF ((ABS(ACZ).GT.ALIM) .AND. (N.GT.2)) GO TO 10
 
51
      KODED = 0
 
52
      AK1 = AK1*CEXP(CZ)
 
53
   10 CONTINUE
 
54
      FDN = 0.0E0
 
55
      IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
 
56
      EZ = Z*CMPLX(8.0E0,0.0E0)
 
57
C-----------------------------------------------------------------------
 
58
C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
 
59
C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
 
60
C     EXPANSION FOR THE IMAGINARY PART.
 
61
C-----------------------------------------------------------------------
 
62
      AEZ = 8.0E0*AZ
 
63
      S = TOL/AEZ
 
64
      JL = RL+RL + 2
 
65
      YY = AIMAG(Z)
 
66
      P1 = CZERO
 
67
      IF (YY.EQ.0.0E0) GO TO 20
 
68
C-----------------------------------------------------------------------
 
69
C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
 
70
C     SIGNIFICANCE WHEN FNU OR N IS LARGE
 
71
C-----------------------------------------------------------------------
 
72
      INU = FNU
 
73
      ARG = (FNU-INU)*PI
 
74
      INU = INU + N - IL
 
75
      AK = -SIN(ARG)
 
76
      BK = COS(ARG)
 
77
      IF (YY.LT.0.0E0) BK = -BK
 
78
      P1 = CMPLX(AK,BK)
 
79
      IF (MOD(INU,2).EQ.1) P1 = -P1
 
80
   20 CONTINUE
 
81
      DO 50 K=1,IL
 
82
        SQK = FDN - 1.0E0
 
83
        ATOL = S*ABS(SQK)
 
84
        SGN = 1.0E0
 
85
        CS1 = CONE
 
86
        CS2 = CONE
 
87
        CK = CONE
 
88
        AK = 0.0E0
 
89
        AA = 1.0E0
 
90
        BB = AEZ
 
91
        DK = EZ
 
92
        DO 30 J=1,JL
 
93
          CK = CK*CMPLX(SQK,0.0E0)/DK
 
94
          CS2 = CS2 + CK
 
95
          SGN = -SGN
 
96
          CS1 = CS1 + CK*CMPLX(SGN,0.0E0)
 
97
          DK = DK + EZ
 
98
          AA = AA*ABS(SQK)/BB
 
99
          BB = BB + AEZ
 
100
          AK = AK + 8.0E0
 
101
          SQK = SQK - AK
 
102
          IF (AA.LE.ATOL) GO TO 40
 
103
   30   CONTINUE
 
104
        GO TO 90
 
105
   40   CONTINUE
 
106
        S2 = CS1
 
107
        IF (X+X.LT.ELIM) S2 = S2 + P1*CS2*CEXP(-Z-Z)
 
108
        FDN = FDN + 8.0E0*DFNU + 4.0E0
 
109
        P1 = -P1
 
110
        M = N - IL + K
 
111
        Y(M) = S2*AK1
 
112
   50 CONTINUE
 
113
      IF (N.LE.2) RETURN
 
114
      NN = N
 
115
      K = NN - 2
 
116
      AK = K
 
117
      RZ = (CONE+CONE)/Z
 
118
      IB = 3
 
119
      DO 60 I=IB,NN
 
120
        Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2)
 
121
        AK = AK - 1.0E0
 
122
        K = K - 1
 
123
   60 CONTINUE
 
124
      IF (KODED.EQ.0) RETURN
 
125
      CK = CEXP(CZ)
 
126
      DO 70 I=1,NN
 
127
        Y(I) = Y(I)*CK
 
128
   70 CONTINUE
 
129
      RETURN
 
130
   80 CONTINUE
 
131
      NZ = -1
 
132
      RETURN
 
133
   90 CONTINUE
 
134
      NZ=-2
 
135
      RETURN
 
136
      END