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 |