Numerical Methods
|
|
Numerical Methods | |
|
|
Fortran libraryFortran has the reputation for delivering the fastest execution speeds and often scientists working on computationally intensive problems use Fortran. Many Fortran libraries containing efficient routines are available. Below is a library of routines that were discussed in this course.
module numeric_lib
! Fortran library containing the routines discussed in the Numerical Methods course at the TU Graz
! Created by Markus Krammer
implicit none
double precision, parameter :: m_pi=4.d0*atan(1.d0) !=3.141592653589793238462643383279502884197169399375105820974944d0
contains
subroutine jacobi(a,n,tmax,eigvec,error)
double precision, intent(inout) :: a(:,:),eigvec(:,:)
integer, intent(in) :: n,tmax
integer, intent(out) :: error
! solves a symmetric real eigenvalue problem a*eigvec=E*eigvec where E are the eigenvalues
! a(:,:): input matrix, eigenvaluses after call in the main diagonal of a(:,:), only upper triangle needed
! n: number of rows / columns of the matrix
! tmax: maximum steps of iteration
! eigvec(:,:): eigenvectors in the columns: eigenvector i = eigvec(:,i)
! error: 0->convergence reached, 1->iteration failed
integer :: i1,i2,i3,t
double precision :: sum_a,limit_a,theta,g,h,tfac,sine,cosine,store_a
error=1
if(n<2) then
eigvec(1,1)=1.0d0
return
end if
eigvec=0.0d0
do i1=1,n
eigvec(i1,i1)=1.0d0
end do
t=0
do while(t.le.tmax)
sum_a=0.0d0
do i1=1,n-1
do i2=i1+1,n
sum_a=sum_a+2.0d0*a(i1,i2)**2
end do
end do
if(sum_a.eq.0.0d0) then
error=0
exit
end if
limit_a=sqrt(sum_a)/n**2
do i1=1,n-1
do i2=i1+1,n
g=1.0d2*abs(a(i1,i2))
if(abs(a(i1,i2)).gt.limit_a) then
h=a(i1,i1)-a(i2,i2)
if((abs(h)+g).eq.abs(h)) then
tfac=a(i1,i2)/h
else
theta=h/(2.0d0*a(i1,i2))
tfac=1.0d0/(abs(theta)+sqrt(1.0d0+theta**2))
if(theta.lt.0.0d0) then
tfac=-tfac
end if
end if
cosine=1.0d0/sqrt(tfac**2+1.0d0)
sine=tfac*cosine
do i3=i1+1,i2-1
store_a=a(i1,i3)
a(i1,i3)=cosine*a(i1,i3)+sine*a(i3,i2)
a(i3,i2)=cosine*a(i3,i2)-sine*store_a
end do
do i3=i2+1,n
store_a=a(i1,i3)
a(i1,i3)=cosine*a(i1,i3)+sine*a(i2,i3)
a(i2,i3)=cosine*a(i2,i3)-sine*store_a
end do
do i3=1,i1-1
store_a=a(i3,i1)
a(i3,i1)=cosine*a(i3,i1)+sine*a(i3,i2)
a(i3,i2)=cosine*a(i3,i2)-sine*store_a
end do
store_a=a(i1,i1)
a(i1,i1)=cosine**2*a(i1,i1)+2.0d0*cosine*sine*a(i1,i2)+sine**2*a(i2,i2)
a(i2,i2)=cosine**2*a(i2,i2)-2.0d0*cosine*sine*a(i1,i2)+sine**2*store_a
a(i1,i2)=0.0d0
do i3=1,n
store_a=eigvec(i3,i2)
eigvec(i3,i2)=cosine*eigvec(i3,i2)-sine*eigvec(i3,i1)
eigvec(i3,i1)=cosine*eigvec(i3,i1)+sine*store_a
end do
end if
end do
end do
t=t+1
end do
if(error.ne.0) then
print '(a,i0,a)','ERROR in jacobi: convergence not reached after ',t,' iterations'
end if
end subroutine jacobi
subroutine tql2(n,d,e,z,error)
integer, intent(in) :: n
double precision, intent(inout) :: d(:),e(:),z(:,:)
integer, intent(out) :: error
! calculate the eigenvalues and eigenvectors of a symmetric tridiagonal matrix by the ql method,
! a full symmetric matrix can also be diagonalized if tred2(...) is used to reduce the problem
! to a tridiagonal problem
! n: number of rows / columns in the matrix
! d(:): conatins the diagonal matrix elements as input and the ascending sorted eigenvalues as output
! e(:): contains the subdiagonal elements, e(1) is arbitraty and e(2:n) holds the (n-1) values,
! e is destroyed during the calculations
! z(:,:): holds the transformation matrix produced by tred2(:) as input and the orthonormal
! eigenvectors as output, sorting is the same as for the eigenvalues,
! if the problem is tridiagonal, z has to be the unity matrix as input
! error: 0->convergence reached, i>0->iteration failed during the calculation of the i-th eigenvalue,
! eigenvalues and -vectors 1 to i-1 are correct but not sorted
integer i1,i2,i3,i4,i5,i41,i42
integer, allocatable :: indx(:)
double precision c1,c2,c3,di41,ei41,f,g,h,p,r,s1,s2,tst1,tst2
integer, parameter :: maxtry=30
logical :: conv
if(n.eq.1) then
error=0
return
end if
e(1:(n-1))=e(2:n)
f=0.0d0
tst1=0.0d0
e(n)=0.0d0
do i4=1,n
error=i4
h=abs(d(i4))+abs(e(i4))
if(tst1.lt.h) tst1=h
do i5=i4,n
tst2=tst1+abs(e(i5))
if(tst2.eq.tst1) exit
end do
if(i5.ne.i4) then
conv=.false.
do i2=1,maxtry
i41=i4+1
i42=i41+1
g=d(i4)
p=(d(i41)-g)/(2.0d0*e(i4))
r=sign(sqrt(p**2+1.0d0),p)+p
d(i4)=e(i4)/r
d(i41)=e(i4)*r
di41=d(i41)
h=g-d(i4)
if(i42.le.n) then
do i1=i42,n
d(i1)=d(i1)-h
end do
end if
f=f+h
p=d(i5)
c1=1.0d0
c2=c1
ei41=e(i41)
s1=0.0d0
do i1=i5-1,i4,-1
c3=c2
c2=c1
s2=s1
g=c1*e(i1)
h=c1*p
r=sqrt(p**2+e(i1)**2)
e(i1+1)=s1*r
s1=e(i1)/r
c1=p/r
p=c1*d(i1)-s1*g
d(i1+1)=h+s1*(c1*g+s1*d(i1))
do i3=1,n
h=z(i3,i1+1)
z(i3,i1+1)=s1*z(i3,i1)+c1*h
z(i3,i1)=c1*z(i3,i1)-s1*h
end do
end do
p=-s1*s2*c3*ei41*e(i4)/di41
e(i4)=s1*p
d(i4)=c1*p
tst2=tst1+abs(e(i4))
if(tst2.le.tst1) then
conv=.true.
exit
end if
end do
if(.not.conv) then
print '(a,i0,a)','ERROR in tql2: only ',error-1,' eigenvalues were determined'
return
end if
end if
d(i4)=d(i4)+f
end do
allocate(indx(n))
call hpsort_indx(n,d,indx)
z=z(:,indx)
deallocate(indx)
end subroutine tql2
subroutine tred2(n,a,d,e)
integer, intent(in) :: n
double precision, intent(inout) :: a(:,:)
double precision, intent(out) :: d(:),e(:)
! reduce a symmetric real eigenvalue problem a*eigvec=E*eigvec (where E are the eigenvalues)
! to a symmetric tridiagonal eigenvalue problem which can be solved by tql2(...) using
! orthogonal similarity transformation
! n: number of rows / columns of the matrix a
! a(:,:): symmetric matrix of the eigenvalueproblem as input and transformation matrix z(:,:)
! for tql2(...) as output, only the lower triangle has to be provided
! d(:): diagonal elements of the tridiagonal matrix
! e(:): subdiagonal elements of the tridiagonal matrix on sites 2 to n
integer :: i1,i2,i3,i4,i2p1
double precision :: f,g,h,hh,scale
d=a(n,:)
if(n.eq.1) then
return
end if
do i1=n,2,-1
i4=i1-1
h=0.0d0
scale=0.0d0
if(i4.ge.2) then
do i3=1,i4
scale=scale+abs(d(i3))
end do
end if
if(scale.eq.0.0d0) then
e(i1)=d(i4)
do i2=1,i4
d(i2)=a(i4,i2)
a(i1,i2)=0.0d0
a(i2,i1)=0.0d0
end do
else
do i3=1,i4
d(i3)=d(i3)/scale
h=h+d(i3)**2
end do
f=d(i4)
g=-sign(sqrt(h),f)
e(i1)=scale*g
h=h-f*g
d(i4)=f-g
do i2=1,i4
e(i2)=0.0d0
end do
do i2=1,i4
f=d(i2)
a(i2,i1)=f
g=e(i2)+a(i2,i2)*f
i2p1=i2+1
if(i4.ge.i2p1) then
do i3=i2p1,i4
g=g+a(i3,i2)*d(i3)
e(i3)=e(i3)+a(i3,i2)*f
end do
end if
e(i2)=g
end do
f=0.0d0
do i2=1,i4
e(i2)=e(i2)/h
f=f+e(i2)*d(i2)
end do
hh=f/(2.0d0*h)
do i2=1,i4
e(i2)=e(i2)-hh*d(i2)
end do
do i2=1,i4
f=d(i2)
g=e(i2)
do i3=i2,i4
a(i3,i2)=a(i3,i2)-f*e(i3)-g*d(i3)
end do
d(i2)=a(i4,i2)
a(i1,i2)=0.0d0
end do
end if
d(i1)=h
end do
do i1=2,n
i4=i1-1
a(n,i4)=a(i4,i4)
a(i4,i4)=1.0d0
h=d(i1)
if(h.ne.0.0d0) then
do i3=1,i4
d(i3)=a(i3,i1)/h
end do
do i2=1,i4
g=0.0d0
do i3=1,i4
g=g+a(i3,i1)*a(i3,i2)
end do
do i3=1,i4
a(i3,i2)=a(i3,i2)-g*d(i3)
end do
end do
end if
do i3=1,i4
a(i3,i1)=0.0d0
end do
end do
do i1=1,n
d(i1)=a(n,i1)
a(n,i1)=0.0d0
end do
a(n,n)=1.0d0
e(1)=0.0d0
end subroutine tred2
subroutine tdiag(n,a,eigval)
integer, intent(in) :: n
double precision, intent(inout) :: a(:,:)
double precision, intent(out) :: eigval(:)
! solve real symmetric eigenvalue problem a*x=eigval*x with tred2(...) and tql2(...)
! n: number of rows / columns of the matrix a
! a(:,:): real symmetric matrix of the eigenvalue problem as input and eigenvectors as output
! eigenvectors are found in the columns: eigenvector i = a(:,i)
! eigval(:): eigenvalues of the initial matrix a
double precision, allocatable :: e(:)
integer :: error
allocate(e(n))
call tred2(n,a,eigval,e)
call tql2(n,eigval,e,a,error)
deallocate(e)
end subroutine tdiag
subroutine reduc(n,a,b,error)
integer, intent(in) :: n
double precision, intent(inout) :: a(:,:),b(:,:)
integer, intent(out) :: error
! reduce the generalized symmetric eigenvalue problem a*x=E*b*x with eigenvectors x and eigenvalues E
! to a standard eigenvalue problem using cholesky factorization of matrix b
! n: number of rows / columns of matrix a and b
! a(:,:): input: real symmetric matrix a(:,:) (only full upper triangle needs to be provided)
! output: symmetric matrix of the standard eigenvalue problem in the full lower triangle, the
! strict upper triangle is not changed
! b(:,:): input: real symmetric and positiv definite matrix b(:,:) (only full upper triangle needs to
! be provided)
! output: cholesky factor l (symmetric matrix) in the full lower diagonal, the strict upper
! triangle is not changed
! error: 0->reduction was successful, 1->b(:,:) is not positive definite
integer :: i1,i2,i3,i1m1,i2m1
double precision :: x,y
double precision, allocatable :: dl(:)
allocate(dl(n))
error=0
do i1=1,n
i1m1=i1-1
do i2=i1,n
x=b(i1,i2)
if(i1.ne.1) then
do i3=1,i1m1
x=x-b(i1,i3)*b(i2,i3)
end do
end if
if(i2.eq.i1) then
if(x.le.0.0d0) then
print '(a)','ERROR in reduc: b is not positive definite'
error=1
return
end if
y=sqrt(x)
dl(i1)=y
else
b(i2,i1)=x/y
end if
end do
end do
do i1=1,n
i1m1=i1-1
y=dl(i1)
do i2=i1,n
x=a(i1,i2)
if(i1.gt.1) then
do i3=1,i1m1
x=x-b(i1,i3)*a(i2,i3)
end do
end if
a(i2,i1)=x/y
end do
end do
do i2=1,n
i2m1=i2-1
do i1=i2,n
x=a(i1,i2)
if(i1.ne.i2) then
do i3=i2,i1-1
x=x-a(i3,i2)*b(i1,i3)
end do
end if
if(i2.gt.1) then
do i3=1,i2m1
x=x-a(i2,i3)*b(i1,i3)
end do
end if
a(i1,i2)=x/dl(i1)
end do
end do
do i1=1,n
b(i1,i1)=dl(i1)
end do
deallocate(dl)
end subroutine reduc
subroutine eigsys(n,a,b,eigval)
integer, intent(in) :: n
double precision, intent(inout) :: a(:,:),b(:,:)
double precision, intent(out) :: eigval(:)
! solves a generalized eigenvalue problem a*x=eigval*b*x with eigenvectors x
! n: number of rows / columns of the matricies a(:,:) and b(:,:)
! a(:,:): input: real symmetric matrix a(:,:)
! output: eigenvectors in the columns: eigenvector i = a(:,i)
! b(:,:): real symmetric positiv defined matrix b(:,:) as input
! values modified during diagonalization
! eigval(:): eigenvalues sorted corresponding to the eigenvectors
integer :: i1,i2,i3,error
double precision :: sum
call reduc(n,a,b,error)
if(error.ne.0) return
call tdiag(n,a,eigval)
do i3=1,n
a(n,i3)=a(n,i3)/b(n,n)
do i1=n-1,1,-1
sum=a(i1,i3)
do i2=i1+1,n
sum=sum-b(i2,i1)*a(i2,i3)
end do
a(i1,i3)=sum/b(i1,i1)
end do
end do
end subroutine eigsys
subroutine gausei(n,ndiag,s,diag,f,tmax,irel,acc,sol,error,w,t0,t_out)
integer, intent(in) :: n,ndiag,s(:),tmax,irel
double precision, intent(in) :: diag(:,:),f(:),acc
double precision, intent(out) :: sol(:)
integer, intent(out) :: error
double precision, optional, intent(inout) :: w
integer, optional, intent(in) :: t0
integer, optional, intent(out) :: t_out
! solves the equation system diag*sol=f for sparse matricies
! n: number of rows / columns of the matrix
! ndiag: number of filled diagonals
! s(:): position of diagonals, 0 is the main diagonal, +/-i is the i-th upper/lower secondary diagonal
! diag(:,:): matrix stored in sparse form
! f(:): inhomogeniety vector
! tmax: maximum number of iterations
! irel: 1->relative accuracy, else->absolute accuracy
! acc: desired accuracy
! sol(:): solution vector
! error: 0->convergence reached, 1->iteration failed
! w: optimizationparameter omega
! t0: step where an omega optimization should take place
! t_out: iterations needed to finish
integer :: main_diag,t,i1,i2,i3,t0_intern
double precision :: s1,dx,sum_dx,old_sum_dx,err_dx,w_intern
if(present(w)) then
w_intern=w
else
w_intern=1.0d0
end if
if(present(t0)) then
t0_intern=t0
w_intern=1.0d0
else
t0_intern=tmax+2
end if
error=1
main_diag=0
do i1=1,ndiag
if(s(i1).eq.0) then
main_diag=i1
end if
end do
if(main_diag.eq.0) then
print '(a)','FATAL ERROR in gausei: no main diagonal found'
return
end if
do i1=1,n
if(diag(i1,main_diag).eq.0.0d0) then
print '(a)','FATAL ERROR in gausei: main diagonal contains zeros'
return
end if
end do
sol=0.0d0
old_sum_dx=0.0d0
sum_dx=0.0d0
t=0
do while((t.lt.tmax).and.(error.eq.1))
error=0
do i1=1,n
s1=0.0d0
do i2=1,ndiag
i3=s(i2)+i1
if((i2.ne.main_diag).and.(i3.le.n)) then
s1=s1+diag(i1,i2)*sol(i3)
end if
end do
dx=w_intern*((s1-f(i1))/diag(i1,main_diag)+sol(i1))
if(t.eq.t0_intern-1) then
old_sum_dx=old_sum_dx+dx**2
else if(t.eq.t0_intern) then
sum_dx=sum_dx+dx**2
end if
sol(i1)=sol(i1)-dx
if(irel.eq.1) then
err_dx=abs(dx/sol(i1))
else
err_dx=abs(dx)
end if
if(err_dx.gt.acc) then
error=1
end if
end do
if(t.eq.t0_intern) then
w_intern=2.0d0/(1.0d0+sqrt(1.0d0-sqrt(sum_dx/old_sum_dx)))
end if
t=t+1
end do
if(error.eq.1) then
print '(a,i0,a)','ERROR in gausei: convergence not reached after ',tmax,' iterations'
end if
if(present(t_out)) then
t_out=t
end if
if(present(w)) then
w=w_intern
end if
end subroutine gausei
subroutine ludcmp(a,n,indx,d,khad)
double precision, intent(inout) :: a(:,:)
integer, intent(in) :: n
integer, intent(out) :: indx(:)
double precision, intent(out) :: d,khad
! perform a LU decomposition of a real quadratic matrix
! a(:,:): input matrix and LU decomposited output matrix
! n: number of rows / columns of the matrix
! indx(:): index of the swaps lines
! d: calculated determinant
! khad: calculated Hadramard condition number
integer :: i1,i2,i3,imax
double precision :: sum,aamax,dum
double precision, parameter :: ctiny=1.0d-20
khad=1.d0
do i1=1,n
sum=0.0d0
do i2=1,n
sum=sum+a(i1,i2)**2
end do
khad=khad*sqrt(sum)
end do
d=1.d0
do i2=1,n
do i1=1,i2-1
sum=a(i1,i2)
do i3=1,i1-1
sum=sum-a(i1,i3)*a(i3,i2)
end do
a(i1,i2)=sum
end do
aamax=0.0d0
do i1=i2,n
sum=a(i1,i2)
do i3=1,i2-1
sum=sum-a(i1,i3)*a(i3,i2)
end do
a(i1,i2)=sum
dum=abs(sum)
if(dum.ge.aamax) then
imax=i1
aamax=dum
end if
end do
if(i2.ne.imax) then
do i3=1,n
dum=a(imax,i3)
a(imax,i3)=a(i2,i3)
a(i2,i3)=dum
end do
d=-d
end if
indx(i2)=imax
if(a(i2,i2).eq.0.0d0) then
a(i2,i2)=ctiny
end if
dum=1.0d0/a(i2,i2)
do i1=i2+1,n
a(i1,i2)=a(i1,i2)*dum
end do
d=d*a(i2,i2)
end do
khad=abs(d)/khad
end subroutine ludcmp
subroutine lubksb(a,n,indx,b)
double precision, intent(in) :: a(:,:)
integer, intent(in) :: n,indx(:)
double precision, intent(inout) :: b(:)
! calculate the solution vector a*x=b for a LU decomposed matrix a
! a(:,:): LU decomposited matrix a
! n: number of rows / columns of the matrix
! indx(:): index of the swaped lines
! b(:): inhomogeniety vector b as input and solutionvector x as output
integer :: i1,i2,ii,ll
double precision :: sum
ii=0
do i1=1,n
ll=indx(i1)
sum=b(ll)
b(ll)=b(i1)
if(ii.ne.0) then
do i2=ii,i1-1
sum=sum-a(i1,i2)*b(i2)
end do
else if(sum.ne.0.0d0) then
ii=i1
end if
b(i1)=sum
end do
do i1=n,1,-1
sum=b(i1)
do i2=i1+1,n
sum=sum-a(i1,i2)*b(i2)
end do
b(i1)=sum/a(i1,i1)
end do
end subroutine lubksb
subroutine luinv(a,n,det,khad,woff)
double precision, intent(inout) :: a(:,:)
integer, intent(in) :: n
double precision, intent(out) :: det,khad
logical, optional, intent(in) :: woff
! invert a real quadratic matrix a with LU decomposition
! a(:,:): input matrix and inverted output matrix
! n: number of rows / columns of the matrix
! det: derterminant of the matris
! khad: Hadramard condition number of the matrix
! woff: .true.->warnings are turned off, .false. or not present->warnings are displayed
integer :: i1
integer, allocatable :: indx(:)
double precision, allocatable :: mat_LU(:,:),b(:)
allocate(indx(n))
allocate(mat_LU(n,n))
allocate(b(n))
mat_LU=a
call ludcmp(mat_LU,n,indx,det,khad)
if((.not.present(woff)).or.(.not.woff)) then
if(khad.lt.1.0d-2) then
print '(a,f12.8)','WARNING in luinv: Hadramard condition number is too low: ',khad
else if(khad.lt.1.0d-1) then
print '(a,f12.8)','WARNING in luinv: Hadramard condition number is low: ',khad
end if
end if
do i1=1,n
b=0.0d0
b(i1)=1.0d0
call lubksb(mat_LU,n,indx,b)
a(:,i1)=b
end do
deallocate(indx)
deallocate(mat_LU)
deallocate(b)
end subroutine luinv
subroutine lusol(a,n,b,x,woff)
double precision, intent(in) :: a(:,:),b(:)
integer, intent(in) :: n
double precision, intent(out) :: x(:)
logical, optional, intent(in) :: woff
! solve the equation system a*x=b of a real quadratic matrix a with LU decomposition
! a(:,:): input matrix
! n: number of rows / columns of the matrix
! b(:): inhomogeniety vector
! x(:): solution vector
! woff: .true.->warnings are turned off, .false. or not present->warnings are displayed
integer, allocatable :: indx(:)
double precision :: det,khad
double precision, allocatable :: mat_LU(:,:)
allocate(indx(n))
allocate(mat_LU(n,n))
mat_LU=a
call ludcmp(mat_LU,n,indx,det,khad)
if((.not.present(woff)).or.(.not.woff)) then
if(khad.lt.1.0d-2) then
print '(a,f12.8)','WARNING in lusol: Hadramard condition number is too low: ',khad
else if(khad.lt.1.0d-1) then
print '(a,f12.8)','WARNING in lusol: Hadramard condition number is low: ',khad
end if
end if
x=b
call lubksb(mat_LU,n,indx,x)
deallocate(indx)
deallocate(mat_LU)
end subroutine lusol
subroutine trid(b,a,c,r,n)
double precision, intent(inout) :: b(:),r(:)
double precision, intent(in) :: a(:),c(:)
integer, intent(in) :: n
! solving an equation system of a tridiagonal matrix
! b(:): main diaglnal of the tridaigonal matrix, changed during process
! a(:): first lower minor diagonal stored in a(2:n)
! c(:): first upper minor diagonal stored in c(1:(n-1))
! r(:): inhomogeniety vector as input and solution vector as output
! n: number of rows / columns of the matrix
double precision :: m
integer :: i1
if(b(1).eq.0.0d0) then
print '(a)','ERROR in trid: matrix is singular'
return
end if
if(n.lt.2) then
r(1)=r(1)/b(1)
return
end if
do i1=2,n
m=a(i1)/b(i1-1)
b(i1)=b(i1)-m*c(i1-1)
if(b(i1).eq.0.0d0) then
print '(a)','ERROR in trid: matrix is singular or close to singularity'
return
end if
r(i1)=r(i1)-m*r(i1-1)
end do
r(n)=r(n)/b(n)
do i1=(n-1),1,-1
r(i1)=(r(i1)-c(i1)*r(i1+1))/b(i1)
end do
end subroutine trid
subroutine hpsort(n,ra)
integer, intent(in) :: n
double precision, intent(inout) :: ra(:)
! sort data in vector ra with heapsort algorithm and return the sorted data in vector ra
! n: number of elements in vector ra
! ra(:): vector to be sorted
integer :: i1,i2,imax,isort
double precision :: rra
if(n.lt.2) then
return
end if
isort=n/2+1
imax=n
do
if(isort.gt.1) then
isort=isort-1
rra=ra(isort)
else
rra=ra(imax)
ra(imax)=ra(1)
imax=imax-1
if(imax.eq.1) then
ra(1)=rra
exit
end if
end if
i1=isort
i2=isort*2
do while(i2.le.imax)
if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
i2=i2+1
end if
if(rra.lt.ra(i2)) then
ra(i1)=ra(i2)
i1=i2
i2=i2*2
else
i2=imax+1
end if
end do
ra(i1)=rra
end do
end subroutine hpsort
subroutine hpsort_indx(n,ra,indx)
integer, intent(in) :: n
double precision, intent(inout) :: ra(:)
integer, intent(out) :: indx(:)
! sort data in vector ra with heapsort algorithm and return the sorted data in vector ra
! n: number of elements in vector ra
! ra(:): vector to be sorted
! indx(:): index where the sorted elements were found in the old array
integer :: i1,i2,imax,isort,iindx
double precision :: rra
do i1=1,n
indx(i1)=i1
end do
if(n.lt.2) then
return
end if
isort=n/2+1
imax=n
do
if(isort.gt.1) then
isort=isort-1
rra=ra(isort)
iindx=indx(isort)
else
rra=ra(imax)
ra(imax)=ra(1)
iindx=indx(imax)
indx(imax)=indx(1)
imax=imax-1
if(imax.eq.1) then
ra(1)=rra
indx(1)=iindx
exit
end if
end if
i1=isort
i2=isort*2
do while(i2.le.imax)
if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
i2=i2+1
end if
if(rra.lt.ra(i2)) then
ra(i1)=ra(i2)
indx(i1)=indx(i2)
i1=i2
i2=i2*2
else
i2=imax+1
end if
end do
ra(i1)=rra
indx(i1)=iindx
end do
end subroutine hpsort_indx
subroutine hpsort_indx_only(n,ra_in,indx)
integer, intent(in) :: n
double precision, intent(in) :: ra_in(:)
integer, intent(out) :: indx(:)
! sort data in vector ra with heapsort algorithm and return the indices of sorting in indx
! n: number of elements in vector ra
! ra_in(:): vector to be sorted
! indx(:): index where the sorted elements were found in the old array
integer :: i1,i2,imax,isort,iindx
double precision :: rra,ra(n)
do i1=1,n
ra(i1)=ra_in(i1)
indx(i1)=i1
end do
if(n.lt.2) then
return
end if
isort=n/2+1
imax=n
do
if(isort.gt.1) then
isort=isort-1
rra=ra(isort)
iindx=indx(isort)
else
rra=ra(imax)
ra(imax)=ra(1)
iindx=indx(imax)
indx(imax)=indx(1)
imax=imax-1
if(imax.eq.1) then
ra(1)=rra
indx(1)=iindx
exit
end if
end if
i1=isort
i2=isort*2
do while(i2.le.imax)
if((i2.lt.imax).and.(ra(i2).lt.ra(i2+1))) then
i2=i2+1
end if
if(rra.lt.ra(i2)) then
ra(i1)=ra(i2)
indx(i1)=indx(i2)
i1=i2
i2=i2*2
else
i2=imax+1
end if
end do
ra(i1)=rra
indx(i1)=iindx
end do
end subroutine hpsort_indx_only
subroutine nestint_sing(fct,x0,x1,acc_x,acc_fct,xzero)
interface
double precision function fct(x)
double precision, intent(in) :: x
end function fct
end interface
double precision, intent(in) :: x0,x1
double precision, intent(in) :: acc_x,acc_fct
double precision, intent(out) :: xzero
! find one root of a function fct in the interval [x0,x1] with nested intervals
! fct: pure function declared by the input
! x0: lower boundary of the search interval
! x1: higher boundary of the search interval
! acc_x: desired absolute accuracy for xzero
! acc_fct: required minimum absolute value of the function
! xzero: found root
double precision :: fct0,fct1,x00,x11
if(x0.ge.x1) then
print '(a)','ERROR in nestint_sing: x0 is not less than x1'
xzero=x0
return
end if
fct0=fct(x0)
fct1=fct(x1)
if((fct0*fct1).gt.0.0d0) then
print '(a)','ERROR in nestint_sing: no root in the given interval'
xzero=x0
return
end if
if((fct0.eq.0.0d0).and.(fct1.eq.0.0d0)) then
print '(a)','WARNING in nestint_sing: the function value on both sides of the search interval are zero'
xzero=(x0+x1)/2.0d0
return
end if
x00=x0
x11=x1
do
xzero=(x00+x11)/2.0d0
fct1=fct(xzero)
if((abs(fct1).lt.acc_fct).and.((x11-x00).lt.acc_x)) then
xzero=(x00+x11)/2.0d0
exit
end if
if((fct0*fct1).le.0.0d0) then
x11=xzero
else
fct0=fct1
x00=xzero
end if
end do
end subroutine nestint_sing
subroutine nestint_mult(fct,xstart,xend,step,acc_x,acc_fct,nzerosmax,xzeros,nzeros)
interface
double precision function fct(x)
double precision, intent(in) :: x
end function fct
end interface
double precision, intent(in) :: xstart,xend,step,acc_x,acc_fct
integer, intent(in) :: nzerosmax
double precision, intent(out) :: xzeros(:)
integer, intent(out) :: nzeros
! find roots of a function fct in the interval [x0,x1] with nested intervals
! fct: pure function declared by the input
! xstart: lower boundary of the search interval
! xend: higher boundary of the search interval
! step: stepsize for scanning the search interval
! acc_x: desired absolute accuracy for xzero
! acc_fct: required minimum absolute value of the function
! nzerosmax: maximum number of roots that can be found
! xzeros(:): found roots
! nzeros: found number of roots
double precision :: x0,x1,fct0,fct1
x0=xstart
fct0=fct(x0)
nzeros=0
do while(x0.lt.xend)
x1=x0+step
fct1=fct(x1)
if((fct0*fct1).le.0.0d0) then
nzeros=nzeros+1
call nestint_sing(fct,x0,x1,acc_x,acc_fct,xzeros(nzeros))
if(nzeros.ge.nzerosmax) then
print '(a)','WARNING in nestint_mult: maximum number of roots found'
exit
end if
end if
x0=x1
fct0=fct1
end do
xzeros(nzeros+1:nzerosmax)=0.0d0
end subroutine nestint_mult
subroutine gauleg(x1,x2,x,w,n)
double precision, intent(in) :: x1,x2
double precision, intent(out) :: x(:),w(:)
integer, intent(in) :: n
! calculating the roots of the gauss legendre polynome and the weighting factors for a gauss quadrature
! x1: start of the interval
! x2: end of the interval
! x(:): gauss legendre points
! w(:): weighting factors
! n: number of points to calculate
integer :: i1,i2,m
double precision :: z1,z,xm,xl,pp,p1,p2,p3
double precision, parameter :: m_eps=3.0d-14
m=(n+1)/2
xm=0.5d0*(x2+x1)
xl=0.5d0*(x2-x1)
do i1=1,m
z=cos(m_pi*(i1-0.25d0)/(n+0.5d0))
do
p1=1.0d0
p2=0.0d0
do i2=1,n
p3=p2
p2=p1
p1=((2.0d0*i2-1.0d0)*z*p2-(i2-1.0d0)*p3)/i2
end do
pp=n*(z*p1-p2)/(z**2-1.0d0)
z1=z
z=z1-p1/pp
if(abs(z-z1).lt.m_eps) then
exit
end if
end do
x(i1)=xm-xl*z
x(n+1-i1)=xm+xl*z
w(i1)=2.0d0*xl/((1.0d0-z**2)*pp**2)
w(n+1-i1)=w(i1)
end do
end subroutine gauleg
subroutine gauint(fct,x0,x1,irel,acc,ostart,oinc,omax,intfct,error,woff)
interface
double precision function fct(x)
double precision, intent(in) :: x
end function fct
end interface
double precision, intent(in) :: x0,x1,acc
integer, intent(in) :: irel,ostart,oinc,omax
double precision, intent(out) :: intfct
integer, intent(out) :: error
logical, optional, intent(in) :: woff
! integrate the function fct from x0 to x1 with gauss legendre quadrature
! fct: function to integrate
! x0: start of integration interval
! x1: end of integration interval
! irel: 1->relative accuracy, else->absolute accuracy
! acc: desired accuracy
! ostart: starting order of the quadrature
! oinc: increment of the order of the quadrature
! omax: maximum order of the quadrature
! intfct: value of the integral
! error: 0->convergence reached, 1->quadrature failed
integer :: order,i1
double precision :: oldintfct,err_acc
double precision, allocatable :: w(:),z(:)
error=1
if(ostart.lt.1) then
print '(a)','ERROR in gauint: starting order < 1'
return
end if
if(oinc.lt.1) then
print '(a)','ERROR in gauint: order increment < 1'
return
end if
if(omax.lt.ostart) then
print '(a)','ERROR in gauint: starting order < maximum order'
return
end if
order=ostart
do while(order.le.omax)
allocate(z(order))
allocate(w(order))
call gauleg(x0,x1,z,w,order)
intfct=0.0d0
do i1=1,order
intfct=intfct+w(i1)*fct(z(i1))
end do
deallocate(z)
deallocate(w)
if(order.ne.ostart) then
if((oldintfct.eq.0.0d0).and.(intfct.eq.0.0d0)) then
if((.not.present(woff)).or.(.not.woff)) print '(a)','WARNING in gauint: integral is zero for every evaluation'
error=0
exit
end if
err_acc=abs(oldintfct-intfct)
if(irel.eq.1) then
err_acc=err_acc/abs(oldintfct)
end if
if(err_acc.lt.acc) then
error=0
exit
end if
end if
oldintfct=intfct
order=order+oinc
end do
if(error.ne.0) then
print '(a)','ERROR in gauint: convergence not reached'
end if
end subroutine gauint
subroutine trapzd(fct,x0,x1,intfct,j)
interface
double precision function fct(x)
double precision, intent(in) :: x
end function fct
end interface
double precision, intent(in) :: x0,x1
double precision, intent(inout) :: intfct
integer, intent(in) :: j
! integrate the function fct from x0 to x1 with trapez methode, j has to start from 1
! fct: function to integrate
! x0: start of integration interval
! x1: end of integration interval
! intfct: value of the integral
! j: 2**(j-2) new values are added to the integral, j=1 is the initialization for the margins
! of the integration interval, j has to increase by 1 every time it is called
integer :: i1,nint
double precision :: del,x,sum
if(j.le.1) then
intfct=0.5d0*(x1-x0)*(fct(x0)+fct(x1))
else
nint=2**(j-2)
del=(x1-x0)/dble(nint)
x=x0+0.5d0*del
sum=0.0d0
do i1=1,nint
sum=sum+fct(x)
x=x+del
end do
intfct=0.5d0*(intfct+sum*(x1-x0)/dble(nint))
end if
end subroutine trapzd
subroutine romb(fct,x0,x1,irel,acc,omin,omax,intfct,error)
interface
double precision function fct(x)
double precision, intent(in) :: x
end function fct
end interface
double precision, intent(in) :: x0,x1,acc
integer, intent(in) :: irel,omin,omax
double precision, intent(out) :: intfct
integer, intent(out) :: error
! integrate the function fct from x0 to x1 with gauss legendre quadrature
! fct: function to integrate
! x0: start of integration interval
! x1: end of integration interval
! irel: 1->relative accuracy, else->absolute accuracy
! acc: desired accuracy
! omin: order where the error calculation starts
! omax: maximum order of the quadrature
! intfct: value of the integral
! error: 0->convergence reached, 1->quadrature failed
integer :: i1,ocurr,omax_intern
double precision :: fac,err_acc
double precision, allocatable :: quad(:),quadold(:)
if(omax.lt.3) then
omax_intern=3
else
omax_intern=omax
end if
allocate(quad(2:omax_intern))
allocate(quadold(2:omax_intern))
error=1
call trapzd(fct,x0,x1,quad(2),1)
do ocurr=3,omax_intern
quadold(2:(ocurr-1))=quad(2:(ocurr-1))
call trapzd(fct,x0,x1,quad(2),ocurr-1)
fac=1.0d0
do i1=2,ocurr-1
fac=fac*ocurr
quad(i1+1)=(dble(fac)*quad(i1)-quadold(i1))/(dble(fac-1))
end do
if(ocurr.ge.omin) then
err_acc=abs(quadold(ocurr-1)-quad(ocurr))
if(irel.eq.1) err_acc=err_acc/abs(quadold(ocurr-1))
if(err_acc.lt.acc) then
intfct=quad(ocurr)
error=0
exit
end if
end if
end do
if(error.ne.0) then
print '(a)','ERROR in romb: convergence not reached'
end if
deallocate(quad)
deallocate(quadold)
end subroutine romb
subroutine mrqcof(fct,x,y,sig,ndata,na,a,alpha,beta,chisq)
interface
subroutine fct(x,a,y,dyda)
double precision, intent(in) :: x,a(:)
double precision, intent(out) :: y,dyda(:)
end subroutine fct
end interface
double precision, intent(in) :: x(:),y(:),sig(:),a(:)
integer, intent(in) :: ndata,na
double precision, intent(out) :: alpha(:,:),beta(:),chisq
! calculate the Marquardt coefficients
! fct: model function provided by the user
! x: x value where the function should be evaluated
! a(:): model parameters
! y: function value
! dyda(:): derivatives of the function with respect to the model parameters
! x(:): measured x data
! y(:): measured y data
! sig(:): uncertainty of the measured y data
! ndata: number of data points
! na: number of fit parameters
! a(:): guessed fitparameters as input and improved fit parameters as output
! alpha(:,:): matrix of the linearized equation system
! beta(:): inhomogeniety vector of the linearized equation system
! chisq: calculated chi-squared
integer :: i1,i2,i3
double precision :: ymod,g,dy
double precision, allocatable :: dyda(:)
allocate(dyda(na))
alpha=0.0d0
beta=0.0d0
chisq=0.0d0
do i1=1,ndata
call fct(x(i1),a,ymod,dyda)
g=1.0d0/(sig(i1)**2)
dy=y(i1)-ymod
do i2=1,na
do i3=1,i2
alpha(i2,i3)=alpha(i2,i3)+g*dyda(i2)*dyda(i3)
end do
beta(i2)=beta(i2)+g*dy*dyda(i2)
end do
chisq=chisq+g*dy**2
end do
do i1=2,na
do i2=1,i1-1
alpha(i2,i1)=alpha(i1,i2)
end do
end do
deallocate(dyda)
end subroutine mrqcof
subroutine mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
interface
subroutine fct(x,a,y,dyda)
double precision, intent(in) :: x,a(:)
double precision, intent(out) :: y,dyda(:)
end subroutine fct
end interface
double precision, intent(in) :: x(:),y(:),sig(:)
integer, intent(in) :: ndata,na
double precision, intent(inout) :: a(:),alpha(:,:),beta(:),ochisq,alambda
double precision, intent(out) :: chisq,covar(:,:)
logical :: succmar
! performs one step of a chi-squared optimization of a non linear fit function
! with a Gauss-Newton-Marquardt algorithm - the function fct must provide the
! functionvalues y and the derivatives with respect to the fit parameters
! fct: model function provided by the user
! x: x value where the function should be evaluated
! a(:): model parameters
! y: function value
! dyda(:): derivatives of the function with respect to the model parameters
! x(:): measured x data
! y(:): measured y data
! sig(:): uncertainty of the measured y data
! ndata: number of data points
! na: number of fit parameters
! a(:): guessed fitparameters as input and improved fit parameters as output
! alpha(:,:): internal matrix, the matrix from the last iteration step is needed
! covar(:,:): covariance matrix calculated if alambda is exactly 0.0d0
! beta(:): internal vector, the vector from the last iteration step is needed
! chisq: current chi-squared
! ochisq: chi-squared from the last iteration step
! alambda: Marquardt parameter lambda
! alambda < 0.0d0: the preparatory step is done
! alambda > 0.0d0: normal gnm iteration is performed
! alambda = 0.0d0: the covariance matrix is calculated
! succmar: true: iteration was succesful - false: iteration failed
integer :: i1
double precision :: det,khad
double precision, allocatable :: atry_da(:),vector(:),normal(:,:)
character(len=50) :: formatstring
if(alambda.lt.0.0d0) then
alambda=0.3d-3
call mrqcof(fct,x,y,sig,ndata,na,a,alpha,beta,chisq)
ochisq=chisq
print '(a)','iteration chi**2 parameters a'
write(formatstring,'(a,i0,a)') '(i9,e13.4,',na,'e13.4)'
print formatstring,0,chisq,a
end if
if(alambda==0.0d0) then
covar=alpha
call luinv(covar,na,det,khad,woff=.true.)
else
allocate(normal(na,na))
normal=alpha
do i1=1,na
normal(i1,i1)=alpha(i1,i1)*(1.0d0+alambda)
end do
allocate(atry_da(na))
call lusol(normal,na,beta,atry_da,woff=.true.)
atry_da=a+atry_da
allocate(vector(na))
call mrqcof(fct,x,y,sig,ndata,na,atry_da,normal,vector,chisq)
if(chisq.le.ochisq) then
alambda=alambda*0.2d0
ochisq=chisq
succmar=.true.
alpha=normal
beta=vector
a=atry_da
else
alambda=alambda*5.0d0
succmar=.false.
end if
deallocate(atry_da)
deallocate(vector)
deallocate(normal)
end if
end subroutine mrqmin
subroutine mrqiter(fct,x,y,sig,acc,tmax,ndata,na,a,da,covar,modvar,modvar_confid,error)
interface
subroutine fct(x,a,y,dyda)
double precision, intent(in) :: x,a(:)
double precision, intent(out) :: y,dyda(:)
end subroutine fct
end interface
double precision, intent(in) :: x(:),y(:),sig(:),acc
integer, intent(in) :: ndata,na,tmax
double precision, intent(inout) :: a(:)
double precision, intent(out) :: da(:),covar(:,:),modvar,modvar_confid(2)
integer, intent(out) :: error
! fit a function fct with model parameters a to data (x,y+/-sig)
! fct: model function provided by the user
! x: x value where the function should be evaluated
! a(:): model parameters
! y: function value
! dyda(:): derivatives of the function with respect to the model parameters
! x(:): x data
! y(:): y data
! sig(:): error of the y data
! acc: desired accuracy of the model parameters
! tmax: maximum number of iterations
! ndata: number of data points
! na: number of model parameters
! a(:): model parameters
! da(:): error of the model parameters
! covar(:,:): normed covariance matrix
! modvar: variance of the model function
! modvar_confid: confindence interval of the variance of the model function
! error: 0->convergence reached, 1->iteration failed
integer :: t,i1,i2
double precision :: chisq,ochisq,alambda,delmax
double precision, allocatable :: alpha(:,:),beta(:),aold(:)
logical :: succmar
character(len=50) :: formatstring
error=1
alambda=-1.0d0
allocate(alpha(na,na))
allocate(beta(na))
allocate(aold(na))
aold=a
do t=1,tmax
call mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
if(succmar) then
write(formatstring,'(a,i0,a)') '(i9,e13.4,',na,'e13.4)'
print formatstring,t,chisq,a
delmax=maxval(abs((a-aold)/a))
if(delmax.lt.acc) then
error=0
alambda=0.0d0
call mrqmin(fct,x,y,sig,ndata,na,a,alpha,covar,beta,chisq,ochisq,alambda,succmar)
modvar=chisq/dble(ndata-na)
modvar_confid(2)=sqrt(2.0d0/dble(ndata-na))
modvar_confid(1)=-modvar_confid(2)
modvar_confid=modvar_confid+1.0d0
do i1=1,na
da(i1)=sqrt(covar(i1,i1))
covar(i1,i1)=1.0d0
do i2=i1+1,na
covar(i1,i2)=covar(i1,i2)/da(i1)
covar(i2,i1)=covar(i1,i2)
end do
end do
print '(/,a,/)','mrqiter converged'
print '(a,f6.3,a,2(f5.3,a))','variance of the model = ',modvar,' (should be between '&
&,modvar_confid(1),' and ',modvar_confid(2),')'
if((modvar.gt.modvar_confid(1)).and.(modvar.lt.modvar_confid(2))) then
print '(a)','good model function'
else
print '(a)','adapt model function!'
end if
exit
else
aold=a
end if
else
print '(i9,a)',t,' iteration was not successful'
end if
end do
if(error.ne.0) then
print '(a)','ERROR in mrqiter: convergence not reached'
end if
deallocate(alpha)
deallocate(beta)
deallocate(aold)
end subroutine mrqiter
subroutine spline(n,x,y,a,b,c,d)
integer, intent(in) :: n
double precision, intent(in) :: x(:),y(:)
double precision, intent(out) :: a(:),b(:),c(:),d(:)
! performs a spline interpolation y=a+b(x-x0)+c(x-x0)**2+d(x-x0)**3 with natural splines
! n: number of datapoints that are interpolated
! x(:): data points that are interpolated
! y(:): function values at the data points x(:)
! a(:): coefficient a
! b(:): coefficient b
! c(:): coefficient c
! d(:): coefficient d
double precision, allocatable :: h(:),r(:),md(:),ld(:),ud(:),inhom(:)
integer :: i1,nm1,nm2
nm1=n-1
nm2=n-2
allocate(h(nm1))
allocate(r(nm1))
allocate(md(nm2))
allocate(ld(nm2))
allocate(ud(nm2))
allocate(inhom(nm2))
do i1=1,nm1
if(x(i1).ge.x(i1+1)) then
print '(a)','ERROR in spline: x values are not strictly monotonically increasing'
return
end if
h(i1)=x(i1+1)-x(i1)
r(i1)=y(i1+1)-y(i1)
end do
do i1=1,n-2
md(i1)=2.0d0*(h(i1)+h(i1+1))
ld(i1)=h(i1)
ud(i1)=h(i1+1)
inhom(i1)=3.0d0*(r(i1+1)/h(i1+1)-r(i1)/h(i1))
end do
call trid(md,ld,ud,inhom,nm2)
c(1)=0.0d0
c(2:nm1)=inhom
do i1=1,nm2
a(i1)=y(i1)
b(i1)=r(i1)/h(i1)-h(i1)*(c(i1+1)+2.0d0*c(i1))/3.0d0
d(i1)=(c(i1+1)-c(i1))/(3.0d0*h(i1))
end do
a(nm1)=y(nm1)
b(nm1)=r(nm1)/h(nm1)-h(nm1)*2.0d0*c(nm1)/3.0d0
d(nm1)=-c(nm1)/(3.0d0*h(nm1))
deallocate(h)
deallocate(r)
deallocate(md)
deallocate(ld)
deallocate(ud)
deallocate(inhom)
end subroutine spline
subroutine splval(n,x,a,b,c,d,xx,yy,dyy,ddyy)
integer, intent(in) :: n
double precision, intent(in) :: x(:),a(:),b(:),c(:),d(:),xx
double precision, intent(out) :: yy
double precision, optional, intent(out) :: dyy,ddyy
! calculates the function value of the spline interpolation y=a+b(x-x0)+c(x-x0)**2+d(x-x0)**3
! n: number of datapoints that were interpolated
! x(:): data points that were interpolated
! a(:): coefficient a
! b(:): coefficient b
! c(:): coefficient c
! d(:): coefficient d
! xx: interpolation point
! yy: function value at the interpolation point xx
! dyy: value of the first derivative of the function at the interpolation point xx
! ddyy: value of the second derivative of the function at the interpolation point xx
integer, parameter :: div_expol=10
integer :: i1,i1l,i1u
double precision :: max_expol,z
max_expol=(x(n)-x(1))/dble(n*div_expol)
if((xx.lt.(x(1)-max_expol)).or.(xx.gt.(x(n)+max_expol))) then
print '(a)','ERROR in splval: xx is not in the interpolation interval'
return
end if
i1u=n+1
i1l=0
do while(i1l.lt.i1u-1)
i1=(i1u+i1l)/2
if(xx.ge.x(i1)) then
i1l=i1
else
i1u=i1
end if
end do
i1=i1l
if(i1.le.0) i1=1
if(i1.ge.n) i1=n-1
z=xx-x(i1)
yy=((d(i1)*z+c(i1))*z+b(i1))*z+a(i1)
if(present(dyy)) dyy=(3.0d0*d(i1)*z+2.0d0*c(i1))*z+b(i1)
if(present(ddyy)) ddyy=6.0d0*d(i1)*z+2.0d0*c(i1)
end subroutine splval
subroutine fft_pow2(a,m,inv)
complex(kind=kind(0.0d0)), intent(inout) :: a(:)
integer, intent(in) :: m,inv
! perform a fast fouriere transformation by the methode of Cooley and Tukey
! a(:): data to transform as input and transformed data as output
! a_new(1)=sum(a_old)
! m: 2**m is the number ot data points
! inv: 1->fouriere transform, else->inverse fouriere transform
integer :: n,nd2,i1,i2,i3,i4,ie1,ie2,ipie
complex(kind=kind(0.0d0)) :: u,w,t
double precision :: ang
n=2**m
nd2=n/2
i2=1
do i1=1,n-1
if(i1.lt.i2) then
t=a(i2)
a(i2)=a(i1)
a(i1)=t
end if
i3=nd2
do while(i3.lt.i2)
i2=i2-i3
i3=i3/2
end do
i2=i2+i3
end do
ie1=1
do i4=1,m
ie2=ie1
ie1=ie1+ie1
u=(1.0d0,0.0d0)
ang=m_pi/dble(ie2)
w=cmplx(cos(ang),-sin(ang),kind(0.0d0))
if(inv.eq.1) w=conjg(w)
do i2=1,ie2
do i1=i2,n,ie1
ipie=i1+ie2
t=a(ipie)*u
a(ipie)=a(i1)-t
a(i1)=a(i1)+t
end do
u=u*w
end do
end do
if(inv.ne.1) then
a=a/dble(n)
end if
end subroutine fft_pow2
subroutine ft_n(a,n,inv)
complex(kind=kind(0.0d0)), intent(inout) :: a(:)
integer, intent(in) :: n,inv
! perform a O(n**2) process fouriere transformation
! a(:): data to transform as input and transformed data as output
! a_new(1)=sum(a_old)
! n: n is the number of data points, it does not have to be a power of 2
! inv: 1->fouriere transform, else->inverse fouriere transform
integer :: i1,i2
complex(kind=kind(0.0d0)), allocatable :: aintern(:),expfac(:)
double precision :: ang
allocate(aintern(n))
aintern=0.0d0
allocate(expfac(n))
ang=2.0d0*m_pi/dble(n)
expfac(1)=(1.0d0,0.0d0)
expfac(2)=cmplx(cos(ang),-sin(ang),kind(0.0d0))
if(inv.eq.1) expfac(2)=conjg(expfac(2))
do i1=3,n
expfac(i1)=expfac(2)*expfac(i1-1)
end do
do i1=1,n
do i2=1,n
aintern(i2)=aintern(i2)+a(i1)*expfac(i1)**(i2-1)
end do
end do
a=aintern
if(inv.ne.1) a=a/dble(n)
end subroutine ft_n
subroutine ft(a,n,inv,info)
complex(kind=kind(0.0d0)), intent(inout) :: a(:)
integer, intent(in) :: n,inv
logical, optional, intent(in) :: info
! perform a fouriere transformation, if n=2**m fast fouriere transformation is done
! a(:): data to transform as input and transformed data as output
! a_new(1)=sum(a_old)
! n: n is the number of data points, it does not have to be a power of 2
! inv: 1->fouriere transform, else->inverse fouriere transform
! info: .true.->information about the used type of fouriere transformation displayed
integer :: m
m=nint(log(dble(n))/log(2.0d0))
if(2**m.eq.n) then
if(present(info).and.info) print '(a)','n is a power of 2 -> use fft'
call fft_pow2(a,m,inv)
else
if(present(info).and.info) print '(a)','n is not a power of 2 -> use slow ft'
call ft_n(a,n,inv)
end if
end subroutine ft
subroutine odeint_ro(derivs,jacobn,ystart,nvar,x1,x2,acc,h1,hmin,nsav,nmax&
&,xfield,yfield,error,dxsav,savstart)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
interface
subroutine jacobn(x,y,dfdx,dfdy)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dfdx(:),dfdy(:,:)
end subroutine jacobn
end interface
double precision, intent(inout) :: ystart(:)
integer, intent(in) :: nvar,nmax
integer, intent(out) :: nsav,error
double precision, intent(in) :: x1,x2,acc,h1,hmin
double precision, intent(out) :: xfield(:),yfield(:,:)
double precision, optional, intent(in) :: dxsav
logical, optional, intent(in) :: savstart
! solve an initial value problem with the method of Rosenbrock
! especially for stiff differential equations
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! jacobn: calculates the second derivatives of the first order differential equation
! x: x value
! y(:): function values
! dfdx(:): derivatives of f = dydx with respect to x
! dfdy(:,:): derivatives of f = dydx with respect to y
! ystart(:): initial values of the initial value problem as input and last values as output
! nvar: number of first order differential equations
! x1: starting point x1
! x2: final point x2
! acc: desired accuracy of one step
! h1: starting step size
! hmin: minimum step size
! nsav: number of saved steps
! nmax: maximum number of steps
! xfield(:): calculated x values
! yfield(:,:): calculated y values, first dimension: different y variables, second dimension: different x values
! error: 0->convergence reached, 1->iteration failed
! dxsav: mimimum distance between two x values that should be saved
! savstart: .true.-> start values are saved in xfield and yfield, else they are not
double precision, parameter :: ctiny=1.0d-30
integer :: nstp,nstpmax
double precision :: h,hdid,hnext,x,xsav,dxsav_intern
double precision, allocatable :: y(:),yscal(:),dydx(:)
error=1
allocate(y(nvar))
allocate(yscal(nvar))
allocate(dydx(nvar))
x=x1
h=sign(h1,x2-x1)
y=ystart
if(present(savstart).and.savstart) then
nsav=1
nstpmax=nmax-1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=ystart
else
nsav=0
nstpmax=nmax
end if
if(.not.present(dxsav)) then
dxsav_intern=0.0d0
else
dxsav_intern=abs(dxsav)
end if
do nstp=1,nstpmax
call derivs(x,y,dydx)
yscal=abs(y)+abs(h*dydx)+ctiny
if((x+h-x2)*(x+h-x1).gt.0.0d0) h=x2-x
call stiff(derivs,jacobn,y,dydx,nvar,x,h,acc,yscal,hdid,hnext)
if(abs(x-xsav).gt.dxsav_intern) then
nsav=nsav+1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=y
end if
if((x-x2)*(x2-x1).ge.0.0d0) then
ystart=y
if(xsav.ne.x) then
nsav=nsav+1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=ystart
end if
error=0
exit
end if
if(abs(hnext).lt.hmin) then
print '(a)','WARNING in odeint_ro: minimum step size underrun'
end if
h=hnext
end do
if(error.ne.0) then
print '(a)','ERROR in odeint_ro: too many steps, final point not reached'
end if
deallocate(y)
deallocate(yscal)
deallocate(dydx)
end subroutine odeint_ro
subroutine stiff(derivs,jacobn,y,dydx,n,x,htry,acc,yscal,hdid,hnext)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
interface
subroutine jacobn(x,y,dfdx,dfdy)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dfdx(:),dfdy(:,:)
end subroutine jacobn
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: n
double precision, intent(in) :: htry,acc,yscal(:)
double precision, intent(out) :: hdid,hnext
! one single step for the method of Rosenbrock
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! jacobn: calculates the second derivatives of the first order differential equation
! x: x value
! y(:): function values
! dfdx(:): derivatives of f = dydx with respect to x
! dfdy(:,:): derivatives of f = dydx with respect to y
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! n: number of first order differential equations
! x: current point as input and next point as output
! htry: suggested step width
! acc: desired accuracy of this step
! yscal: scaled y values needed for the error estimation
! hdid: performed step size
! hnext: suggested step size for the next step
double precision, parameter :: safety=0.9d0,shrnk=0.4d0,pshrnk=-1.0d0/4.0d0&
&,grow=4.0d0,pgrow=-1.0d0/5.0d0,errcongrow=(grow/safety)**(1.0d0/pgrow)&
&,errconshrnk=(shrnk/safety)**(1.0d0/pshrnk),gam=1.0d0/2.0d0&
&,a21=2.0d0&
&,a31=48.0d0/25.0d0,a32=6.0d0/25.0d0&
&,c21=-8.0d0&
&,c31=372.0d0/25.0d0,c32=12.0d0/5.0d0&
&,c41=-112.0d0/125.0d0,c42=-54.0d0/125.0d0,c43=-2.0d0/5.0d0&
&,b1=19.0d0/9.0d0,b2=1.0d0/2.0d0,b3=25.0d0/108.0d0,b4=125.0d0/108.0d0&
&,e1=17.0d0/54.0d0,e2=7.0d0/36.0d0,e3=0.0d0,e4=125.0d0/108.0d0&
&,c1x=1.0d0/2.0d0,c2x=-3.0d0/2.0d0,c3x=121.0d0/50.0d0,c4x=29.0d0/250.0d0&
&,a2x=1.0d0,a3x=3.0d0/5.0d0
integer :: i1
integer, allocatable :: indx(:)
double precision :: det,khad,errmax,h,xsav
double precision, allocatable :: a(:,:),dfdx(:),dfdy(:,:),dysav(:),err(:)&
&,g1(:),g2(:),g3(:),g4(:),ysav(:)
logical :: stepok
allocate(indx(n))
allocate(a(n,n))
allocate(dfdx(n))
allocate(dfdy(n,n))
allocate(dysav(n))
allocate(err(n))
allocate(g1(n))
allocate(g2(n))
allocate(g3(n))
allocate(g4(n))
allocate(ysav(n))
xsav=x
ysav=y
dysav=dydx
call jacobn(xsav,ysav,dfdx,dfdy)
h=htry
stepok=.false.
do while(.not.stepok)
a=-dfdy
do i1=1,n
a(i1,i1)=1.0d0/(gam*h)+a(i1,i1)
end do
call ludcmp(a,n,indx,det,khad)
g1=dysav+h*c1x*dfdx
call lubksb(a,n,indx,g1)
y=ysav+a21*g1
x=xsav+a2x*h
call derivs(x,y,dydx)
g2=dydx+h*c2x*dfdx+c21*g1/h
call lubksb(a,n,indx,g2)
y=ysav+a31*g1+a32*g2
x=xsav+a3x*h
call derivs(x,y,dydx)
g3=dydx+h*c3x*dfdx+(c31*g1+c32*g2)/h
call lubksb(a,n,indx,g3)
g4=dydx+h*c4x*dfdx+(c41*g1+c42*g2+c43*g3)/h
call lubksb(a,n,indx,g4)
y=ysav+b1*g1+b2*g2+b3*g3+b4*g4
err=e1*g1+e2*g2+e3*g3+e4*g4
x=xsav+h
if(x.eq.xsav)print '(a)','ERROR in stiff: stepsize not significant'
errmax=maxval(abs(err/yscal))
errmax=errmax/acc
if(errmax.le.1.0d0) then
hdid=h
if(errmax.gt.errcongrow) then
hnext=safety*h*errmax**pgrow
else
hnext=grow*h
endif
stepok=.true.
else
if(errmax.lt.errconshrnk) then
hnext=safety*h*errmax**pshrnk
else
hnext=shrnk*h
endif
h=hnext
endif
end do
deallocate(indx)
deallocate(a)
deallocate(dfdx)
deallocate(dfdy)
deallocate(dysav)
deallocate(err)
deallocate(g1)
deallocate(g2)
deallocate(g3)
deallocate(g4)
deallocate(ysav)
end subroutine stiff
subroutine odeint_rk(derivs,ystart,nvar,x1,x2,acc,h1,hmin,nsav,nmax,xfield,yfield,error&
&,dxsav,savstart,rkchoose)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: ystart(:)
integer, intent(in) :: nvar,nmax
integer, intent(out) :: nsav,error
double precision, intent(in) :: x1,x2,acc,h1,hmin
double precision, intent(out) :: xfield(:),yfield(:,:)
double precision, optional, intent(in) :: dxsav
logical, optional, intent(in) :: savstart
integer, optional, intent(in) :: rkchoose
! solve an initial value problem with the Runge-Kutta method, different methods can be choosen
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! ystart(:): initial values of the initial value problem as input and last values as output
! nvar: number of first order differential equations
! x1: starting point x1
! x2: final point x2
! acc: desired accuracy of one step
! h1: starting step size
! hmin: minimum step size
! nsav: number of saved steps
! nmax: maximum number of steps
! xfield(:): calculated x values
! yfield(:,:): calculated y values, first dimension: function values, second dimension: different x values
! error: 0->convergence reached, 1->iteration failed
! dxsav: mimimum distance between two x values that should be saved
! savstart: .true.-> start values are saved in xfield and yfield, else they are not
! rkchoose: choose the desired Runge-Kutta method
! 1->Runge-Kutta-Fehlberg
! 2->Cash-Karp
! 3->Dormand-Prince
! 4->classical Runge-Kutta 4th order
! else or not present->Runge-Kutta-Fehlberg
double precision, parameter :: ctiny=1.0d-30
integer :: nstp,nstpmax,rk_intern
double precision :: h,hnext,x,xsav,dxsav_intern
double precision, allocatable :: y(:),yscal(:),dydx(:)
error=1
if(.not.present(rkchoose)) then
rk_intern=1
else if(rkchoose.eq.2) then
rk_intern=2
else if(rkchoose.eq.3) then
rk_intern=3
else if(rkchoose.eq.4) then
rk_intern=4
else
rk_intern=1
end if
allocate(y(nvar))
allocate(yscal(nvar))
allocate(dydx(nvar))
x=x1
h=sign(h1,x2-x1)
y=ystart
if(present(savstart).and.savstart) then
nsav=1
nstpmax=nmax-1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=ystart
else
nsav=0
nstpmax=nmax
end if
if(.not.present(dxsav)) then
dxsav_intern=0.0d0
else
dxsav_intern=abs(dxsav)
end if
do nstp=1,nstpmax
call derivs(x,y,dydx)
yscal=abs(y)+abs(h*dydx)+ctiny
if((x+h-x2)*(x+h-x1).gt.0.0d0) h=x2-x
call rkqc(derivs,y,dydx,nvar,x,h,acc,yscal,hnext,rk_intern)
if(abs(x-xsav).gt.dxsav_intern) then
nsav=nsav+1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=y
end if
if((x-x2)*(x2-x1).ge.0.0d0) then
ystart=y
if(xsav.ne.x) then
nsav=nsav+1
xfield(nsav)=x
xsav=x
yfield(:,nsav)=ystart
end if
error=0
exit
end if
if(abs(hnext).lt.hmin) then
print '(a,3e16.6)','WARNING in odeint_rk: minimum step size underrun',h,hnext,x-x2
end if
h=hnext
end do
if(error.ne.0) then
print '(a)','ERROR in odeint_rk: too many steps, final point not reached'
end if
deallocate(y)
deallocate(yscal)
deallocate(dydx)
end subroutine odeint_rk
subroutine rkqc(derivs,y,dydx,nvar,x,htry,acc,yscal,hnext,rkchoose)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: htry,acc,yscal(:)
double precision, intent(out) :: hnext
integer, intent(in) :: rkchoose
! one single step of a Runge-Kutta method, different methods can be choosen
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! htry: suggested step width
! acc: desired accuracy of this step
! yscal: scaled y values needed for the error estimation
! hnext: suggested step size for the next step
! rkchoose: which Runge-Kutta method is chosen
! 1->Runge-Kutta-Fehlberg
! 2->Cash-Karp
! 3->Dormand-Prince
! 4->classical Runge-Kutta 4th order
double precision, parameter :: safety=0.9d0,shrnk=0.4d0,pshrnk=-1.0d0/4.0d0,grow=4.0d0,pgrow=-1.0d0/5.0d0
double precision, parameter :: errcongrow=(grow/safety)**(1.0d0/pgrow),errconshrnk=(shrnk/safety)**(1.0d0/pshrnk)
double precision :: h,xsav,errmax
double precision, allocatable :: ysav(:),dysav(:),yerr(:)
logical :: stepok
xsav=x
allocate(ysav(nvar))
allocate(dysav(nvar))
allocate(yerr(nvar))
ysav=y
dysav=dydx
h=htry
stepok=.false.
do while(.not.stepok)
select case(rkchoose)
case(1)
call rkf(derivs,y,dydx,nvar,x,h,yerr)
case(2)
call rkck(derivs,y,dydx,nvar,x,h,yerr)
case(3)
call rkdp(derivs,y,dydx,nvar,x,h,yerr)
case(4)
call rkc4err(derivs,y,dydx,nvar,x,h,yerr)
end select
errmax=maxval(abs(yerr/yscal))
errmax=errmax/acc
if(errmax.gt.1.0d0) then
if(errmax.gt.errconshrnk) then
hnext=shrnk*h
else
hnext=safety*h*errmax**pshrnk
end if
h=hnext
x=xsav
y=ysav
dydx=dysav
else
stepok=.true.
if(errmax.lt.errcongrow) then
hnext=grow*h
else
hnext=safety*h*errmax**pgrow
end if
end if
end do
deallocate(ysav)
deallocate(dysav)
deallocate(yerr)
end subroutine rkqc
subroutine rkf(derivs,y,dydx,nvar,x,h,yerr)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: h
double precision, intent(out) :: yerr(:)
! one single step of a Runge-Kutta methode in the Fehlberg version of error calculation
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! h: step size
! yerr(:): calculated estimated error
double precision :: xsav
double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ytemp(:)
double precision, parameter :: a2=0.25d0,a3=3.0d0/8.0d0,a4=12.0d0/13.0d0,a5=1.0d0,a6=0.5d0&
&,b21=0.25d0&
&,b31=3.0d0/32.0d0,b32=9.0d0/32.0d0&
&,b41=1932.0d0/2197.0d0,b42=-7200.0d0/2197.0d0,b43=7296.0d0/2197.0d0&
&,b51=439.0d0/216.0d0,b52=-8.0d0,b53=3680.0d0/513.0d0,b54=-845.0d0/4104.0d0&
&,b61=-8.0d0/27.0d0,b62=2.0d0,b63=-3544.0d0/2565.0d0,b64=1859.0d0/4104.0d0,b65=-11.0d0/40.0d0&
&,c1=16.0d0/135.0d0,c3=6656.0d0/12825.0d0,c4=28561.0d0/56430.0d0,c5=-9.0d0/50.0d0,c6=2.0d0/55.0d0&
&,dc1=c1-25.0d0/216.0d0,dc3=c3-1408.0d0/2565.0d0,dc4=c4-2197.0d0/4104.0d0,dc5=c5+1.0d0/5.0d0,dc6=c6
allocate(ak2(nvar))
allocate(ak3(nvar))
allocate(ak4(nvar))
allocate(ak5(nvar))
allocate(ak6(nvar))
allocate(ytemp(nvar))
xsav=x
ytemp=y+b21*h*dydx
x=xsav+a2*h
call derivs(x,ytemp,ak2)
ytemp=y+h*(b31*dydx+b32*ak2)
x=xsav+a3*h
call derivs(x,ytemp,ak3)
ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
x=xsav+a4*h
call derivs(x,ytemp,ak4)
ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
x=xsav+a5*h
call derivs(x,ytemp,ak5)
ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
x=xsav+a6*h
call derivs(x,ytemp,ak6)
y=y+h*(c1*dydx+c3*ak3+c4*ak4+c5*ak5+c6*ak6)
yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6)
x=xsav+h
deallocate(ak2)
deallocate(ak3)
deallocate(ak4)
deallocate(ak5)
deallocate(ak6)
deallocate(ytemp)
end subroutine rkf
subroutine rkck(derivs,y,dydx,nvar,x,h,yerr)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: h
double precision, intent(out) :: yerr(:)
! one single step of a Runge-Kutta methode in the Cash-Karp version of error calculation
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! h: step size
! yerr(:): calculated estimated error
double precision :: xsav
double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ytemp(:)
double precision, parameter :: a2=0.2d0,a3=0.3d0,a4=0.6d0,a5=1.0d0,a6=0.875d0&
&,b21=0.2d0&
&,b31=3.0d0/40.0d0,b32=9.0d0/40.0d0&
&,b41=0.3d0,b42=-0.9d0,b43=1.2d0&
&,b51=-11.0d0/54.0d0,b52=2.50d0,b53=-70.0d0/27.0d0,b54=35.0d0/27.0d0&
&,b61=1631.0d0/55296.0d0,b62=175.0d0/512.0d0,b63=575.0d0/13824.0d0,b64=44275.0d0/110592.0d0,b65=253.0d0/4096.0d0&
&,c1=37.0d0/378.0d0,c3=250.0d0/621.0d0,c4=125.0d0/594.0d0,c6=512.0d0/1771.0d0&
&,dc1=c1-2825.0d0/27648.0d0,dc3=c3-18575.0d0/48384.0d0,dc4=c4-13525.0d0/55296.0d0,dc5=-277.0d0/14336.0d0,dc6=c6-0.25d0
allocate(ak2(nvar))
allocate(ak3(nvar))
allocate(ak4(nvar))
allocate(ak5(nvar))
allocate(ak6(nvar))
allocate(ytemp(nvar))
xsav=x
ytemp=y+b21*h*dydx
x=xsav+a2*h
call derivs(x,ytemp,ak2)
ytemp=y+h*(b31*dydx+b32*ak2)
x=xsav+a3*h
call derivs(x,ytemp,ak3)
ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
x=xsav+a4*h
call derivs(x,ytemp,ak4)
ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
x=xsav+a5*h
call derivs(x,ytemp,ak5)
ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
x=xsav+a6*h
call derivs(x,ytemp,ak6)
y=y+h*(c1*dydx+c3*ak3+c4*ak4+c6*ak6)
yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6)
x=xsav+h
deallocate(ak2)
deallocate(ak3)
deallocate(ak4)
deallocate(ak5)
deallocate(ak6)
deallocate(ytemp)
end subroutine rkck
subroutine rkdp(derivs,y,dydx,nvar,x,h,yerr)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: h
double precision, intent(out) :: yerr(:)
! one single step of a Runge-Kutta methode in the Dormand-Prince version of error calculation
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! h: step size
! yerr(:): calculated estimated error
double precision :: xsav
double precision, allocatable :: ak2(:),ak3(:),ak4(:),ak5(:),ak6(:),ak7(:),ytemp(:)
double precision, parameter :: a2=0.2d0,a3=0.3d0,a4=0.8d0,a5=8.0d0/9.0d0,a6=1.0d0,a7=1.0d0&
&,b21=0.2d0&
&,b31=3.0d0/40.0d0,b32=9.0d0/40.0d0&
&,b41=44.0d0/45.0d0,b42=-56.0d0/15.0d0,b43=32.0d0/9.0d0&
&,b51=19372.0d0/6561.0d0,b52=-25360.0d0/2187.0d0,b53=64448.0d0/6561.0d0,b54=-212.0d0/729.0d0&
&,b61=9017.0d0/3168.0d0,b62=-355.0d0/33.0d0,b63=46732.0d0/5247.0d0,b64=49.0d0/176.0d0,b65=-5103.0d0/18656.0d0&
&,b71=35.0d0/384.0d0,b73=500.0d0/1113.0d0,b74=125.0d0/192.0d0,b75=-2187.0d0/6784.0d0,b76=11.0d0/84.0d0&
&,c1=5179.0d0/57600.0d0,c3=7571.0d0/16695.0d0,c4=393.0d0/640.0d0,c5=-92097.0d0/339200.0d0,c6=187.0d0/2100.0d0&
&,c7=1.0d0/40.0d0&
&,dc1=c1-35.0d0/384.0d0,dc3=c3-500.0d0/1113.0d0,dc4=c4-125.0d0/192.0d0,dc5=c5+2187.0d0/6784.0d0&
&,dc6=c6-11.0d0/84.0d0,dc7=c7
allocate(ak2(nvar))
allocate(ak3(nvar))
allocate(ak4(nvar))
allocate(ak5(nvar))
allocate(ak6(nvar))
allocate(ak7(nvar))
allocate(ytemp(nvar))
xsav=x
ytemp=y+b21*h*dydx
x=xsav+a2*h
call derivs(x,ytemp,ak2)
ytemp=y+h*(b31*dydx+b32*ak2)
x=xsav+a3*h
call derivs(x,ytemp,ak3)
ytemp=y+h*(b41*dydx+b42*ak2+b43*ak3)
x=xsav+a4*h
call derivs(x,ytemp,ak4)
ytemp=y+h*(b51*dydx+b52*ak2+b53*ak3+b54*ak4)
x=xsav+a5*h
call derivs(x,ytemp,ak5)
ytemp=y+h*(b61*dydx+b62*ak2+b63*ak3+b64*ak4+b65*ak5)
x=xsav+a6*h
call derivs(x,ytemp,ak6)
ytemp=y+h*(b71*dydx+b73*ak3+b74*ak4+b75*ak5+b76*ak6)
x=xsav+a7*h
call derivs(x,ytemp,ak7)
y=y+h*(c1*dydx+c3*ak3+c4*ak4+c5*ak5+c6*ak6+c7*ak7)
yerr=h*(dc1*dydx+dc3*ak3+dc4*ak4+dc5*ak5+dc6*ak6+dc7*ak7)
x=xsav+h
deallocate(ak2)
deallocate(ak3)
deallocate(ak4)
deallocate(ak5)
deallocate(ak6)
deallocate(ak7)
deallocate(ytemp)
end subroutine rkdp
subroutine rkc4err(derivs,y,dydx,nvar,x,h,yerr)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: h
double precision, intent(out) :: yerr(:)
! calssical Runge-Kutte method 4th order, error calculation by h/2 method
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! h: step size
! yerr(:): calculated estimated error
double precision :: hhalf,xsav
double precision, allocatable :: ysav(:),dysav(:)
xsav=x
allocate(ysav(nvar))
allocate(dysav(nvar))
ysav=y
dysav=dydx
hhalf=0.5d0*h
call rkc4(derivs,y,dydx,nvar,x,hhalf)
call derivs(x,y,dydx)
call rkc4(derivs,y,dydx,nvar,x,hhalf)
if(x.eq.xsav) then
print '(a)','WARNING in rkc4err: half stepsize not significant'
x=xsav+h
if(x.eq.xsav) then
print '(a)','ERROR in rkc4err: stepsize not significant'
end if
end if
call rkc4(derivs,ysav,dysav,nvar,xsav,h)
yerr=y-ysav
deallocate(ysav)
deallocate(dysav)
end subroutine rkc4err
subroutine rkc4(derivs,y,dydx,nvar,x,h)
interface
subroutine derivs(x,y,dydx)
double precision, intent(in) :: x,y(:)
double precision, intent(out) :: dydx(:)
end subroutine derivs
end interface
double precision, intent(inout) :: y(:),dydx(:),x
integer, intent(in) :: nvar
double precision, intent(in) :: h
! classical Runge-Kutte method 4th order
! derivs: calculates the derivatives of the first order differential equation
! x: x value
! y(:): function values
! dydx(:): derivatives
! y(:): initial values of the initial value problem as input and next values as output
! dydx(:): derivatives according to the differntial equation
! nvar: number of first order differential equations
! x: current point as input and next point as output
! h: step size
double precision, allocatable :: ysav(:),yki(:)
allocate(ysav(nvar))
allocate(yki(nvar))
ysav=y
y=y+h/6.0d0*dydx
yki=ysav+h/2.0d0*dydx
x=x+h/2.0d0
call derivs(x,yki,dydx)
y=y+h/3.0d0*dydx
yki=ysav+h/2.0d0*dydx
call derivs(x,yki,dydx)
y=y+h/3.0d0*dydx
yki=ysav+h*dydx
x=x+h/2.0d0
call derivs(x,yki,dydx)
y=y+h/6.0d0*dydx
deallocate(ysav)
deallocate(yki)
end subroutine rkc4
real function getnan()
! returns the indefinite value nan (not a number) for every data type
! getnan: returns nan
real :: dum
dum=0.0
getnan=dum/dum
end function getnan
real function getinf()
! returns the indefinite value nan (not a number) for every data type
! getnan: returns nan
real :: dum
dum=0.0
getinf=1.0/dum
end function getinf
integer function factorial(n)
integer, intent(in) :: n
! calculate n! (n factorial)
! n: integer number as input
! factorial: n! as output
integer :: i1
i1=1
factorial=product((/(i1,i1=1,n)/))
end function factorial
integer function nCr(n,k)
integer, intent(in) :: n,k
! calculate the binomial coefficient (n,k)=n!/(k!(n-k)!)
! n: integer n as input
! k: integer k as input
! nCr: binomial coefficient as output
integer :: i1
nCr=1
if(k.gt.n-k) then
do i1=1,n-k
nCr=nCr*(i1+k)/i1
end do
else
do i1=1,k
nCr=nCr*(i1+n-k)/i1
end do
end if
end function nCr
subroutine randinit(seed)
integer, intent(inout) :: seed
! initialize the random number generator rand() with an arbitrary seed derived out of the system time
! seed: arbitrary integer number out of which the seed is generated, odd number recommended
! as output the used seed is handed over
integer :: seed1
real :: seedout
call system_clock(count=seed1)
seed1=seed1+time()
seed=seed1*2+seed
if(seed.lt.0) seed=-seed
seedout=rand(seed)
end subroutine randinit
subroutine randinit_const(seed)
integer, intent(in) :: seed
! initialize the random number generator rand() with a certain seed
! seed: seed used for the initialization
real :: seedout
seedout=rand(seed)
end subroutine randinit_const
real function rand_0()
! create a random number from a uniform distribution in the interval (0,1]
rand_0=1.0-rand()
end function rand_0
real function randn()
! creates a gaussian distributed random number with mean 0 and sigma 1 (normal distribution)
! code from http://www.netlib.org/random/ downloaded on 06.10.2014
! part of the original comment from the downloaded code:
! Some of the code is adapted from Dagpunar's book:
! Dagpunar, J. 'Principles of random variate generation'
! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9
! Version 1.13, 2 October 2000
! Changes from version 1.01
! 1. The random_order, random_Poisson & random_binomial routines have been
! replaced with more efficient routines.
! 2. A routine, seed_random_number, has been added to seed the uniform random
! number generator. This requires input of the required number of seeds
! for the particular compiler from a specified I/O unit such as a keyboard.
! 3. Made compatible with Lahey's ELF90.
! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
! 5. INTENT for array f corrected in random_mvnorm.
!
! Author: Alan Miller
! CSIRO Division of Mathematical & Information Sciences
! Private Bag 10, Clayton South MDC
! Clayton 3169, Victoria, Australia
! Phone: (+61) 3 9545-8016 Fax: (+61) 3 9545-8080
! e-mail: amiller @ bigpond.net.au
!
! Adapted from the following Fortran 77 code
! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
!
! The function random_normal() returns a normally distributed pseudo-random
! number with zero mean and unit variance.
!
! The algorithm uses the ratio of uniforms method of A.J. Kinderman
! and J.F. Monahan augmented with quadratic bounding curves.
real, parameter :: s=0.449871,t=-0.386595,a=0.19600,b=0.25472,r1=0.27597,r2=0.27846
real :: u,v,x,y,q
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
do
u=rand()
v=rand()
v=1.7156*(v-0.5)
! Evaluate the quadratic form
x=u-s
y=abs(v)-t
q=x**2+y*(a*y-b*x)
! Accept P if inside inner ellipse
if(q.lt.r1) exit
! Reject P if outside outer ellipse
if(q.gt.r2) cycle
! Reject P if outside acceptance region or except if inside acceptance region
if(v**2.lt.-4.0*LOG(u)*u**2) exit
end do
! Return ratio of P's coordinates as the normal deviate
randn=v/u
end function randn
real function randfermihot()
! get a random number of the x > 0 part of 1/(1+exp(x)) by inversion method
real, parameter :: log2=log(2.0)
randfermihot=-log(exp(log2*(1.0-rand()))-1.0)
end function randfermihot
subroutine random_init(seed)
integer, intent(in) :: seed
! initialize the random number generator random_number() with an arbitrary seed derived out of the system time
! seed: arbitrary integer number out of which the seed is generated, odd number recommended
integer :: i1,n,clock
integer, allocatable :: seed1(:)
call random_seed(size=n)
allocate(seed1(n))
call system_clock(count=clock)
seed1=seed+clock+37*(/(i1-1,i1=1,n)/)
call random_seed(put=seed1)
deallocate(seed1)
end subroutine random_init
subroutine random_init_const(seed)
integer, intent(in) :: seed
! initialize the random number generator random_number() with a certain seed
! seed: seed used for the initialization
integer :: i1,n
integer, allocatable :: seed1(:)
call random_seed(size=n)
allocate(seed1(n))
seed1=seed+37*(/(i1-1,i1=1,n)/)
call random_seed(put=seed1)
deallocate(seed1)
end subroutine random_init_const
real function random_0()
! create a random number from a uniform distribution in the interval (0,1]
call random_number(random_0)
random_0=1.0-random_0
end function random_0
real function random_1()
! create a random number from a uniform distribution in the interval [0,1)
call random_number(random_1)
end function random_1
double precision function random_exp(rate)
! create a random number from an exponential distribution with specified rate parameter
! in the interval (0,\infty)
! rate: rate parameter p(random_exp)=1/rate*exp(-rate*random_exp)
double precision :: rate
real :: random_01
do
call random_number(random_01)
if(random_01.ne.0.0) exit
end do
random_exp=-log(real(random_01,8))/rate
end function random_exp
real function random_n()
! creates a gaussian distributed random number with mean 0 and sigma 1 (normal distribution)
! code from http://www.netlib.org/random/ downloaded on 06.10.2014
! part of the original comment from the downloaded code:
! Some of the code is adapted from Dagpunar's book:
! Dagpunar, J. 'Principles of random variate generation'
! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9
! Version 1.13, 2 October 2000
! Changes from version 1.01
! 1. The random_order, random_Poisson & random_binomial routines have been
! replaced with more efficient routines.
! 2. A routine, seed_random_number, has been added to seed the uniform random
! number generator. This requires input of the required number of seeds
! for the particular compiler from a specified I/O unit such as a keyboard.
! 3. Made compatible with Lahey's ELF90.
! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.
! 5. INTENT for array f corrected in random_mvnorm.
!
! Author: Alan Miller
! CSIRO Division of Mathematical & Information Sciences
! Private Bag 10, Clayton South MDC
! Clayton 3169, Victoria, Australia
! Phone: (+61) 3 9545-8016 Fax: (+61) 3 9545-8080
! e-mail: amiller @ bigpond.net.au
!
! Adapted from the following Fortran 77 code
! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
!
! The function random_normal() returns a normally distributed pseudo-random
! number with zero mean and unit variance.
!
! The algorithm uses the ratio of uniforms method of A.J. Kinderman
! and J.F. Monahan augmented with quadratic bounding curves.
real, parameter :: s=0.449871,t=-0.386595,a=0.19600,b=0.25472,r1=0.27597,r2=0.27846
real :: u,v,x,y,q
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
do
call random_number(u)
call random_number(v)
v=1.7156*(v-0.5)
! Evaluate the quadratic form
x=u-s
y=abs(v)-t
q=x**2+y*(a*y-b*x)
! Accept P if inside inner ellipse
if(q.lt.r1) exit
! Reject P if outside outer ellipse
if(q.gt.r2) cycle
! Reject P if outside acceptance region or except if inside acceptance region
if(v**2.lt.-4.0*LOG(u)*u**2) exit
end do
! Return ratio of P's coordinates as the normal deviate
random_n=v/u
end function random_n
real function random_fermihot()
! get a random number of the x > 0 part of 1/(1+exp(x)) by inversion method
real, parameter :: log2=log(2.0)
random_fermihot=-log(exp(log2*(random_0()))-1.0)
end function random_fermihot
subroutine linfit(x,y,n,k,d)
double precision, intent(in) :: x(:),y(:)
integer, intent(in) :: n
double precision, optional, intent(out) :: k,d
! perform a linear fit of the given data y=kx+d
! x(:): x data for the fit
! y(:): y data for the fit
! n: number of data points
! k: gradient of the linear fit
! d: shift of the linear fit
double precision :: S_x,S_y,S_xx,S_xy,S_yy,Delta
S_x=sum(x)
S_y=sum(y)
S_xx=sum(x*x)
S_xy=sum(x*y)
S_yy=sum(y*y)
Delta=dble(n)*S_xx-S_x*S_x
if(present(k)) k=(dble(n)*S_xy-S_x*S_y)/Delta
if(present(d)) d=(S_xx*S_y-S_xy*S_x)/Delta
end subroutine linfit
subroutine linfit_error(x,y,sig,n,a,da,covar,modvar,modvar_confid)
double precision, intent(in) :: x(:),y(:),sig(:)
integer, intent(in) :: n
double precision, intent(out) :: a(2)
double precision, optional, intent(out) :: da(2),covar(2,2),modvar,modvar_confid(2)
! perform a linear fit of the given data y=a(1)x+a(2)
! x(:): x data
! y(:): y data
! sig(:): error of the y data
! n: number of data points
! a(2): model parameters
! da(2): error of the model parameters
! covar(2,2): normed covariance matrix
! modvar: variance of the model function
! modvar_confid: confindence interval of the variance of the model function
double precision :: ysig,ysig_2,x_m,xx_m,xy_m,y_m,C_xx,C_xy,chisq
integer :: i1
ysig_2=1/sum(1/sig**2)
ysig=sqrt(ysig_2)
x_m=sum(x/sig**2)*ysig_2
xx_m=sum(x**2/sig**2)*ysig_2
C_xx=xx_m-x_m**2
y_m=sum(y/sig**2)*ysig_2
xy_m=sum(x*y/sig**2)*ysig_2
C_xy=xy_m-x_m*y_m
a(1)=C_xy/C_xx
a(2)=y_m-x_m*a(1)
if(present(da)) then
da(1)=ysig_2/C_xx
da(2)=sqrt(ysig_2*xx_m/C_xx)
end if
if(present(covar)) then
covar=1.0d0
covar(1,2)=-ysig_2*x_m/C_xx/sqrt(da(1)*da(2))
covar(2,1)=covar(1,2)
end if
if(present(modvar)) then
chisq=0.0d0
do i1=1,n
chisq=chisq+((y(i1)-a(1)*x(i1)-a(2))/sig(i1))**2
end do
modvar=chisq/dble(n-2)
end if
if(present(modvar_confid)) then
modvar_confid(2)=sqrt(2.0d0/dble(n-2))
modvar_confid(1)=-modvar_confid(2)
modvar_confid=modvar_confid+1.0d0
end if
end subroutine linfit_error
subroutine norm_vector(vector,n,norm)
double precision, intent(inout) :: vector(:)
integer, intent(in) :: n
double precision, optional, intent(out) :: norm
! norm the vector
! vector(:): vector to norm
! n: size of the vector
! norm: calculated norm
integer :: i1
double precision :: norm_intern
norm_intern=0.0d0
do i1=1,n
norm_intern=norm_intern+vector(i1)*vector(i1)
end do
norm_intern=sqrt(norm_intern)
vector=vector/norm_intern
if(present(norm)) norm=norm_intern
end subroutine norm_vector
subroutine lanczos(sparse_mat,index_mat,n_mat,n_tot,n_eigval,eigval,acc,error&
&,eigvec,skip_diag,max_steps,ghost_state_level,startvec)
double precision, intent(in) :: sparse_mat(:),acc
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
integer, intent(inout) :: n_eigval
integer, intent(out) :: error
double precision, intent(out) :: eigval(:)
double precision, optional, intent(out) :: eigvec(:,:)
integer, optional, intent(in) :: skip_diag,max_steps
double precision, optional, intent(in) :: ghost_state_level,startvec(:)
! find maximum absolute value eigenvalues of a huge symmetric matrix by the method of Lanczos
! sparse_mat(:): matrix to diagonalize in spares form, hermiticity considered
! index_mat(:,1): row indices of the spares matrix
! index_mat(:,2): column indices of the spares matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! n_eigval: number of eigenvalues that should be converged and calculated as input
! and number of calculated eigenvalues as output
! eigval(:): calculated eigenvalues
! acc: desired accuracy
! error: -1: convergence not reached in max_steps Lanczos steps
! 1: possibly a ghost state appeared
! 2: possibly a ghost state appeared in the first step
! 0: convergence reached
! eigvec(:,:): calculated eigenvectors, row: basis dimension, column: different energies
! skip_diag: skipping Lanczos steps between two convergence checks, > 0
! max_steps: maximum number of Lanczos steps
! ghost_state_level: minimum level for the coupling constant in the Lanczos matrix before subspace is diagonalized
! startvec(:): start vector for Lanczos method
integer :: i1,lanczos_step,skip_diag_intern,max_steps_intern,last_diag,seed
double precision, allocatable :: x_n_1(:),x_n_2(:),mat_diag(:),mat_coupling(:)
double precision, allocatable :: eigval_lanczos(:),eigvec_lanczos(:,:)
double precision :: ghost_state_level_intern
if(present(skip_diag)) then
skip_diag_intern=skip_diag
else
skip_diag_intern=4
end if
if(present(max_steps)) then
max_steps_intern=max_steps
else
max_steps_intern=20
end if
if(present(ghost_state_level)) then
ghost_state_level_intern=ghost_state_level
else
ghost_state_level_intern=1.0d-9
end if
error=-1
print '(/,a)','calculating Lanczos eigenvalues'
allocate(x_n_1(n_mat))
if(present(startvec)) then
x_n_1=startvec
else
seed=34123451
call randinit(seed)
call get_lanczos_vec_x_0(x_n_1,n_mat)
end if
allocate(x_n_2(n_mat))
allocate(mat_diag(max_steps_intern+1))
allocate(mat_coupling(max_steps_intern+1))
mat_coupling(1)=0.0d0
call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
if(mat_coupling(2).lt.ghost_state_level_intern) error=2
lanczos_step=1
last_diag=lanczos_step/skip_diag_intern
allocate(eigval_lanczos(n_eigval))
allocate(eigvec_lanczos(max_steps_intern+1,n_eigval))
eigval_lanczos=0.0d0
do while((lanczos_step.lt.max_steps_intern).and.(error.lt.0))
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
error=1
exit
end if
if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
if(error.ge.0) exit
last_diag=lanczos_step/skip_diag_intern
end if
if(lanczos_step.ge.max_steps_intern) exit
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
error=1
exit
end if
if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
if(error.ge.0) exit
last_diag=lanczos_step/skip_diag_intern
end if
end do
if(error.eq.0) then
print '(/,a)','lanczos converged'
else if(error.eq.1) then
print '(/,a)','WARNING in lanczos: possible ghost state appeared'
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step-1,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
error=1
if(lanczos_step.lt.n_eigval) then
eigval_lanczos(lanczos_step+1:n_eigval)=getnan()
if(present(eigvec)) then
eigvec(:,lanczos_step+1:n_eigval)=getnan()
end if
n_eigval=lanczos_step
end if
lanczos_step=lanczos_step-1
else if(error.eq.2) then
print '(/,a)','WARNING in lanczos: possible ghost state appeared in the first step'
eigval_lanczos(1)=mat_diag(1)
if(n_eigval.gt.1) then
eigval_lanczos(2:n_eigval)=getnan()
if(present(eigvec)) then
eigvec(:,2:n_eigval)=getnan()
end if
n_eigval=1
end if
if(present(eigvec)) then
eigvec(:,1)=x_n_1
end if
else
print '(/,a,i0,a)','ERROR in lanczos: no convergence after ',max_steps_intern,' iterations'
end if
eigval=eigval_lanczos
deallocate(eigval_lanczos)
if(present(eigvec).and.((error.eq.0).or.(error.eq.1))) then
print '(/,a)','calculating Lanczos eigenvectors'
if(present(startvec)) then
x_n_1=startvec
else
call randinit_const(seed)
call get_lanczos_vec_x_0(x_n_1,n_mat)
end if
do i1=1,n_eigval
eigvec(:,i1)=x_n_1*eigvec_lanczos(1,i1)
end do
call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(2,i1)
end do
max_steps_intern=lanczos_step
lanczos_step=1
do while(lanczos_step.lt.max_steps_intern)
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_1*eigvec_lanczos(lanczos_step+1,i1)
end do
if(lanczos_step.ge.max_steps_intern) exit
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(lanczos_step+1,i1)
end do
end do
end if
deallocate(eigvec_lanczos)
deallocate(x_n_1)
deallocate(x_n_2)
deallocate(mat_diag)
deallocate(mat_coupling)
contains
subroutine get_lanczos_vec_x_0(x_0,n)
double precision, intent(out) :: x_0(:)
integer, intent(in) :: n
! initialize starting vector for Lanczos method with random numbers
! x_0(:): starting vector
! n: size of the vector
integer :: i1
do i1=1,n
x_0(i1)=rand()
end do
call norm_vector(x_0,n)
end subroutine get_lanczos_vec_x_0
subroutine get_lanczos_vec_x_1(x_0,sparse_mat,index_mat,n_mat,n_tot,x_1,eigval_0,eigval_1,coupling_1)
double precision, intent(in) :: x_0(:),sparse_mat(:)
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
double precision, intent(out) :: x_1(:),eigval_0,eigval_1,coupling_1
! perform the first Lanczos step
! x_0(:): starting vector
! sparse_mat(:): matrix of the eigenvalue problem
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! x_1(:): calculated first Lanczos vector
! eigval_0: first diagonal term in the Lanczos matrix
! eigval_0: second diagonal term in the Lanczos matrix
! coupling_1: first offdiagonal term in the Lanczos matrix
integer :: i1
call get_expval(sparse_mat,index_mat,n_tot,x_0,eigval_0)
x_1=-eigval_0*x_0
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,1))
else
x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,2))
x_1(index_mat(i1,2))=x_1(index_mat(i1,2))+sparse_mat(i1)*x_0(index_mat(i1,1))
end if
end do
call norm_vector(x_1,n_mat,coupling_1)
call get_expval(sparse_mat,index_mat,n_tot,x_1,eigval_1)
end subroutine get_lanczos_vec_x_1
subroutine get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n,x_n1,eigval_n,coupling_n&
&,eigval_n1,coupling_n1)
double precision, intent(in) :: sparse_mat(:),eigval_n,coupling_n
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
double precision, intent(inout) :: x_n(:),x_n1(:)
double precision, intent(out) :: eigval_n1,coupling_n1
! perform a general Lanczos step
! sparse_mat(:): matrix of the eigenvalue problem
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! x_n(:): vector for the n-th interation
! x_n1(:): vector for the (n-1)-th iteration as input and (n+1)-th interation as output
! eigval_n: n-th diagonal term in the Lanczos matrix
! coupling_n: n-th offdiagonal term in the Lanczos matrix
! eigval_n1: (n+1)-th diagonal term in the Lanczos matrix
! coupling_n1: (n+1)-th offdiagonal term in the Lanczos matrix
integer :: i1
x_n1=-eigval_n*x_n-coupling_n*x_n1
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,1))
else
x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,2))
x_n1(index_mat(i1,2))=x_n1(index_mat(i1,2))+sparse_mat(i1)*x_n(index_mat(i1,1))
end if
end do
call norm_vector(x_n1,n_mat,coupling_n1)
call get_expval(sparse_mat,index_mat,n_tot,x_n1,eigval_n1)
end subroutine get_lanczos_vec_x_n_plus_1
subroutine get_expval(sparse_mat,index_mat,n_tot,vector,expval)
double precision, intent(in) :: sparse_mat(:),vector(:)
integer, intent(in) :: index_mat(:,:),n_tot
double precision, intent(out) :: expval
! calculate the expectation value of a vector and the sparse matrix (hermiticity used)
! sparse_mat(:): matrix to calculate the expectation value
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! vector(:): vector to calculate the expectation value
! expval: expectation value
integer :: i1
expval=0.0d0
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
expval=expval+vector(index_mat(i1,1))*sparse_mat(i1)*vector(index_mat(i1,2))
else
expval=expval+2.0d0*vector(index_mat(i1,1))*sparse_mat(i1)*vector(index_mat(i1,2))
end if
end do
end subroutine get_expval
subroutine diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step&
&,n_eigval,eigval_lanczos,eigvec_lanczos,error,acc)
double precision, intent(in) :: mat_diag(:),mat_coupling(:),acc
integer, intent(in) :: lanczos_step,n_eigval
double precision, intent(inout) :: eigval_lanczos(:)
double precision, intent(out) :: eigvec_lanczos(:,:)
integer, intent(out) :: error
! calculate the lowest eigenvalues in the spaned subspace and check if they converged
! mat_diag(:): main diagonal of the Lanczos matrix
! mat_coupling(:): first minor diagonal of the Lanczos matrix
! lanczos_step: current number of executed Lanczos steps
! n_eigval: number of eigenvalues that should be calculated
! eigval_lanczos(:): calculated lowest eigenvalues of the Lanczos matrix
! eigvec_lanczos(:,:): eigenvectors of the calculated lowest eigenvalues of the Lanczos matrix
! error: unchanged: no convergence for this check
! 0: convergence reached
! acc: desired accuracy
double precision, allocatable :: lanczos_diag(:),lanczos_offdiag(:),eigvec(:,:)
double precision :: error_eigval
integer :: i1,error_diag
allocate(eigvec(lanczos_step+1,lanczos_step+1))
allocate(lanczos_diag(lanczos_step+1))
allocate(lanczos_offdiag(lanczos_step+1))
eigvec=0.0d0
do i1=1,lanczos_step+1
eigvec(i1,i1)=1.0d0
end do
lanczos_diag=mat_diag
lanczos_offdiag=mat_coupling
call tql2(lanczos_step+1,lanczos_diag,lanczos_offdiag,eigvec,error_diag)
deallocate(lanczos_offdiag)
error_eigval=maxval(abs((eigval_lanczos-lanczos_diag(1:n_eigval))/eigval_lanczos))
eigval_lanczos=lanczos_diag(1:n_eigval)
eigvec_lanczos=eigvec(:,1:n_eigval)
if(error_eigval.lt.acc) error=0
print '(a,i3,a,f22.15)','Iteration ',lanczos_step,', Lanczos lowest eigenvalue: ',eigval_lanczos(1)
deallocate(lanczos_diag)
deallocate(eigvec)
end subroutine diagonalize_mat_and_check_finished
end subroutine lanczos
subroutine norm_vector_complex(vector,n,norm)
complex(kind=kind(0.0d0)), intent(inout) :: vector(:)
integer, intent(in) :: n
double precision, optional, intent(out) :: norm
! norm the vector, complex version
! vector(:): vector to norm
! n: size of the vector
! norm: calculated norm
integer :: i1
double precision :: norm_intern
norm_intern=0.0d0
do i1=1,n
norm_intern=norm_intern+abs(vector(i1))**2
end do
norm_intern=sqrt(norm_intern)
vector=vector/norm_intern
if(present(norm)) norm=norm_intern
end subroutine norm_vector_complex
subroutine lanczos_complex(sparse_mat,index_mat,n_mat,n_tot,n_eigval,eigval,acc,error&
&,eigvec,skip_diag,max_steps,ghost_state_level,startvec)
complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:)
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
integer, intent(inout) :: n_eigval
double precision, intent(in) :: acc
integer, intent(out) :: error
double precision, intent(out) :: eigval(:)
complex(kind=kind(0.0d0)), optional, intent(out) :: eigvec(:,:)
integer, optional, intent(in) :: skip_diag,max_steps
double precision, optional, intent(in) :: ghost_state_level
complex(kind=kind(0.0d0)), optional, intent(in) :: startvec(:)
! find maximum absolute value eigenvalues of a huge Hermitian symmetric matrix by the method of Lanczos, complex version
! sparse_mat(:): matrix to diagonalize in spares form, hermiticity considered
! index_mat(:,1): row indices of the spares matrix
! index_mat(:,2): column indices of the spares matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! n_eigval: number of eigenvalues that should be converged and calculated as input
! and number of calculated eigenvalues as output
! eigval(:): calculated eigenvalues
! acc: desired accuracy
! error: -1: convergence not reached in max_steps Lanczos steps
! 1: possibly a ghost state appeared
! 2: possibly a ghost state appeared in the first step
! 0: convergence reached
! eigvec(:,:): calculated eigenvectors, row: basis dimension, column: different energies
! skip_diag: skipping Lanczos steps between two convergence checks, > 0
! max_steps: maximum number of Lanczos steps
! ghost_state_level: minimum level for the coupling constant in the Lanczos matrix before subspace is diagonalized
! startvec(:): start vector for Lanczos method
integer :: i1,lanczos_step,skip_diag_intern,max_steps_intern,last_diag,seed
complex(kind=kind(0.0d0)), allocatable :: x_n_1(:),x_n_2(:)
double precision, allocatable :: mat_diag(:),mat_coupling(:),eigval_lanczos(:),eigvec_lanczos(:,:)
double precision :: ghost_state_level_intern
if(present(skip_diag)) then
skip_diag_intern=skip_diag
else
skip_diag_intern=4
end if
if(present(max_steps)) then
max_steps_intern=max_steps
else
max_steps_intern=20
end if
if(present(ghost_state_level)) then
ghost_state_level_intern=ghost_state_level
else
ghost_state_level_intern=1.0d-9
end if
error=-1
print '(/,a)','calculating Lanczos eigenvalues'
allocate(x_n_1(n_mat))
if(present(startvec)) then
x_n_1=startvec
else
seed=34123451
call randinit(seed)
call get_lanczos_vec_x_0(x_n_1,n_mat)
end if
allocate(x_n_2(n_mat))
allocate(mat_diag(max_steps_intern+1))
allocate(mat_coupling(max_steps_intern+1))
mat_coupling(1)=0.0d0
call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
if(mat_coupling(2).lt.ghost_state_level_intern) error=2
lanczos_step=1
last_diag=lanczos_step/skip_diag_intern
allocate(eigval_lanczos(n_eigval))
allocate(eigvec_lanczos(max_steps_intern+1,n_eigval))
eigval_lanczos=0.0d0
do while((lanczos_step.lt.max_steps_intern).and.(error.lt.0))
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
error=1
exit
end if
if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
if(error.ge.0) exit
last_diag=lanczos_step/skip_diag_intern
end if
if(lanczos_step.ge.max_steps_intern) exit
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
if(mat_coupling(lanczos_step+1).lt.ghost_state_level_intern) then
error=1
exit
end if
if((last_diag.ne.lanczos_step/skip_diag_intern).and.(lanczos_step.ge.n_eigval)) then
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
if(error.ge.0) exit
last_diag=lanczos_step/skip_diag_intern
end if
end do
if(error.eq.0) then
print '(/,a)','lanczos converged'
else if(error.eq.1) then
print '(/,a)','WARNING in lanczos: possible ghost state appeared'
call diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step-1,n_eigval&
&,eigval_lanczos,eigvec_lanczos,error,acc)
error=1
if(lanczos_step.lt.n_eigval) then
eigval_lanczos(lanczos_step+1:n_eigval)=getnan()
if(present(eigvec)) then
eigvec(:,lanczos_step+1:n_eigval)=getnan()
end if
n_eigval=lanczos_step
end if
lanczos_step=lanczos_step-1
else if(error.eq.2) then
print '(/,a)','WARNING in lanczos: possible ghost state appeared in the first step'
eigval_lanczos(1)=mat_diag(1)
if(n_eigval.gt.1) then
eigval_lanczos(2:n_eigval)=getnan()
if(present(eigvec)) then
eigvec(:,2:n_eigval)=getnan()
end if
n_eigval=1
end if
if(present(eigvec)) then
eigvec(:,1)=x_n_1
end if
else
print '(/,a,i0,a)','ERROR in lanczos: no convergence after ',max_steps_intern,' iterations'
end if
eigval=eigval_lanczos
deallocate(eigval_lanczos)
if(present(eigvec).and.((error.eq.0).or.(error.eq.1))) then
print '(/,a)','calculating Lanczos eigenvectors'
if(present(startvec)) then
x_n_1=startvec
else
call randinit_const(seed)
call get_lanczos_vec_x_0(x_n_1,n_mat)
end if
do i1=1,n_eigval
eigvec(:,i1)=x_n_1*eigvec_lanczos(1,i1)
end do
call get_lanczos_vec_x_1(x_n_1,sparse_mat,index_mat,n_mat,n_tot,x_n_2,mat_diag(1),mat_diag(2),mat_coupling(2))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(2,i1)
end do
max_steps_intern=lanczos_step
lanczos_step=1
do while(lanczos_step.lt.max_steps_intern)
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_2,x_n_1,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_1*eigvec_lanczos(lanczos_step+1,i1)
end do
if(lanczos_step.ge.max_steps_intern) exit
lanczos_step=lanczos_step+1
call get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n_1,x_n_2,mat_diag(lanczos_step)&
&,mat_coupling(lanczos_step),mat_diag(lanczos_step+1),mat_coupling(lanczos_step+1))
do i1=1,n_eigval
eigvec(:,i1)=eigvec(:,i1)+x_n_2*eigvec_lanczos(lanczos_step+1,i1)
end do
end do
end if
deallocate(eigvec_lanczos)
deallocate(x_n_1)
deallocate(x_n_2)
deallocate(mat_diag)
deallocate(mat_coupling)
contains
subroutine get_lanczos_vec_x_0(x_0,n)
complex(kind=kind(0.0d0)), intent(out) :: x_0(:)
integer, intent(in) :: n
! initialize starting vector for Lanczos method with random numbers
! x_0(:): starting vector
! n: size of the vector
integer :: i1
do i1=1,n
x_0(i1)=cmplx(rand(),rand())
end do
call norm_vector_complex(x_0,n)
end subroutine get_lanczos_vec_x_0
subroutine get_lanczos_vec_x_1(x_0,sparse_mat,index_mat,n_mat,n_tot,x_1,eigval_0,eigval_1,coupling_1)
complex(kind=kind(0.0d0)), intent(in) :: x_0(:),sparse_mat(:)
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
complex(kind=kind(0.0d0)), intent(out) :: x_1(:)
double precision, intent(out) :: eigval_0,eigval_1,coupling_1
! perform the first Lanczos step
! x_0(:): starting vector
! sparse_mat(:): matrix of the eigenvalue problem
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! x_1(:): calculated first Lanczos vector
! eigval_0: first diagonal term in the Lanczos matrix
! eigval_0: second diagonal term in the Lanczos matrix
! coupling_1: first offdiagonal term in the Lanczos matrix
integer :: i1
call get_expval(sparse_mat,index_mat,n_tot,x_0,eigval_0)
x_1=-eigval_0*x_0
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,1))
else
x_1(index_mat(i1,1))=x_1(index_mat(i1,1))+sparse_mat(i1)*x_0(index_mat(i1,2))
x_1(index_mat(i1,2))=x_1(index_mat(i1,2))+conjg(sparse_mat(i1))*x_0(index_mat(i1,1))
end if
end do
call norm_vector_complex(x_1,n_mat,coupling_1)
call get_expval(sparse_mat,index_mat,n_tot,x_1,eigval_1)
end subroutine get_lanczos_vec_x_1
subroutine get_lanczos_vec_x_n_plus_1(sparse_mat,index_mat,n_mat,n_tot,x_n,x_n1,eigval_n,coupling_n&
&,eigval_n1,coupling_n1)
complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:)
double precision, intent(in) :: eigval_n,coupling_n
integer, intent(in) :: index_mat(:,:),n_mat,n_tot
complex(kind=kind(0.0d0)), intent(inout) :: x_n(:),x_n1(:)
double precision, intent(out) :: eigval_n1,coupling_n1
! perform a general Lanczos step
! sparse_mat(:): matrix of the eigenvalue problem
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_mat: number of rows/columns of the matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! x_n(:): vector for the n-th interation
! x_n1(:): vector for the (n-1)-th iteration as input and (n+1)-th interation as output
! eigval_n: n-th diagonal term in the Lanczos matrix
! coupling_n: n-th offdiagonal term in the Lanczos matrix
! eigval_n1: (n+1)-th diagonal term in the Lanczos matrix
! coupling_n1: (n+1)-th offdiagonal term in the Lanczos matrix
integer :: i1
x_n1=-eigval_n*x_n-coupling_n*x_n1
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,1))
else
x_n1(index_mat(i1,1))=x_n1(index_mat(i1,1))+sparse_mat(i1)*x_n(index_mat(i1,2))
x_n1(index_mat(i1,2))=x_n1(index_mat(i1,2))+conjg(sparse_mat(i1))*x_n(index_mat(i1,1))
end if
end do
call norm_vector_complex(x_n1,n_mat,coupling_n1)
call get_expval(sparse_mat,index_mat,n_tot,x_n1,eigval_n1)
end subroutine get_lanczos_vec_x_n_plus_1
subroutine get_expval(sparse_mat,index_mat,n_tot,vector,expval)
complex(kind=kind(0.0d0)), intent(in) :: sparse_mat(:),vector(:)
integer, intent(in) :: index_mat(:,:),n_tot
double precision, intent(out) :: expval
! calculate the expectation value of a vector and the sparse matrix (hermiticity used)
! sparse_mat(:): matrix to calculate the expectation value
! index_mat(:,1): row indices of the sparse matrix
! index_mat(:,2): column indices of the sparse matrix
! n_tot: total number of non vanishing elements in the sparse matrix
! vector(:): vector to calculate the expectation value
! expval: expectation value
integer :: i1
expval=0.0d0
do i1=1,n_tot
if(index_mat(i1,1).eq.index_mat(i1,2)) then
expval=expval+real(sparse_mat(i1))*abs(vector(index_mat(i1,1)))**2
else
expval=expval+2.0d0*real(conjg(vector(index_mat(i1,1)))*sparse_mat(i1)*vector(index_mat(i1,2)))
end if
end do
end subroutine get_expval
subroutine diagonalize_mat_and_check_finished(mat_diag,mat_coupling,lanczos_step&
&,n_eigval,eigval_lanczos,eigvec_lanczos,error,acc)
double precision, intent(in) :: mat_diag(:),mat_coupling(:),acc
integer, intent(in) :: lanczos_step,n_eigval
double precision, intent(inout) :: eigval_lanczos(:)
double precision, intent(out) :: eigvec_lanczos(:,:)
integer, intent(out) :: error
! calculate the lowest eigenvalues in the spaned subspace and check if they converged
! mat_diag(:): main diagonal of the Lanczos matrix
! mat_coupling(:): first minor diagonal of the Lanczos matrix
! lanczos_step: current number of executed Lanczos steps
! n_eigval: number of eigenvalues that should be calculated
! eigval_lanczos(:): calculated lowest eigenvalues of the Lanczos matrix
! eigvec_lanczos(:,:): eigenvectors of the calculated lowest eigenvalues of the Lanczos matrix
! error: unchanged: no convergence for this check
! 0: convergence reached
! acc: desired accuracy
double precision, allocatable :: lanczos_diag(:),lanczos_offdiag(:),eigvec(:,:)
double precision :: error_eigval
integer :: i1,error_diag
allocate(eigvec(lanczos_step+1,lanczos_step+1))
allocate(lanczos_diag(lanczos_step+1))
allocate(lanczos_offdiag(lanczos_step+1))
eigvec=0.0d0
do i1=1,lanczos_step+1
eigvec(i1,i1)=1.0d0
end do
lanczos_diag=mat_diag
lanczos_offdiag=mat_coupling
call tql2(lanczos_step+1,lanczos_diag,lanczos_offdiag,eigvec,error_diag)
deallocate(lanczos_offdiag)
error_eigval=maxval(abs((eigval_lanczos-lanczos_diag(1:n_eigval))/eigval_lanczos))
eigval_lanczos=lanczos_diag(1:n_eigval)
eigvec_lanczos=eigvec(:,1:n_eigval)
if(error_eigval.lt.acc) error=0
print '(a,i3,a,f22.15)','Iteration ',lanczos_step,', Lanczos lowest eigenvalue: ',eigval_lanczos(1)
deallocate(lanczos_diag)
deallocate(eigvec)
end subroutine diagonalize_mat_and_check_finished
end subroutine lanczos_complex
subroutine gram_schmidt_onb(vectors,n_vec,n_basis)
double precision, intent(inout) :: vectors(:,:)
integer, intent(in) :: n_vec,n_basis
! build an orthonormal basis (ONB) out of a set of linear independent vectors
! vectors(:,:): linear independent vectors as input and ONB as output
! one vector is stored as a column vector v_n=vectors(:,n)
! n_vec: number of vectors that are handed over
! n_basis: number of components of one vector
double precision :: overlap_v
integer :: i1,i2,i3
call norm_vector(vectors(:,1),n_basis)
do i1=2,n_vec
do i2=1,i1-1
overlap_v=0.0d0
do i3=1,n_basis
overlap_v=overlap_v+vectors(i3,i1)*vectors(i3,i2)
end do
vectors(:,i1)=vectors(:,i1)-overlap_v*vectors(:,i2)
end do
call norm_vector(vectors(:,i1),n_basis)
end do
end subroutine gram_schmidt_onb
subroutine gram_schmidt_on_single_vector(onb,vector,n_vec,n_basis)
double precision, intent(inout) :: onb(:,:),vector(:)
integer, intent(in) :: n_vec,n_basis
! add one vector to an orthonormal basis (ONB)
! onb(:,:): orthonormal basis with the vectors in the columns
! vector(:): vector to add to the ONB
! n_vec: number of vectors that are handed over in the ONB
! n_basis: number of components of one vector
double precision :: overlap_v
integer :: i1,i2
do i1=1,n_vec
overlap_v=0.0d0
do i2=1,n_basis
overlap_v=overlap_v+onb(i2,i1)*vector(i2)
end do
vector=vector-overlap_v*onb(:,i1)
end do
call norm_vector(vector,n_basis)
end subroutine gram_schmidt_on_single_vector
end module numeric_lib
|