1      module m_linpack
2      public :: zgedi, zgefa
3
4      integer, parameter, private :: dp = kind(1.d0)  ! to make it consistent
5
6c     Updated (somewhat) by A. Garcia, Sep 2016
7c
8      CONTAINS
9
10      subroutine zgedi(a,lda,n,ipvt,det,work,job)
11      integer lda,n,ipvt(n),job
12      complex(dp) a(lda,n),det(2),work(n)
13c
14c     zgedi computes the determinant and inverse of a matrix
15c     using the factors computed by zgeco or zgefa.
16c
17c     on entry
18c
19c        a       complex(dp)(lda, n)
20c                the output from zgeco or zgefa.
21c
22c        lda     integer
23c                the leading dimension of the array  a .
24c
25c        n       integer
26c                the order of the matrix  a .
27c
28c        ipvt    integer(n)
29c                the pivot vector from zgeco or zgefa.
30c
31c        work    complex(dp)(n)
32c                work vector.  contents destroyed.
33c
34c        job     integer
35c                = 11   both determinant and inverse.
36c                = 01   inverse only.
37c                = 10   determinant only.
38c
39c     on return
40c
41c        a       inverse of original matrix if requested.
42c                otherwise unchanged.
43c
44c        det     complex(dp)(2)
45c                determinant of original matrix if requested.
46c                otherwise not referenced.
47c                determinant = det(1) * 10.0**det(2)
48c                with  1.0 .le. cabs1(det(1)) .lt. 10.0
49c                or  det(1) .eq. 0.0 .
50c
51c     error condition
52c
53c        a division by zero will occur if the input factor contains
54c        a zero on the diagonal and the inverse is requested.
55c        it will not occur if the subroutines are called correctly
56c        and if zgeco has set rcond .gt. 0.0 or zgefa has set
57c        info .eq. 0 .
58c
59c     linpack. this version dated 08/14/78 .
60c     cleve moler, university of new mexico, argonne national lab.
61c
62c     subroutines and functions
63c
64c     blas zaxpy,zscal,zswap
65c     fortran abs,cmplx,mod
66c
67c     internal variables
68c
69      complex(dp) t
70      real(dp) ten
71      integer i,j,k,kb,kp1,l,nm1
72c
73      complex(dp) zdum
74      real(dp) cabs1
75      real(dp) real,aimag
76      complex(dp) zdumr,zdumi
77      real(dp) :: local_real, local_aimag
78      local_real(zdumr) = zdumr
79      local_aimag(zdumi) = (0.0d0,-1.0d0)*zdumi
80      cabs1(zdum) = abs(local_real(zdum)) + abs(local_aimag(zdum))
81c
82c     compute determinant
83c
84      if (job/10 .eq. 0) go to 70
85         det(1) = (1.0d0,0.0d0)
86         det(2) = (0.0d0,0.0d0)
87         ten = 10.0d0
88         do 50 i = 1, n
89            if (ipvt(i) .ne. i) det(1) = -det(1)
90            det(1) = a(i,i)*det(1)
91c        ...exit
92            if (cabs1(det(1)) .eq. 0.0d0) go to 60
93   10       if (cabs1(det(1)) .ge. 1.0d0) go to 20
94               det(1) = dcmplx(ten,0.0d0)*det(1)
95               det(2) = det(2) - (1.0d0,0.0d0)
96            go to 10
97   20       continue
98   30       if (cabs1(det(1)) .lt. ten) go to 40
99               det(1) = det(1)/dcmplx(ten,0.0d0)
100               det(2) = det(2) + (1.0d0,0.0d0)
101            go to 30
102   40       continue
103   50    continue
104   60    continue
105   70 continue
106c
107c     compute inverse(u)
108c
109      if (mod(job,10) .eq. 0) go to 150
110         do 100 k = 1, n
111            a(k,k) = (1.0d0,0.0d0)/a(k,k)
112            t = -a(k,k)
113#ifdef OLD_CRAY
114            call cscal(k-1,t,a(1,k),1)
115#else
116            call zscal(k-1,t,a(1,k),1)
117#endif
118            kp1 = k + 1
119            if (n .lt. kp1) go to 90
120            do 80 j = kp1, n
121               t = a(k,j)
122               a(k,j) = (0.0d0,0.0d0)
123#ifdef OLD_CRAY
124               call caxpy(k,t,a(1,k),1,a(1,j),1)
125#else
126               call zaxpy(k,t,a(1,k),1,a(1,j),1)
127#endif
128   80       continue
129   90       continue
130  100    continue
131c
132c        form inverse(u)*inverse(l)
133c
134         nm1 = n - 1
135         if (nm1 .lt. 1) go to 140
136         do 130 kb = 1, nm1
137            k = n - kb
138            kp1 = k + 1
139            do 110 i = kp1, n
140               work(i) = a(i,k)
141               a(i,k) = (0.0d0,0.0d0)
142  110       continue
143            do 120 j = kp1, n
144               t = work(j)
145#ifdef OLD_CRAY
146               call caxpy(n,t,a(1,j),1,a(1,k),1)
147#else
148               call zaxpy(n,t,a(1,j),1,a(1,k),1)
149#endif
150  120       continue
151            l = ipvt(k)
152#ifdef OLD_CRAY
153            if (l .ne. k) call cswap(n,a(1,k),1,a(1,l),1)
154#else
155            if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1)
156#endif
157  130    continue
158  140    continue
159  150 continue
160      return
161      end subroutine zgedi
162
163      subroutine zgefa(a,lda,n,ipvt,info)
164      integer lda,n,ipvt(n),info
165      complex(dp) a(lda,n)
166c
167c     zgefa factors a complex(dp) matrix by gaussian elimination.
168c
169c     zgefa is usually called by zgeco, but it can be called
170c     directly with a saving in time if  rcond  is not needed.
171c     (time for zgeco) = (1 + 9/n)*(time for zgefa) .
172c
173c     on entry
174c
175c        a       complex(dp)(lda, n)
176c                the matrix to be factored.
177c
178c        lda     integer
179c                the leading dimension of the array  a .
180c
181c        n       integer
182c                the order of the matrix  a .
183c
184c     on return
185c
186c        a       an upper triangular matrix and the multipliers
187c                which were used to obtain it.
188c                the factorization can be written  a = l*u  where
189c                l  is a product of permutation and unit lower
190c                triangular matrices and  u  is upper triangular.
191c
192c        ipvt    integer(n)
193c                an integer vector of pivot indices.
194c
195c        info    integer
196c                = 0  normal value.
197c                = k  if  u(k,k) .eq. 0.0 .  this is not an error
198c                     condition for this subroutine, but it does
199c                     indicate that zgesl or zgedi will divide by zero
200c                     if called.  use  rcond  in zgeco for a reliable
201c                     indication of singularity.
202c
203c     linpack. this version dated 08/14/78 .
204c     cleve moler, university of new mexico, argonne national lab.
205c
206c     subroutines and functions
207c
208c     blas zaxpy,zscal,izamax_l
209c     fortran abs
210c
211c     internal variables
212c
213      complex(dp) t
214      integer j,k,kp1,l,nm1
215c
216c!    Statement functions...
217c!
218      complex(dp) zdum
219      real(dp) cabs1
220      real(dp) real,aimag
221      complex(dp) zdumr,zdumi
222      real(zdumr) = zdumr
223      aimag(zdumi) = (0.0d0,-1.0d0)*zdumi
224      cabs1(zdum) = abs(real(zdum)) + abs(aimag(zdum))
225c
226c     gaussian elimination with partial pivoting
227c
228      info = 0
229      nm1 = n - 1
230      if (nm1 .lt. 1) go to 70
231      do 60 k = 1, nm1
232         kp1 = k + 1
233c
234c        find l = pivot index
235c
236         l = izamax_l(n-k+1,a(k,k),1) + k - 1
237         ipvt(k) = l
238c
239c        zero pivot implies this column already triangularized
240c
241         if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
242c
243c           interchange if necessary
244c
245            if (l .eq. k) go to 10
246               t = a(l,k)
247               a(l,k) = a(k,k)
248               a(k,k) = t
249   10       continue
250c
251c           compute multipliers
252c
253            t = -(1.0d0,0.0d0)/a(k,k)
254#ifdef OLD_CRAY
255            call cscal(n-k,t,a(k+1,k),1)
256#else
257            call zscal(n-k,t,a(k+1,k),1)
258#endif
259c
260c           row elimination with column indexing
261c
262            do 30 j = kp1, n
263               t = a(l,j)
264               if (l .eq. k) go to 20
265                  a(l,j) = a(k,j)
266                  a(k,j) = t
267   20          continue
268#ifdef OLD_CRAY
269               call caxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
270#else
271               call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
272#endif
273   30       continue
274         go to 50
275   40    continue
276            info = k
277   50    continue
278   60 continue
279   70 continue
280      ipvt(n) = n
281      if (cabs1(a(n,n)) .eq. 0.0d0) info = n
282      return
283      end subroutine zgefa
284
285      integer function izamax_l(n,zx,incx)
286c
287c     finds the index of element having max. absolute value.
288c     jack dongarra, 3/11/78.
289c
290      integer n,incx
291      integer :: ix, i
292      complex(dp) zx(*)
293      real(dp) smax
294c
295      izamax_l = 1
296      if(n.le.1)return
297      if(incx.eq.1)go to 20
298c
299c        code for increments not equal to 1
300c
301      ix = 1
302      if(incx.lt.0)ix = (-n+1)*incx + 1
303      smax = sdcabs1(zx(1))
304      ix = ix + incx
305      do 10 i = 2,n
306         if(sdcabs1(zx(ix)).le.smax) go to 5
307         izamax_l = ix
308         smax = sdcabs1(zx(ix))
309    5    ix = ix + incx
310   10 continue
311      return
312c
313c        code for increments equal to 1
314c
315   20 smax = sdcabs1(zx(1))
316      do 30 i = 2,n
317         if(sdcabs1(zx(i)).le.smax) go to 30
318         izamax_l = i
319         smax = sdcabs1(zx(i))
320   30 continue
321      return
322      end function izamax_l
323C
324C  The code below is a modified form of the routine sdcabs1
325C  changed to avoid problems on the Cray when compiled
326C  with optimisation turned on. JDG, November 1999
327C
328      real(dp) function sdcabs1(z)
329      complex(dp) z
330      real(dp) t(2)
331      t(1) = real(z)
332      t(2) = aimag(z)
333      sdcabs1 = abs(t(1)) + abs(t(2))
334      return
335      end function sdcabs1
336
337      end module m_linpack
338
339