Δημοσιεύτηκε: 27 Μάιος 2010, 21:26
Για να μην σας μπερδεύω ο τελικός κώδικας στην fortran είναι ο παρακάτω:
- Κώδικας: Επιλογή όλων
module all_external_subroutines
interface eigen_interface
subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info)
character(1)::jobz,uplo
integer::n,lda,lwork,info
real(8),dimension(1:n)::w
real(8),dimension(1:lwork)::work
real(8),dimension(1:n,1:n)::a
end subroutine dsyev
end interface
end module all_external_subroutines
!=======================================================================
program p_bvp
use all_external_subroutines
implicit none
real(8)::t,dt,a=0.d0,b=1.d0,tol1,tol2,pi,ts,tf
real(8),allocatable,dimension(:)::f,rhs,w1,work1
real(8),allocatable,dimension(:,:)::jacobian,jac
integer::k,i,j,iter,gauss_error,step,lda1,lwork1,info
pi=4.d0*datan(1.d0)
k=6
tol2=0.5e-6
step=2**k; print*, "Total points = " ,step
dt=(b-a)/step
lda1=step+1
lwork1=3*step+2
call cpu_time(ts)
iter=0
allocate(f(0:step),jacobian(0:step,0:step),rhs(0:step),jac(0:step,0:step))
f=0.d0
iter=iter+1
do while (iter<=1000)
do i=1,step-1
t=a+i*dt
rhs(i)=((f(i+1)-2.d0*f(i)+f(i-1))/dt**2)+f(i)*((f(i+1)-f(i-1))/(2.d0*dt))
end do
rhs(0)= f(0)-1.d0
jacobian(0,:) = 0.d0
jacobian(0,0) = 1.d0
do i=1,step-1
do j=0,step
if (j==i-1)then
jacobian(i,i-1)=1.d0/dt**2-f(i)/(2.d0*dt)
else if (j==i+1) then
jacobian(i,i+1)=1.d0/dt**2+f(i)/(2.d0*dt)
else if (i==j) then
jacobian(i,i)=-(2.d0/dt**2)+((f(i+1)-f(i-1))/(2.d0*dt))
else
jacobian(i,j)=0.d0
end if
end do
end do
rhs(step)= f(step)-0.d0
jacobian(step,:)=0.d0
jacobian(step,step)=1.d0
jac=jacobian
call gauss_elimination( matrix=jacobian, rhs=rhs, n=step+1, error=gauss_error )
if( gauss_error /= 0 ) then
write(*,*), "Gauss elimantion error: ", gauss_error
stop "Program can not continue"
end if
f = f - rhs
tol1=maxval(dabs(rhs)); print*, "Max error = ", tol1
if(tol1<tol2) exit
iter=iter+1
end do
open (10,file="out_k=6.dat")
do i=0,step
t=a+i*dt
write (10,'(2f16.10)')t,f(i)
end do
close (10)
deallocate(f,rhs,jacobian)
allocate(w1(1:step+1),work1(1:lwork1))
call dsyev(jobz='v',uplo='u',n=step+1,a=jac,lda=lda1,w=w1,work=work1,lwork=lwork1,info=info)
open(14,file="com_eig_u_k=6.dat")
PRINT 1050
PRINT 1060, -w1
do i=1,size(w1)
write(14,'(F22.16)')-w1(i)
end do
print 1010,"w1 size =",size(w1)
print 1010,"work1 size =",size(work1)
close(14)
deallocate(w1,work1,jac)
call cpu_time(tf)
print 1100,"time passed =",tf-ts,"sec"
print 1020,"info =",info
1100 format(a14,f9.4,a4)
1020 format(a8,i3.1)
1010 format(a14,i5.1)
1050 FORMAT (/6x, 'Eigenvalues')
1060 FORMAT (1x, f22.16)
contains
pure subroutine gauss_elimination( matrix, rhs, n, error, log_det, det_sign )
implicit none
integer, intent(in) :: n
real(8), intent(inout), dimension (1:n,1:n) :: matrix
real(8), intent(inout), dimension (1:n) :: rhs
integer, intent(out) :: error
real(8), intent(out), optional :: log_det
integer, intent(out), optional :: det_sign
integer :: p, r, i, j, counter
real(8) :: maxv, m, eps
real(8), dimension (1:n) :: temp
error = 0; counter = 0; eps = 10.d0*epsilon(1.d0)
if( n == 1 ) then
if( matrix(1,1) /= 0.d0 ) then
rhs(1) = rhs(1) / matrix(1,1)
return
else
error = + 2
endif
endif
all_steps: do r = 1, n-1
maxv = dabs(matrix(r,r)); p = r
do i = r+1,n
if( dabs(matrix(i,r)) > maxv ) then
maxv = dabs(matrix(i,r))
p = i
endif
enddo
if(dabs(maxv)<eps)then; error = +1; return; endif
if( p /= r ) then
counter = counter + 1
temp(r:n) = matrix(r,r:n); m = rhs(r)
matrix(r,r:n) = matrix(p,r:n); rhs(r) = rhs(p)
matrix(p,r:n) = temp(r:n); rhs(p) = m
endif
rows: do i = r+1, n
m = matrix(i,r) / matrix(r,r)
columns: do j = r+1, n
matrix(i,j) = matrix(i,j) - m*matrix(r,j)
enddo columns
rhs(i) = rhs(i) - m*rhs(r)
enddo rows
enddo all_steps
! ... if required, find the log10 and the sign of the determinant ........
if(present(log_det).and.present(det_sign))then
log_det = 0.d0
det_sign = 1
do i = 1, n
log_det = log_det + dlog10(dabs(matrix(i,i)))
if( matrix(i,i) < 0 ) det_sign = - det_sign
enddo
if(mod(counter,2)==1) det_sign = - det_sign
endif
! ..... back substitution ....................
rhs(n) = rhs(n) / matrix(n,n)
do i = n-1, 1, -1
m = 0.d0
do j = n, i+1, -1
m = m + matrix(i,j)*rhs(j)
enddo
rhs(i) = ( rhs(i) - m ) / matrix(i,i)
enddo
end subroutine gauss_elimination
end program p_bvp