Ubuntu-gr Forum

Δεν ξέρω αν έπρεπε να ανοίξω νέο θεμα αλλά νομίζω οτί η ερώτηση που έχω "κολλάει" στην συγκέκριμενη ενότητα.Προσπαθώ να υπολογίσω τις ιδιοτιμές του Ιακωβιανού πίνακα όπως αυτός προκύπτει σε ένα πρόβλημα συνοριακών τιμών καλώντας την υπορουτινά DSEYV της lapack γράφοντας την παρακάτω εντολή: gfortran -lblas -llapack -o bvp bvp.f90. Τρέχοντας το εκτελέσιμο αρχείο μου βγάζει το εξής μήνυμα:** On entry to DSYEV parameter number 8 had an illegal value. Αντιλάμβανομαι ότι μάλλον έχω ορίσει λάθος κάποια/κάποιες μεταβλητές αλλά ειλικρινά δεν μπορώ να καταλάβω που...
Επιπλέον στο εργαστήριο της σχολής χρησιμοποιώντας τον compiler της Compaq (CVF 6.5) και την αντίστοιχη υπορουτινά που περιλαμβάνεται στην βιβλιοθήκη της IMSL (DEVCRG) δεν έχω κάποια πρόβλημα.Η τελευταία "εκδοχή" του κώδικα είναι η παρακάτω

Κώδικας: Επιλογή όλων


module all_external_subroutines
interface eigen_interface

SUBROUTINE DSYEV(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)

CHARACTER(LEN=1)::JOBZ, UPLO
INTEGER::N,LDA,LWORK,INFO
COMPLEX(8),DIMENSION(:)::W,WORK
REAL(8),DIMENSION(:,:) :: 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,tmp
real(8),dimension(:),allocatable::f,rhs
real(8),dimension(:,:),allocatable::jacobian, jac
character(len=1)1::jobz,uplo
complex(8),dimension(:,:),allocatable::eigen_vectors
complex(8),dimension(:), allocatable::eigen_values,w
integer::k,i,j,iter,gauss_error,step,err,ierror,info

k=6
tol2=0.5e-6
step=2**k; print*, "Total points = " ,step
dt=(b-a)/step

iter=0

allocate(f(0:step),jacobian(0:step,0:step),rhs(0:step),jac(0:step,0:step))
allocate(w(step))
f=0.d0

iter=iter+1

do while (iter<=100)
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.dat")
do i=0,step
t=a+i*dt
write (10,'(2f16.10)')t,f(i)
end do

close (10)

deallocate(jacobian)
allocate( eigen_values(1:step+1), eigen_vectors(1:step+1,1:step+1), stat = err ); if( err /= 0 ) stop

call dsyev(JOBZ='V' ,UPLO='L' ,N=step+1 ,A=jac ,LDA=step+1 ,W=w ,WORK=eigen_values ,LWORK=step+1 ,INFO=info)

print*, -eigen_values

print *,' info =', info

do i = 1, step+1
tmp = cdabs( dot_product( eigen_vectors(:,i), eigen_vectors(:,i) ) )
if( dabs( tmp - 1.d0 ) > 1.d-13 ) print*, "Error for i = ", i
end do

deallocate(f,rhs,jac)

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


Υ.Γ. :Αν είμαι σε λάθος ενότητα ζητάω συγγνώμη από τους διαχειριστές