TABLE OF CONTENTS


ABINIT/dzegdi [ Functions ]

[ Top ] [ Functions ]

NAME

  dzgedi

FUNCTION

  This routine is the clone of zgefa.F90 using real*8 a(2) instead of complex*16
  for the purpose of ABINIT

COPYRIGHT

  Copyright (C) 2014-2017 ABINIT group (XG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      berryphase,dfptnl_mv,qmatrix,relaxpol,smatrix,uderiv

CHILDREN

SOURCE

 31 #if defined HAVE_CONFIG_H
 32 #include "config.h"
 33 #endif
 34 
 35 #include "abi_common.h"
 36 
 37 subroutine dzgedi(a,lda,n,ipvt,det,work,job)
 38 
 39 
 40 !This section has been created automatically by the script Abilint (TD).
 41 !Do not modify the following lines by hand.
 42 #undef ABI_FUNC
 43 #define ABI_FUNC 'dzgedi'
 44 !End of the abilint section
 45 
 46       implicit none
 47 
 48       integer :: lda,n,ipvt(n),job
 49       real*8 :: a(2,lda,n),det(2,2),work(2,n)
 50 !
 51 !     zgedi computes the determinant and inverse of a matrix
 52 !     using the factors computed by zgeco or zgefa.
 53 !
 54 !     on entry
 55 !
 56 !        a       complex*16(lda, n)
 57 !                the output from zgeco or zgefa.
 58 !
 59 !        lda     integer
 60 !                the leading dimension of the array  a .
 61 !
 62 !        n       integer
 63 !                the order of the matrix  a .
 64 !
 65 !        ipvt    integer(n)
 66 !                the pivot vector from zgeco or zgefa.
 67 !
 68 !        work    complex*16(n)
 69 !                work vector.  contents destroyed.
 70 !
 71 !        job     integer
 72 !                = 11   both determinant and inverse.
 73 !                = 01   inverse only.
 74 !                = 10   determinant only.
 75 !
 76 !     on return
 77 !
 78 !        a       inverse of original matrix if requested.
 79 !                otherwise unchanged.
 80 !
 81 !        det     complex*16(2)
 82 !                determinant of original matrix if requested.
 83 !                otherwise not referenced.
 84 !                determinant = det(1) * 10.0**det(2)
 85 !                with  1.0 .le. cabs1(det(1)) .lt. 10.0
 86 !                or  det(1) .eq. 0.0 .
 87 !
 88 !     error condition
 89 !
 90 !        a division by zero will occur if the input factor contains
 91 !        a zero on the diagonal and the inverse is requested.
 92 !        it will not occur if the subroutines are called correctly
 93 !        and if zgeco has set rcond .gt. 0.0 or zgefa has set
 94 !        info .eq. 0 .
 95 !
 96 !     linpack. this version dated 08/14/78 .
 97 !     cleve moler, university of new mexico, argonne national lab.
 98 !
 99 !     subroutines and functions
100 !
101 !     internal variables
102 !
103       double precision :: r(2),rk(2),rkj(2)
104       double precision :: ten,rinv2,rabs
105       integer :: i,j,k,kb,kp1,l,nm1
106 !
107 !     compute determinant
108 !
109       if (job/10 .eq. 0) go to 70
110          det(1,1) = 1.0d0; det(2,1) = 0.0d0
111          det(1,2) = 0.0d0; det(2,2) = 0.0d0
112          ten = 10.0d0
113          do i = 1, n
114             if (ipvt(i) .ne. i) then
115                 det(1,1) = -det(1,1)
116                 det(2,1) = -det(2,1)
117             end if
118             r(1)=det(1,1); r(2)=det(2,1)
119             det(1,1) = r(1)*a(1,i,i)-r(2)*a(2,i,i)
120             det(2,1) = r(2)*a(1,i,i)+r(1)*a(2,i,i)
121 !        ...exit
122             rabs = abs(det(1,1))+abs(det(2,1))
123             if (rabs .eq. 0.0d0) go to 60
124    10       continue
125             rabs = abs(det(1,1))+abs(det(2,1))
126             if (rabs .ge. 1.0d0) go to 20
127                det(1,1) = ten*det(1,1); det(2,1) = ten*det(2,1)
128                det(1,2) = det(1,2) - 1.0d0
129             go to 10
130    20       continue
131    30       continue
132             rabs = abs(det(1,1))+abs(det(2,1))
133             if (rabs .lt. ten) go to 40
134                det(1,1) = det(1,1)/ten; det(2,1) = det(2,1)/ten
135                det(1,2) = det(1,2) + 1.0d0
136             go to 30
137    40       continue
138          end do
139    60    continue
140    70 continue
141 !
142 !     compute inverse(u)
143 !
144       if (mod(job,10) .eq. 0) go to 150
145          do 100 k = 1, n
146             !a(k,k) = (1.0d0,0.0d0)/a(k,k)
147             !t = -a(k,k)
148             !call zscal(k-1,t,a(1,k),1)
149             rinv2 = 1.d0/(a(1,k,k)**2+a(2,k,k)**2)
150             a(1,k,k) =  rinv2*a(1,k,k)
151             a(2,k,k) = -rinv2*a(2,k,k)
152             rk(1) = -a(1,k,k); rk(2) = -a(2,k,k)
153             do i=1,k-1
154                r(1)=a(1,i,k)
155                r(2)=a(2,i,k)
156                a(1,i,k)=rk(1)*r(1)-rk(2)*r(2)
157                a(2,i,k)=rk(1)*r(2)+rk(2)*r(1)
158             end do
159             kp1 = k + 1
160             if (n .lt. kp1) go to 90
161             do 80 j = kp1, n
162                !t = a(k,j)
163                !a(k,j) = (0.0d0,0.0d0)
164                !call zaxpy(k,t,a(1,k),1,a(1,j),1)
165                rkj(1) = a(1,k,j); rkj(2) = a(2,k,j)
166                a(1,k,j) = 0.d0; a(2,k,j) = 0.d0
167                do i=1,k
168                   a(1,i,j)=rkj(1)*a(1,i,k)-rkj(2)*a(2,i,k)+a(1,i,j)
169                   a(2,i,j)=rkj(2)*a(1,i,k)+rkj(1)*a(2,i,k)+a(2,i,j)
170                end do
171    80       continue
172    90       continue
173   100    continue
174   do i=1,n
175   end do
176 !
177 !        form inverse(u)*inverse(l)
178 !
179          nm1 = n - 1
180          if (nm1 .lt. 1) go to 140
181          do 130 kb = 1, nm1
182             k = n - kb
183             kp1 = k + 1
184             do 110 i = kp1, n
185                work(1,i) = a(1,i,k); work(2,i) = a(2,i,k)
186                a(1,i,k) = 0.0d0; a(2,i,k) = 0.d0
187   110       continue
188             do 120 j = kp1, n
189                r(1) = work(1,j); r(2) = work(2,j)
190                !call zaxpy(n,t,a(1,j),1,a(1,k),1)
191                do i=1,n
192                   a(1,i,k)=r(1)*a(1,i,j)-r(2)*a(2,i,j)+a(1,i,k)
193                   a(2,i,k)=r(2)*a(1,i,j)+r(1)*a(2,i,j)+a(2,i,k)
194                end do
195   120       continue
196             l = ipvt(k)
197             if (l .ne. k) then
198                !call zswap(n,a(1,k),1,a(1,l),1)
199                do i=1,n
200                   r(1) = a(1,i,k); r(2) = a(2,i,k)
201                   a(1,i,k) = a(1,i,l); a(2,i,k) = a(2,i,l)
202                   a(1,i,l) = r(1); a(2,i,l) = r(2)
203                end do
204             end if
205   130    continue
206   140    continue
207   150 continue
208 
209 end subroutine dzgedi