← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/integrate/linpack_lite/dgefa.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:
Lib/integrate/linpack_lite/dgefa.f
1
 
      subroutine dgefa(a,lda,n,ipvt,info)
2
 
      integer lda,n,ipvt(1),info
3
 
      double precision a(lda,1)
4
 
c
5
 
c     dgefa factors a double precision matrix by gaussian elimination.
6
 
c
7
 
c     dgefa is usually called by dgeco, but it can be called
8
 
c     directly with a saving in time if  rcond  is not needed.
9
 
c     (time for dgeco) = (1 + 9/n)*(time for dgefa) .
10
 
c
11
 
c     on entry
12
 
c
13
 
c        a       double precision(lda, n)
14
 
c                the matrix to be factored.
15
 
c
16
 
c        lda     integer
17
 
c                the leading dimension of the array  a .
18
 
c
19
 
c        n       integer
20
 
c                the order of the matrix  a .
21
 
c
22
 
c     on return
23
 
c
24
 
c        a       an upper triangular matrix and the multipliers
25
 
c                which were used to obtain it.
26
 
c                the factorization can be written  a = l*u  where
27
 
c                l  is a product of permutation and unit lower
28
 
c                triangular matrices and  u  is upper triangular.
29
 
c
30
 
c        ipvt    integer(n)
31
 
c                an integer vector of pivot indices.
32
 
c
33
 
c        info    integer
34
 
c                = 0  normal value.
35
 
c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
36
 
c                     condition for this subroutine, but it does
37
 
c                     indicate that dgesl or dgedi will divide by zero
38
 
c                     if called.  use  rcond  in dgeco for a reliable
39
 
c                     indication of singularity.
40
 
c
41
 
c     linpack. this version dated 08/14/78 .
42
 
c     cleve moler, university of new mexico, argonne national lab.
43
 
c
44
 
c     subroutines and functions
45
 
c
46
 
c     blas daxpy,dscal,idamax
47
 
c
48
 
c     internal variables
49
 
c
50
 
      double precision t
51
 
      integer idamax,j,k,kp1,l,nm1
52
 
c
53
 
c
54
 
c     gaussian elimination with partial pivoting
55
 
c
56
 
      info = 0
57
 
      nm1 = n - 1
58
 
      if (nm1 .lt. 1) go to 70
59
 
      do 60 k = 1, nm1
60
 
         kp1 = k + 1
61
 
c
62
 
c        find l = pivot index
63
 
c
64
 
         l = idamax(n-k+1,a(k,k),1) + k - 1
65
 
         ipvt(k) = l
66
 
c
67
 
c        zero pivot implies this column already triangularized
68
 
c
69
 
         if (a(l,k) .eq. 0.0d0) go to 40
70
 
c
71
 
c           interchange if necessary
72
 
c
73
 
            if (l .eq. k) go to 10
74
 
               t = a(l,k)
75
 
               a(l,k) = a(k,k)
76
 
               a(k,k) = t
77
 
   10       continue
78
 
c
79
 
c           compute multipliers
80
 
c
81
 
            t = -1.0d0/a(k,k)
82
 
            call dscal(n-k,t,a(k+1,k),1)
83
 
c
84
 
c           row elimination with column indexing
85
 
c
86
 
            do 30 j = kp1, n
87
 
               t = a(l,j)
88
 
               if (l .eq. k) go to 20
89
 
                  a(l,j) = a(k,j)
90
 
                  a(k,j) = t
91
 
   20          continue
92
 
               call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
93
 
   30       continue
94
 
         go to 50
95
 
   40    continue
96
 
            info = k
97
 
   50    continue
98
 
   60 continue
99
 
   70 continue
100
 
      ipvt(n) = n
101
 
      if (a(n,n) .eq. 0.0d0) info = n
102
 
      return
103
 
      end