← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/integrate/linpack_lite/dgtsl.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/integrate/linpack_lite/dgtsl.f
 
1
      subroutine dgtsl(n,c,d,e,b,info)
 
2
      integer n,info
 
3
      double precision c(1),d(1),e(1),b(1)
 
4
c
 
5
c     dgtsl given a general tridiagonal matrix and a right hand
 
6
c     side will find the solution.
 
7
c
 
8
c     on entry
 
9
c
 
10
c        n       integer
 
11
c                is the order of the tridiagonal matrix.
 
12
c
 
13
c        c       double precision(n)
 
14
c                is the subdiagonal of the tridiagonal matrix.
 
15
c                c(2) through c(n) should contain the subdiagonal.
 
16
c                on output c is destroyed.
 
17
c
 
18
c        d       double precision(n)
 
19
c                is the diagonal of the tridiagonal matrix.
 
20
c                on output d is destroyed.
 
21
c
 
22
c        e       double precision(n)
 
23
c                is the superdiagonal of the tridiagonal matrix.
 
24
c                e(1) through e(n-1) should contain the superdiagonal.
 
25
c                on output e is destroyed.
 
26
c
 
27
c        b       double precision(n)
 
28
c                is the right hand side vector.
 
29
c
 
30
c     on return
 
31
c
 
32
c        b       is the solution vector.
 
33
c
 
34
c        info    integer
 
35
c                = 0 normal value.
 
36
c                = k if the k-th element of the diagonal becomes
 
37
c                    exactly zero.  the subroutine returns when
 
38
c                    this is detected.
 
39
c
 
40
c     linpack. this version dated 08/14/78 .
 
41
c     jack dongarra, argonne national laboratory.
 
42
c
 
43
c     no externals
 
44
c     fortran dabs
 
45
c
 
46
c     internal variables
 
47
c
 
48
      integer k,kb,kp1,nm1,nm2
 
49
      double precision t
 
50
c     begin block permitting ...exits to 100
 
51
c
 
52
         info = 0
 
53
         c(1) = d(1)
 
54
         nm1 = n - 1
 
55
         if (nm1 .lt. 1) go to 40
 
56
            d(1) = e(1)
 
57
            e(1) = 0.0d0
 
58
            e(n) = 0.0d0
 
59
c
 
60
            do 30 k = 1, nm1
 
61
               kp1 = k + 1
 
62
c
 
63
c              find the largest of the two rows
 
64
c
 
65
               if (dabs(c(kp1)) .lt. dabs(c(k))) go to 10
 
66
c
 
67
c                 interchange row
 
68
c
 
69
                  t = c(kp1)
 
70
                  c(kp1) = c(k)
 
71
                  c(k) = t
 
72
                  t = d(kp1)
 
73
                  d(kp1) = d(k)
 
74
                  d(k) = t
 
75
                  t = e(kp1)
 
76
                  e(kp1) = e(k)
 
77
                  e(k) = t
 
78
                  t = b(kp1)
 
79
                  b(kp1) = b(k)
 
80
                  b(k) = t
 
81
   10          continue
 
82
c
 
83
c              zero elements
 
84
c
 
85
               if (c(k) .ne. 0.0d0) go to 20
 
86
                  info = k
 
87
c     ............exit
 
88
                  go to 100
 
89
   20          continue
 
90
               t = -c(kp1)/c(k)
 
91
               c(kp1) = d(kp1) + t*d(k)
 
92
               d(kp1) = e(kp1) + t*e(k)
 
93
               e(kp1) = 0.0d0
 
94
               b(kp1) = b(kp1) + t*b(k)
 
95
   30       continue
 
96
   40    continue
 
97
         if (c(n) .ne. 0.0d0) go to 50
 
98
            info = n
 
99
         go to 90
 
100
   50    continue
 
101
c
 
102
c           back solve
 
103
c
 
104
            nm2 = n - 2
 
105
            b(n) = b(n)/c(n)
 
106
            if (n .eq. 1) go to 80
 
107
               b(nm1) = (b(nm1) - d(nm1)*b(n))/c(nm1)
 
108
               if (nm2 .lt. 1) go to 70
 
109
               do 60 kb = 1, nm2
 
110
                  k = nm2 - kb + 1
 
111
                  b(k) = (b(k) - d(k)*b(k+1) - e(k)*b(k+2))/c(k)
 
112
   60          continue
 
113
   70          continue
 
114
   80       continue
 
115
   90    continue
 
116
  100 continue
 
117
c
 
118
      return
 
119
      end