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