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

[idomeneas]

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

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

Τα

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

W,WORK

είναι REAL(8). Δες και εδώ πριν χρησιμοποιείς κάτι από LAPACK http://www.netlib.org/lapack/double/dsyev.f Μετά από αυτό θα έχεις ένα LWORK όπως το είπαμε, και το allocate θα γίνει για το WORK(LWORK). Τέλος σαν παρατήρηση. Για να χρησιμοποιείς DSYEV, θα πρέπει αυτό που λύνεις να είναι όντως ένα συμμετρικό πρόβλημα και τότε δεν υπάρχουν Complex μεταβλητές γιατί τα συμμετρικά μητρώα ως γνωστόν έχουν πραγματικές ιδιοτιμές-ιδιοδιανύσματα
Έτσι στο interface θα αλλάξεις το COMPLEX(8) σε REAL(8), και το ίδιο στο main program. Το WORK είναι το workspace που χρειάζεται η lapack και θα γίνει allocate με την τιμή του LWORK, και το αποτέλεσμα (eigenvalues) θα γίνει allocate για την τιμή step+1. Σε αυτό θα περιέχονται οι ιδιοτιμές (W στο interface).Κάνε αυτά και μετά πες αν δούλεψε.

[yallou]

Το πρόβλημα λύθηκε και επιτέλους έχω τις ιδιοτιμές του πίνακα .Το 1ο λάθος ήταν στο interface τα w,work έπρεπε όντως για το πρόβλημα αυτό να είναι real(8) και όχι complex(8),και επίσης έπρεπε να τα δηλώσω έτσι:

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


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


για να μην έχω Segmentation fault και τo 2o λάθος ήταν στο allocation που έκανα.
Dimitris και idomeneas ευχαριστώ πάρα πολύ για πολλοστή φορά

[yallou]

Μια ετεροχρονισμένη ερώτηση:Λύνοντας το ίδιο πρόβλημα στην Octave και καλώντας την εντολή j_eig=eig(jac) οι ιδιοτιμές που υπολογίζονται είναι διαφορετικές απο αυτές που έχω υπολογίσει με την gfortran/Lapack. Γιατί συμβαίνει αυτό.Αν έχω καταλάβει καλά η Octave χρησιμοποιεί τις ρουτινές της Lapack για αυτούς τους υπολογισμούς.Ενδεικτικά οι τιμές των ιδιοτιμών με την Octave είναι:

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


16373.817
16343.292
16292.502
16221.573
16130.681
16020.052
15889.962
15740.733
15572.737
15386.391
15182.159
14960.548
14722.109
14467.436
14197.161
13911.956
13612.53
13299.628
12974.028
12636.539
12288
11929.278
11561.264
11184.874
10801.043
10410.725
10014.891
9614.5259
9210.6236
8804.1889
8396.2322
7987.7678
7579.8111
7173.3764
6769.4741
6369.1085
5973.2749
5582.9574
5199.1263
4822.7361
4454.7223
10.183301
40.707888
91.497872
162.42698
253.31887
363.94757
494.03805
643.26687
811.26303
997.60887
1201.8411
1423.452
1661.8905
1916.5639
2186.8391
2472.044
2771.4696
3084.3715
3409.9718
3747.461
4096
-1.0
-1.0


και με την fortran:

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


16408.44313358766158
16374.36435669779348
16325.90236972711682
16257.95101013651401
16170.61500033698030
16064.18017674126349
15938.93573956556429
15795.19616984482855
15633.31140972829053
15453.67050850839041
15256.70265539303728
15042.87708480910260
14812.70243109937837
14566.72577411609927
14305.53148520734612
14029.73992740050562
13740.00603858639442
13437.01781467786532
13121.49470387844849
12794.18592022473058
12455.86868303007759
12107.34638806607654
11749.44671592963641
11383.01968287524869
11008.93563935048769
10628.08322150335880
10241.36726100321175
9849.70665861724046
9454.03222710077353
9055.28450908944251
8654.41157582096275
8252.36681266608139
7850.10669761243025
7448.58857902604723
7048.76845921635777
6651.59879055659985
6258.02629116909702
5868.98978747757974
5485.41809126300632
5108.22791923642671
4738.32186356064994
4376.58642220119236
4023.89009844151133
3681.08157931005462
3348.98800295331148
3028.41332501155193
2720.13679359562502
2424.91154119934117
2143.46329935270342
1876.48923742124362
1624.65691998182024
1388.60336704601173
1168.93418806815271
966.22274606003941
781.00929895507397
613.80008014964028
465.06636705391105
335.24387317190332
224.73364008856848
133.90790149036312
63.12622798152661
12.58800685344813
-1.00000000000000
-22.31511396855747
-2049.10568108890038


[logari81]

καταρχήν για να τις συγκρινεις κανετες ενα sorting πρώτα.

[yallou]

η 1η στήλη είναι οι τιμές απο την Octave και η 2η στήλη είναι οι τιμές απο την Fortran

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


16375.529 16408.443
16345.97 16374.364
16296.752 16325.902
16228.009 16257.951
16139.909 16170.615
16032.663 16064.18
15906.53 15938.936
15761.815 15795.196
15598.866 15633.311
15418.075 15453.671
15219.879 15256.703
15004.754 15042.877
14773.219 14812.702
14525.831 14566.726
14263.188 14305.531
13985.92 14029.74
13694.697 13740.006
13390.22 13437.018
13073.223 13121.495
12744.468 12794.186
12404.749 12455.869
12054.883 12107.346
11695.714 11749.447
11328.106 11383.02
10952.945 11008.936
10571.136 10628.083
10183.597 10241.367
9791.2632 9849.7067
9395.0787 9454.0322
8995.9982 9055.2845
8594.9833 8654.4116
8193 8252.3668
7791.0167 7850.1067
7390.0018 7448.5886
6990.9213 7048.7685
6594.7368 6651.5988
6202.4027 6258.0263
5814.864 5868.9898
5433.0545 5485.4181
5057.894 5108.2279
4690.2862 4738.3219
4331.1167 4376.5864
3981.2509 4023.8901
3641.5315 3681.0816
3312.777 3348.988
2995.7794 3028.4133
2691.3024 2720.1368
2400.0794 2424.9115
2122.8121 2143.4633
1860.1683 1876.4892
1612.781 1624.6569
1381.2459 1388.6034
1166.1209 1168.9342
967.92434 966.22275
787.13355 781.0093
624.18411 613.80008
479.46856 465.06637
353.3355 335.24387
246.0887 224.73364
157.98637 133.9079
89.24023 63.126228
40.013854 12.588007
10.411118 -1
-1 -22.315114
-1 -2049.1057


[Dimitris]

Ποιος είναι ο πινακας σου για να το τσεκάρουμε κι εμείς;

[yallou]

Ο πίνακας είναι o "jac". Τον θέτω ίσο με τον jacobian-δες αν θέλεις τον κώδικα που έχω ανεβάσει στην προηγούμενη σελίδα-πριν ξεκινήσει η απαλοιφή Gauss.

[yallou]

Για να μην σας μπερδεύω ο τελικός κώδικας στην 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



[idomeneas]

Χρησιμοποίησε την dgeev συνάρτηση της Lapack. Αν είναι σωστά αυτά που έχεις κάνει με την συμμετρική, dsyev, θα πρέπει να βγάλεις τα ίδια αποτελέσματα μεταξύ αυτών. Τέλος, αν δε με απατά η μνήμη μου το Octave χρησιμοποιεί τη dgeev διότι δεν ξέρει αν εσύ έχεις συμμετρικό ή μη, μητρώο. Τέλος σύγκρινε όλες αυτές τις περιπτώσεις μεταξύ τους και πες τα συμπεράσματα σου. Για εγχειρίδιο της dgeev γράψε

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

man dgeev

αν έχεις εγκαταστήσει το manual της lapack ή ψάξε στο δίκτυο στο netlib.org

[yallou]

Idomeneas για μια ακόμη φορά είχες δίκιο . Χρησιμοποιώντας την dgeev έχω τις ίδιες ιδιοτιμές με την Octave.
Στην 1η στήλη οι ιδιοτιμες απο την Octave και στην 2η απο την Fortran

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


16375.529 16375.529
16345.97 16345.97
16296.752 16296.752
16228.009 16228.009
16139.909 16139.909
16032.663 16032.663
15906.53 15906.53
15761.815 15761.815
15598.866 15598.866
15418.075 15418.075
15219.879 15219.879
15004.754 15004.754
14773.219 14773.219
14525.831 14525.831
14263.188 14263.188
13985.92 13985.92
13694.697 13694.697
13390.22 13390.22
13073.223 13073.223
12744.468 12744.468
12404.749 12404.749
12054.883 12054.883
11695.714 11695.714
11328.106 11328.106
10952.945 10952.945
10571.136 10571.136
10183.597 10183.597
9791.2632 9791.2632
9395.0787 9395.0787
8995.9982 8995.9982
8594.9833 8594.9833
8193.0000 8193.0000
7791.0167 7791.0167
7390.0018 7390.0018
6990.9213 6990.9213
6594.7368 6594.7368
6202.4027 6202.4027
5814.864 5814.864
5433.0545 5433.0545
5057.894 5057.894
4690.2862 4690.2862
4331.1167 4331.1167
3981.2509 3981.2509
3641.5315 3641.5315
3312.777 3312.777
2995.7794 2995.7794
2691.3024 2691.3024
2400.0794 2400.0794
2122.8121 2122.8121
1860.1683 1860.1683
1612.781 1612.781
1381.2459 1381.2459
1166.1209 1166.1209
967.92434 967.92434
787.13355 787.13355
624.18411 624.18411
479.46856 479.46856
353.3355 353.3355
246.0887 246.0887
157.98637 157.98637
89.24023 89.24023
40.013854 40.013854
10.411118 10.411118
-1.00000 -1.00000
-1.00000 -1.00000


Ο πίνακας Jac είναι συμμετρικός και τριδιαγώνιος δεν θα έπρεπε οι 2 υπορουτίνες να δίνουν τις ίδιες ιδιοτιμές

Off topic:
Η αντίστοιχη υπορουτίνα με την devcrg της IMSL που καλούσα στο εργαστήριο της σχολής τελικά είναι η dgeev και όχι η dseyv...