Permalink
Find file
Fetching contributors…
Cannot retrieve contributors at this time
454 lines (447 sloc) 13 KB
SUBROUTINE inv1(A,Ainv,N)
IMPLICIT NONE
INTEGER N,i,j,info,ipvt(N),job
DOUBLE COMPLEX A(N,N),Ainv(N,N)
DOUBLE COMPLEX det(2),work(N)
******************
**** The determinant is det(1)*10.d0**det(2)
****
****job: 11 both determinant and inverse
**** 01 inverse only
*** 10 determinant only
****
*****************
job=01
DO i=1,N
DO j=1,N
Ainv(i,j)=A(i,j)
ENDDO
ENDDO
CALL zgefa(Ainv,N,N,ipvt,info)
CALL zgedi(Ainv,N,N,ipvt,det,work,job)
DO i=1,N
DO j=1,N
IF (ABS(DBLE(Ainv(i,j))).LT.1.d-16) THEN
Ainv(i,j)=-(0.d0,1.d0)*DBLE((0.d0,1.d0)*Ainv(i,j))
ENDIF
IF (ABS(DBLE((0.d0,1.d0)*Ainv(i,j))).LT.1.d-16) THEN
Ainv(i,j)=DBLE(Ainv(i,j))
ENDIF
ENDDO
ENDDO
RETURN
END
***********************************************************************
***********************************************************************
integer function izamax(n,zx,incx)
c
c finds the index of element having max. absolute value.
c jack dongarra, 1/15/85.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex zx(*)
double precision smax
integer i,incx,ix,n
double precision dcabs1
c
izamax = 0
if( n.lt.1 .or. incx.le.0 )return
izamax = 1
if(n.eq.1)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
smax = dcabs1(zx(1))
ix = ix + incx
do 10 i = 2,n
if(dcabs1(zx(ix)).le.smax) go to 5
izamax = i
smax = dcabs1(zx(ix))
5 ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 smax = dcabs1(zx(1))
do 30 i = 2,n
if(dcabs1(zx(i)).le.smax) go to 30
izamax = i
smax = dcabs1(zx(i))
30 continue
return
end
***********************************************************************
***********************************************************************
subroutine zaxpy(n,za,zx,incx,zy,incy)
c
c constant times a vector plus a vector.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex zx(*),zy(*),za
integer i,incx,incy,ix,iy,n
double precision dcabs1
if(n.le.0)return
if (dcabs1(za) .eq. 0.0d0) return
if (incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
zy(iy) = zy(iy) + za*zx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
20 do 30 i = 1,n
zy(i) = zy(i) + za*zx(i)
30 continue
return
end
***********************************************************************
***********************************************************************
subroutine zgefa(a,lda,n,ipvt,info)
integer lda,n,ipvt(1),info
Double Complex a(lda,1)
c
c zgefa factors a Double Complex matrix by gaussian elimination.
c
c zgefa is usually called by zgeco, but it can be called
c directly with a saving in time if rcond is not needed.
c (time for zgeco) = (1 + 9/n)*(time for zgefa) .
c
c on entry
c
c a Double Complex(lda, n)
c the matrix to be factored.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c on return
c
c a an upper triangular matrix and the multipliers
c which were used to obtain it.
c the factorization can be written a = l*u where
c l is a product of permutation and unit lower
c triangular matrices and u is upper triangular.
c
c ipvt integer(n)
c an integer vector of pivot indices.
c
c info integer
c = 0 normal value.
c = k if u(k,k) .eq. 0.0 . this is not an error
c condition for this subroutine, but it does
c indicate that zgesl or zgedi will divide by zero
c if called. use rcond in zgeco for a reliable
c indication of singularity.
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas zaxpy,zscal,izamax
c fortran dabs
c
c internal variables
c
Double Complex t
integer izamax,j,k,kp1,l,nm1
c
Double Complex zdum
double precision cabs1
double precision dreal,dimag
Double Complex zdumr,zdumi
dreal(zdumr) = zdumr
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c
c gaussian elimination with partial pivoting
c
info = 0
nm1 = n - 1
if (nm1 .lt. 1) go to 70
do 60 k = 1, nm1
kp1 = k + 1
c
c find l = pivot index
c
l = izamax(n-k+1,a(k,k),1) + k - 1
ipvt(k) = l
c
c zero pivot implies this column already triangularized
c
if (cabs1(a(l,k)) .eq. 0.0d0) go to 40
c
c interchange if necessary
c
if (l .eq. k) go to 10
t = a(l,k)
a(l,k) = a(k,k)
a(k,k) = t
10 continue
c
c compute multipliers
c
t = -(1.0d0,0.0d0)/a(k,k)
call zscal(n-k,t,a(k+1,k),1)
c
c row elimination with column indexing
c
do 30 j = kp1, n
t = a(l,j)
if (l .eq. k) go to 20
a(l,j) = a(k,j)
a(k,j) = t
20 continue
call zaxpy(n-k,t,a(k+1,k),1,a(k+1,j),1)
30 continue
go to 50
40 continue
info = k
50 continue
60 continue
70 continue
ipvt(n) = n
if (cabs1(a(n,n)) .eq. 0.0d0) info = n
return
end
***********************************************************************
***********************************************************************
subroutine zswap (n,zx,incx,zy,incy)
c
c interchanges two vectors.
c jack dongarra, 3/11/78.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex zx(*),zy(*),ztemp
integer i,incx,incy,ix,iy,n
c
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments not equal
c to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
ztemp = zx(ix)
zx(ix) = zy(iy)
zy(iy) = ztemp
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
20 do 30 i = 1,n
ztemp = zx(i)
zx(i) = zy(i)
zy(i) = ztemp
30 continue
return
end
***********************************************************************
***********************************************************************
double precision function dcabs1(z)
double complex z,zz
double precision t(2)
equivalence (zz,t(1))
zz = z
dcabs1 = dabs(t(1)) + dabs(t(2))
return
end
***********************************************************************
***********************************************************************
subroutine zgedi(a,lda,n,ipvt,det,work,job)
integer lda,n,ipvt(1),job
Double Complex a(lda,1),det(2),work(1)
c
c zgedi computes the determinant and inverse of a matrix
c using the factors computed by zgeco or zgefa.
c
c on entry
c
c a Double Complex(lda, n)
c the output from zgeco or zgefa.
c
c lda integer
c the leading dimension of the array a .
c
c n integer
c the order of the matrix a .
c
c ipvt integer(n)
c the pivot vector from zgeco or zgefa.
c
c work Double Complex(n)
c work vector. contents destroyed.
c
c job integer
c = 11 both determinant and inverse.
c = 01 inverse only.
c = 10 determinant only.
c
c on return
c
c a inverse of original matrix if requested.
c otherwise unchanged.
c
c det Double Complex(2)
c determinant of original matrix if requested.
c otherwise not referenced.
c determinant = det(1) * 10.0**det(2)
c with 1.0 .le. cabs1(det(1)) .lt. 10.0
c or det(1) .eq. 0.0 .
c
c error condition
c
c a division by zero will occur if the input factor contains
c a zero on the diagonal and the inverse is requested.
c it will not occur if the subroutines are called correctly
c and if zgeco has set rcond .gt. 0.0 or zgefa has set
c info .eq. 0 .
c
c linpack. this version dated 08/14/78 .
c cleve moler, university of new mexico, argonne national lab.
c
c subroutines and functions
c
c blas zaxpy,zscal,zswap
c fortran dabs,dcmplx,mod
c
c internal variables
c
Double Complex t
double precision ten
integer i,j,k,kb,kp1,l,nm1
c
Double Complex zdum
double precision cabs1
double precision dreal,dimag
Double Complex zdumr,zdumi
dreal(zdumr) = zdumr
dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
c
c compute determinant
c
if (job/10 .eq. 0) go to 70
det(1) = (1.0d0,0.0d0)
det(2) = (0.0d0,0.0d0)
ten = 10.0d0
do 50 i = 1, n
if (ipvt(i) .ne. i) det(1) = -det(1)
det(1) = a(i,i)*det(1)
c ...exit
if (cabs1(det(1)) .eq. 0.0d0) go to 60
10 if (cabs1(det(1)) .ge. 1.0d0) go to 20
det(1) = dcmplx(ten,0.0d0)*det(1)
det(2) = det(2) - (1.0d0,0.0d0)
go to 10
20 continue
30 if (cabs1(det(1)) .lt. ten) go to 40
det(1) = det(1)/dcmplx(ten,0.0d0)
det(2) = det(2) + (1.0d0,0.0d0)
go to 30
40 continue
50 continue
60 continue
70 continue
c
c compute inverse(u)
c
if (mod(job,10) .eq. 0) go to 150
do 100 k = 1, n
a(k,k) = (1.0d0,0.0d0)/a(k,k)
t = -a(k,k)
call zscal(k-1,t,a(1,k),1)
kp1 = k + 1
if (n .lt. kp1) go to 90
do 80 j = kp1, n
t = a(k,j)
a(k,j) = (0.0d0,0.0d0)
call zaxpy(k,t,a(1,k),1,a(1,j),1)
80 continue
90 continue
100 continue
c
c form inverse(u)*inverse(l)
c
nm1 = n - 1
if (nm1 .lt. 1) go to 140
do 130 kb = 1, nm1
k = n - kb
kp1 = k + 1
do 110 i = kp1, n
work(i) = a(i,k)
a(i,k) = (0.0d0,0.0d0)
110 continue
do 120 j = kp1, n
t = work(j)
call zaxpy(n,t,a(1,j),1,a(1,k),1)
120 continue
l = ipvt(k)
if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1)
130 continue
140 continue
150 continue
return
end
***********************************************************************
***********************************************************************
subroutine zscal(n,za,zx,incx)
c
c scales a vector by a constant.
c jack dongarra, 3/11/78.
c modified 3/93 to return if incx .le. 0.
c modified 12/3/93, array(1) declarations changed to array(*)
c
double complex za,zx(*)
integer i,incx,ix,n
c
if( n.le.0 .or. incx.le.0 )return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
do 10 i = 1,n
zx(ix) = za*zx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
20 do 30 i = 1,n
zx(i) = za*zx(i)
30 continue
return
end
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************
***********************************************************************