← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/integrate/linpack_lite/dgesl.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/integrate/linpack_lite/dgesl.f
 
1
      subroutine dgesl(a,lda,n,ipvt,b,job)
 
2
      integer lda,n,ipvt(1),job
 
3
      double precision a(lda,1),b(1)
 
4
c
 
5
c     dgesl solves the double precision system
 
6
c     a * x = b  or  trans(a) * x = b
 
7
c     using the factors computed by dgeco or dgefa.
 
8
c
 
9
c     on entry
 
10
c
 
11
c        a       double precision(lda, n)
 
12
c                the output from dgeco or dgefa.
 
13
c
 
14
c        lda     integer
 
15
c                the leading dimension of the array  a .
 
16
c
 
17
c        n       integer
 
18
c                the order of the matrix  a .
 
19
c
 
20
c        ipvt    integer(n)
 
21
c                the pivot vector from dgeco or dgefa.
 
22
c
 
23
c        b       double precision(n)
 
24
c                the right hand side vector.
 
25
c
 
26
c        job     integer
 
27
c                = 0         to solve  a*x = b ,
 
28
c                = nonzero   to solve  trans(a)*x = b  where
 
29
c                            trans(a)  is the transpose.
 
30
c
 
31
c     on return
 
32
c
 
33
c        b       the solution vector  x .
 
34
c
 
35
c     error condition
 
36
c
 
37
c        a division by zero will occur if the input factor contains a
 
38
c        zero on the diagonal.  technically this indicates singularity
 
39
c        but it is often caused by improper arguments or improper
 
40
c        setting of lda .  it will not occur if the subroutines are
 
41
c        called correctly and if dgeco has set rcond .gt. 0.0
 
42
c        or dgefa has set info .eq. 0 .
 
43
c
 
44
c     to compute  inverse(a) * c  where  c  is a matrix
 
45
c     with  p  columns
 
46
c           call dgeco(a,lda,n,ipvt,rcond,z)
 
47
c           if (rcond is too small) go to ...
 
48
c           do 10 j = 1, p
 
49
c              call dgesl(a,lda,n,ipvt,c(1,j),0)
 
50
c        10 continue
 
51
c
 
52
c     linpack. this version dated 08/14/78 .
 
53
c     cleve moler, university of new mexico, argonne national lab.
 
54
c
 
55
c     subroutines and functions
 
56
c
 
57
c     blas daxpy,ddot
 
58
c
 
59
c     internal variables
 
60
c
 
61
      double precision ddot,t
 
62
      integer k,kb,l,nm1
 
63
c
 
64
      nm1 = n - 1
 
65
      if (job .ne. 0) go to 50
 
66
c
 
67
c        job = 0 , solve  a * x = b
 
68
c        first solve  l*y = b
 
69
c
 
70
         if (nm1 .lt. 1) go to 30
 
71
         do 20 k = 1, nm1
 
72
            l = ipvt(k)
 
73
            t = b(l)
 
74
            if (l .eq. k) go to 10
 
75
               b(l) = b(k)
 
76
               b(k) = t
 
77
   10       continue
 
78
            call daxpy(n-k,t,a(k+1,k),1,b(k+1),1)
 
79
   20    continue
 
80
   30    continue
 
81
c
 
82
c        now solve  u*x = y
 
83
c
 
84
         do 40 kb = 1, n
 
85
            k = n + 1 - kb
 
86
            b(k) = b(k)/a(k,k)
 
87
            t = -b(k)
 
88
            call daxpy(k-1,t,a(1,k),1,b(1),1)
 
89
   40    continue
 
90
      go to 100
 
91
   50 continue
 
92
c
 
93
c        job = nonzero, solve  trans(a) * x = b
 
94
c        first solve  trans(u)*y = b
 
95
c
 
96
         do 60 k = 1, n
 
97
            t = ddot(k-1,a(1,k),1,b(1),1)
 
98
            b(k) = (b(k) - t)/a(k,k)
 
99
   60    continue
 
100
c
 
101
c        now solve trans(l)*x = y
 
102
c
 
103
         if (nm1 .lt. 1) go to 90
 
104
         do 80 kb = 1, nm1
 
105
            k = n - kb
 
106
            b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1)
 
107
            l = ipvt(k)
 
108
            if (l .eq. k) go to 70
 
109
               t = b(l)
 
110
               b(l) = b(k)
 
111
               b(k) = t
 
112
   70       continue
 
113
   80    continue
 
114
   90    continue
 
115
  100 continue
 
116
      return
 
117
      end