TABLE OF CONTENTS
ABINIT/dzegdi [ 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