← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/linalg/src/lu.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/linalg/src/lu.f
 
1
 
 
2
c     Calculate LU decomposition of a matrix
 
3
c     Author: Pearu Peterson, March 2002
 
4
c
 
5
c     prefixes: d,z,s,c   (double,complex double,float,complex float)
 
6
c     suffixes: _c,_r     (column major order,row major order)
 
7
 
 
8
      subroutine dlu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
 
9
      integer m,n,piv(k),i,j,k,permute_l,m1
 
10
      double precision l(m,k),u(k,n),a(m,n)
 
11
      double precision p(m1,m1)
 
12
 
 
13
cf2py intent(in,copy) :: a
 
14
cf2py intent(out) :: info
 
15
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
 
16
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
 
17
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
 
18
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
 
19
cf2py intent(out) :: p,l,u
 
20
cf2py integer optional,intent(in):: permute_l = 0
 
21
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
 
22
cf2py depend(m1) :: p
 
23
 
 
24
cf2py callprotoargument double*,double*,double*,double*,int*,int*,int*,int*,int*,int*,int*
 
25
 
 
26
      external dgetrf,dlaswp
 
27
      call dgetrf(m,n,a,m,piv,info)
 
28
      if (info.lt.0) then
 
29
         return
 
30
      endif
 
31
      do 20 i=1,m
 
32
         do 10, j=1,n
 
33
            if (j.le.k) then
 
34
               if (i.eq.j) then
 
35
                  l(i,j) = 1d0
 
36
               elseif (i.gt.j) then
 
37
                  l(i,j) = a(i,j)
 
38
               endif
 
39
            endif
 
40
            if (i.le.k.and.i.le.j) then
 
41
               u(i,j) = a(i,j)
 
42
            endif
 
43
 10      continue
 
44
 20   continue
 
45
      if (permute_l.ne.0) then
 
46
         call dlaswp(n,l,m,1,k,piv,1)
 
47
      else
 
48
         do 25 i=1,m
 
49
            p(i,i)=1d0
 
50
 25       continue
 
51
         call dlaswp(m,p,m,1,k,piv,1)
 
52
      endif
 
53
      end
 
54
 
 
55
      subroutine zlu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
 
56
      integer m,n,piv(k),i,j,k,permute_l,m1
 
57
      complex*16 l(m,k),u(k,n),a(m,n)
 
58
      double precision p(m1,m1)
 
59
 
 
60
cf2py intent(in,copy) :: a
 
61
cf2py intent(out) :: info
 
62
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
 
63
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
 
64
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
 
65
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
 
66
cf2py intent(out) :: p,l,u
 
67
cf2py integer optional,intent(in):: permute_l = 0
 
68
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
 
69
cf2py depend(m1) :: p
 
70
 
 
71
cf2py callprotoargument double*,complex_double*,complex_double*,complex_double*,int*,int*,int*,int*,int*,int*,int*
 
72
 
 
73
      external zgetrf,zlaswp,dlaswp
 
74
      call zgetrf(m,n,a,m,piv,info)
 
75
      if (info.lt.0) then
 
76
         return
 
77
      endif
 
78
      do 20 i=1,m
 
79
         do 10, j=1,n
 
80
            if (j.le.k) then
 
81
               if (i.eq.j) then
 
82
                  l(i,j) = 1d0
 
83
               elseif (i.gt.j) then
 
84
                  l(i,j) = a(i,j)
 
85
               endif
 
86
            endif
 
87
            if (i.le.k.and.i.le.j) then
 
88
               u(i,j) = a(i,j)
 
89
            endif
 
90
 10      continue
 
91
 20   continue
 
92
      if (permute_l.ne.0) then
 
93
         call zlaswp(n,l,m,1,k,piv,1)
 
94
      else
 
95
         do 25 i=1,m
 
96
            p(i,i)=1d0
 
97
 25       continue
 
98
         call dlaswp(m,p,m,1,k,piv,1)
 
99
      endif
 
100
      end
 
101
 
 
102
      subroutine slu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
 
103
      integer m,n,piv(k),i,j,k,permute_l,m1
 
104
      real l(m,k),u(k,n),a(m,n)
 
105
      real p(m1,m1)
 
106
 
 
107
cf2py intent(in,copy) :: a
 
108
cf2py intent(out) :: info
 
109
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
 
110
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
 
111
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
 
112
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
 
113
cf2py intent(out) :: p,l,u
 
114
cf2py integer optional,intent(in):: permute_l = 0
 
115
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
 
116
cf2py depend(m1) :: p
 
117
 
 
118
cf2py callprotoargument float*,float*,float*,float*,int*,int*,int*,int*,int*,int*,int*
 
119
 
 
120
      external sgetrf,slaswp
 
121
      call sgetrf(m,n,a,m,piv,info)
 
122
      if (info.lt.0) then
 
123
         return
 
124
      endif
 
125
      do 20 i=1,m
 
126
         do 10, j=1,n
 
127
            if (j.le.k) then
 
128
               if (i.eq.j) then
 
129
                  l(i,j) = 1e0
 
130
               elseif (i.gt.j) then
 
131
                  l(i,j) = a(i,j)
 
132
               endif
 
133
            endif
 
134
            if (i.le.k.and.i.le.j) then
 
135
               u(i,j) = a(i,j)
 
136
            endif
 
137
 10      continue
 
138
 20   continue
 
139
      if (permute_l.ne.0) then
 
140
         call slaswp(n,l,m,1,k,piv,1)
 
141
      else
 
142
         do 25 i=1,m
 
143
            p(i,i)=1e0
 
144
 25       continue
 
145
         call slaswp(m,p,m,1,k,piv,1)
 
146
      endif
 
147
      end
 
148
 
 
149
      subroutine clu_c(p,l,u,a,m,n,k,piv,info,permute_l,m1)
 
150
      integer m,n,piv(k),i,j,k,permute_l,m1
 
151
      complex l(m,k),u(k,n),a(m,n)
 
152
      real p(m1,m1)
 
153
 
 
154
cf2py intent(in,copy) :: a
 
155
cf2py intent(out) :: info
 
156
cf2py integer intent(hide,cache),depend(k),dimension(k) :: piv
 
157
cf2py integer intent(hide),depend(a) :: m = shape(a,0)
 
158
cf2py integer intent(hide),depend(a) :: n = shape(a,1)
 
159
cf2py integer intent(hide),depend(m,n) :: k = (m<n?m:n)
 
160
cf2py intent(out) :: p,l,u
 
161
cf2py integer optional,intent(in):: permute_l = 0
 
162
cf2py integer intent(hide),depend(permute_l,m) :: m1 = (permute_l?1:m)
 
163
cf2py depend(m1) :: p
 
164
 
 
165
cf2py callprotoargument float*,complex_float*,complex_float*,complex_float*,int*,int*,int*,int*,int*,int*,int*
 
166
 
 
167
      external cgetrf,claswp,slaswp
 
168
      call cgetrf(m,n,a,m,piv,info)
 
169
      if (info.lt.0) then
 
170
         return
 
171
      endif
 
172
      do 20 i=1,m
 
173
         do 10, j=1,n
 
174
            if (j.le.k) then
 
175
               if (i.eq.j) then
 
176
                  l(i,j) = 1e0
 
177
               elseif (i.gt.j) then
 
178
                  l(i,j) = a(i,j)
 
179
               endif
 
180
            endif
 
181
            if (i.le.k.and.i.le.j) then
 
182
               u(i,j) = a(i,j)
 
183
            endif
 
184
 10      continue
 
185
 20   continue
 
186
      if (permute_l.ne.0) then
 
187
         call claswp(n,l,m,1,k,piv,1)
 
188
      else
 
189
         do 25 i=1,m
 
190
            p(i,i)=1e0
 
191
 25       continue
 
192
         call slaswp(m,p,m,1,k,piv,1)
 
193
      endif
 
194
      end