Προγραμματισμός σε fortran

[gman]

Off topic:
Δεν έχετε και άδικο ρε παιδιά. Μας μαθαίνου Fortran και C, έχουν απαιτήσεις, αλλά ποτέ κανένας δεν μας εξήγησε τι είναι ένα IDE ή ένα πρότυπο γλώσσας προγραμματισμού. Προσπαθώ κι εγώ μέσω ίντερνετ να μάθω κάποια πράγματα, έχω κατεβάσει και μερικά e-books και κάτι γίνεται. Δεν είναι όμως πάντα εύκολο. Ευτυχώς υπάρχουν άτομα που βοηθούν.

[Dimitris]

Μη νομίζεις κι εμείς έτσι μαθαμε, από μόνοι μας. Γι'αυτό έγραψα κι εδώ viewtopic.php?f=9&t=10534 μερικά βασικά.

[sotos_ubu]

Και εγώ στο πανεπιστήμιο διδάχθηκα στο πρώτο εξάμηνο Fortran και στο δεύτερο C. Θυμάμαι ότι ο καθηγητής στο μάθημα αλλά και οι διδακτορικοί φοιτητές στο εργαστήριο (στα πρώτα μαθήματα) μας εξήγησαν τι είναι και τι κάνουν ακριβώς αυτά (compiler, linker, IDE, πρότυπα της Fortran). Απλά το αναφέρω επειδή διαβάζοντας αυτό το θέμα, μου δημιουργήθηκε μια αίσθηση μηδενισμού κάποιων καταστάσεων..

[logari81]

μου δημιουργήθηκε μια αίσθηση μηδενισμού κάποιων καταστάσεων

οι απόψεις που αναφέρονται ειναι προσωπικές και αφορούν τις εμπειρίες του καθενός. Προσωπικά όχι μονο για προγραμματισμο αλλα γενικότερα έχω τυχει εξαίρετους καθηγητες αλλα και ανεκδιήγητους καθηγητες. Δεν μπορουμε να μηδενιζουμε τα πάντα αλλα απο εκει και περα δεν μπορουμε να κλεινουμε τα ματια σε μια πραγματικότητα που θελει πολλους καθηγητές στην θεση τους λόγω πολιτικών ή αλλων διασυνδεσεων κι όχι ικανοτήτων.

[Dimitris]

Και εγώ στο πανεπιστήμιο διδάχθηκα στο πρώτο εξάμηνο Fortran και στο δεύτερο C. Θυμάμαι ότι ο καθηγητής στο μάθημα αλλά και οι διδακτορικοί φοιτητές στο εργαστήριο (στα πρώτα μαθήματα) μας εξήγησαν τι είναι και τι κάνουν ακριβώς αυτά (compiler, linker, IDE, πρότυπα της Fortran). Απλά το αναφέρω επειδή διαβάζοντας αυτό το θέμα, μου δημιουργήθηκε μια αίσθηση μηδενισμού κάποιων καταστάσεων..

Off topic:
Φυσικά προσωπικές εμπειρίες περιγράφω. Απλά είμαι λιγάκι απογοητευμένος με την ελληνική εκπαίδευση. Τουλάχιστον προσπαθώ να τη διορθώση με τον τρόπο μου (δες στο link που έδωσα παραπάνω) Και λίγη εποικοδομητική κριτική δε βλάπτει. Και μιας και λέμε για εποικοδομητική ριξτε μια ματιά σ'αυτό όσοι ενδιαφέρονται για fortran:

Στη fortran 90/95 υποστηρίζονται δομές ή νέοι τύποι δεδομένων και υπερφόρτωση τελεστών. Ας δούμε μια απλή εφαρμογή τους. Από quaternions πώς τα πάμε; Quaternions είναι η επέκταση των μιγαδικών στις 4 διαστάσεις (Σημείωση ότι δεν ορίζονται τριαδικοί αριθμοί όσο κι αν προσπάθησε ο καημένος ο Hamilton) Επίσης οι τετραδικοί αριθμοί δεν είναι αντιμεταθετικοί! Για περισσότερες λεπτομέρειες σχετικά με το πώς ορίζονται οι πράξεις μπορείτε να βρείτε εδώ http://en.wikipedia.org/wiki/Quaternion Οι τετραδικοί αριθμοί όσο παράξενοι κι αν φαίνεται κάνουν τις περιστροφές στις 3 διαστάσεις πολύ πιο γρήγορα απ'ότι οι γωνίες Euler (logari81, αυτό ίσως είναι για σένα. Μια συνάδελφος τους εφάρμοσε στις περιστροφές του πλέγματός της.)

Επί της ουσίας. Ένα μικρό παράδειγμα που μας επιτρέπει να γράψουμε q = w + v, όπου όλοι είναι τετραδικοί. Το module που κάνει τη δουλειά:

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

module structure
implicit none

intrinsic :: conjg, abs

private
public :: quaternion_print
public :: operator(+), operator(-), operator(*)
public :: conjg, abs

!###########################################################################################################

type, public :: quaternion
real(8) :: t
real(8) :: x
real(8) :: y
real(8) :: z
end type quaternion

!###########################################################################################################

interface operator(+)
module procedure quat_add, quat_integer_add, quat_real_add, quat_real8_add
end interface

interface operator(-)
module procedure quat_subt
end interface

interface operator(*)
module procedure quat_integer_left_mult, quat_real_left_mult, quat_real8_left_mult, &
quat_integer_right_mult, quat_real_right_mult, quat_real8_right_mult, &
quat_mult
end interface

interface conjg
module procedure quat_conjg
end interface

interface abs
module procedure quat_abs
end interface

!###########################################################################################################

contains

!###########################################################################################################

subroutine quaternion_print(q)
type(quaternion), intent(in) :: q
real :: tp, xp, yp, zp
tp = q % t
xp = q % x
yp = q % y
zp = q % z

print "(a1, 4f32.12, a1)", "(", tp, xp, yp, zp, ")"
print *
end subroutine quaternion_print

!###########################################################################################################

type(quaternion) function quat_add(a,b) result (res)
type(quaternion), intent(in) :: a, b
!type(quaternion) :: res
res % t = a % t + b % t
res % x = a % x + b % x
res % y = a % y + b % y
res % z = a % z + b % z
end function quat_add

type(quaternion) function quat_integer_add(i,q) result (res)
integer, intent(in) :: i
type(quaternion), intent(in) :: q
!type(quaternion) :: res
res % t = q % t + i
res % x = q % x
res % y = q % y
res % z = q % z
end function quat_integer_add

type(quaternion) function quat_real_add(r,q) result (res)
real, intent(in) :: r
type(quaternion), intent(in) :: q
!type(quaternion) :: res
res % t = q % t + r
res % x = q % x
res % y = q % y
res % z = q % z
end function quat_real_add

type(quaternion) function quat_real8_add(r,q) result (res)
real(8), intent(in) :: r
type(quaternion), intent(in) :: q
!type(quaternion) :: res
res % t = q % t + r
res % x = q % x
res % y = q % y
res % z = q % z
end function quat_real8_add

!###########################################################################################################

type(quaternion) function quat_subt(a,b) result (res)
type(quaternion), intent(in) :: a, b
!type(quaternion) :: res
res % t = a % t - b % t
res % x = a % x - b % x
res % y = a % y - b % y
res % z = a % z - b % z
end function quat_subt

!###########################################################################################################

type(quaternion) function quat_integer_left_mult(r,a) result (res)
integer, intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_integer_left_mult

type(quaternion) function quat_real_left_mult(r,a) result (res)
real, intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_real_left_mult

type(quaternion) function quat_real8_left_mult(r,a) result (res)
real(8), intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_real8_left_mult

type(quaternion) function quat_integer_right_mult(a,r) result (res)
integer, intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_integer_right_mult

type(quaternion) function quat_real_right_mult(a,r) result (res)
real, intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_real_right_mult

type(quaternion) function quat_real8_right_mult(a,r) result (res)
real(8), intent(in) :: r
type(quaternion), intent(in) :: a
!type(quaternion) :: res
res % t = r * a % t
res % x = r * a % x
res % y = r * a % y
res % z = r * a % z
end function quat_real8_right_mult

type(quaternion) function quat_mult(u,v) result (res)
type(quaternion), intent(in) :: u, v
!type(quaternion) :: res
res % t = u % t * v % t - u % x * v % x - u % y * v % y - u % z * v % z
res % x = u % t * v % x + u % x * v % t + u % y * v % z - u % z * v % y
res % y = u % t * v % y + u % y * v % t + u % z * v % x - u % x * v % z
res % z = u % t * v % z + u % z * v % t + u % x * v % y - u % y * v % x
end function quat_mult

!###########################################################################################################

type(quaternion) function quat_conjg(u) result (res)
type(quaternion), intent(in) :: u
!type(quaternion) :: res
res % t = u % t
res % x = - (u % x)
res % y = - (u % y)
res % z = - (u % z)
end function quat_conjg

!###########################################################################################################

function quat_abs(u) result (res)
type(quaternion), intent(in) :: u
real(8) :: res
res = sqrt((u % t) **2 + (u % x) **2 + (u % y) **2 + (u % z) **2)
end function quat_abs

!###########################################################################################################

end module





και το πρόγραμμα που το χρησιμοποιεί:

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

program main
use structure
type (quaternion) :: u, v, w
integer :: i
integer :: a
real :: b
real(8) :: c
complex :: z

! open(10,file="outpout")


a = 5
b = 0.35
c = -176.0987
u = quaternion(1,2,3,4)
v = quaternion(5,6,7,8)

w=u+v
call quaternion_print(w)
w=u-v
call quaternion_print(w)
w = a * u
call quaternion_print(w)
w = b * u
call quaternion_print(w)
w = c * u
call quaternion_print(w)
w = u * a
call quaternion_print(w)
w = u * b
call quaternion_print(w)
w = u * c
call quaternion_print(w)
w = u * v
call quaternion_print(w)
w = v * u
call quaternion_print(w)

z = (1,2)
print *,z
print *, abs(z)
print *, abs(w)
print*, conjg(w)

w = c + w
call quaternion_print(w)

end program


[yallou]

Δεν ξέρω αν έπρεπε να ανοίξω νέο θεμα αλλά νομίζω οτί η ερώτηση που έχω "κολλάει" στην συγκέκριμενη ενότητα.Προσπαθώ να υπολογίσω τις ιδιοτιμές του Ιακωβιανού πίνακα όπως αυτός προκύπτει σε ένα πρόβλημα συνοριακών τιμών καλώντας την υπορουτινά 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


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

[Dimitris]

Σου το λέει το σφάλμα. Η 8η είναι λάθος, δηλαδή η LWORK. Για τις διαστάσεις διάβασε τα σχόλια της DSYEV.

[idomeneas]

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

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

Το LWORK πρέπει να είναι >=max(1,3*N-1). Οχι step+1 αλλά 3*step+2

[yallou]

Ευχαριστώ για το άμεσο ενδοαφέρον (Dimitris) και την προτεινόμενη λύση (idomeneas) αλλά δυστύχως μέχρι στιγμής το πρόβλημα παραμένει.Αλλάζοντας το Lwork σε 3*step+2 έχω Segmentation fault και βάζοντας την τιμή 3*step-1 παίρνω πάλι το μήνυμα ** On entry to DSYEV parameter number 8 had an illegal value.Θα επιχειρήσω να ξαναγράψω το interface γιατί δεν είμαι πολύ σίγουρος για το κατά πόσο είναι σώστο.

[Dimitris]

Και το allocate πρέπει να αλλάξεις όχι μόνο το όρισμα της συνάρτησης.