- #1
dilasluis
- 32
- 0
I'm testing a part of a code and I have a problem. After the subroutine call at line 44 (CALL ZEIGSUB), nothing is written in UNIT NRES. The other units are written correctly. Here's the code:
CONTAINS
THANKS
Fortran:
PROGRAM TEST
IMPLICIT NONE
complex *16, allocatable :: KSTAR(:,:), VR(:,:), CM2(:,:)
real *8, allocatable :: M2(:,:), FREQ(:), MLF(:)
integer *4 :: DOFMECH2, TTT, I, MOI
REAL *8 :: OM, OMM1, TESTOM
REAL(8), PARAMETER :: PI = 4.*datan(1.d+0)
integer *4, parameter :: NRES = 9
open(unit=NRES,file="/Users/ndilokelwaluis/Documents/Fortran/Vibration/MODEL_V1-0/moiresul.dat",status='unknown')
OM = 60
MOI = 3
DOFMECH2 = 4
allocate(KSTAR(4,4))
allocate(VR(4,4))
allocate(M2(4,4))
ALLOCATE(CM2(4,4))
allocate(FREQ(4))
allocate(MLF(4))
KSTAR(1,1) = (-21.10,-22.50)
KSTAR(1,2) = ( 53.50,-50.50)
KSTAR(1,3) = (-34.50,127.50)
KSTAR(1,4) = ( 7.50, 0.50)
KSTAR(2,1) = ( -0.46, -7.78)
KSTAR(2,2) = ( -3.50,-37.50)
KSTAR(2,3) = (-15.50, 58.50)
KSTAR(2,4) = (-10.50, -1.50)
KSTAR(3,1) = ( 4.30, -5.50)
KSTAR(3,2) = ( 39.70,-17.10)
KSTAR(3,3) = (68.50, 12.50)
KSTAR(3,4) = ( -7.50, -3.50)
KSTAR(4,1) = ( 5.50, 4.40)
KSTAR(4,2) = ( 14.40, 43.30)
KSTAR(4,3) = (32.50,46.00)
KSTAR(4,4) = (19.00,32.50)
M2 = reshape( (/ 1.0,1.6,3.0,0.0,0.8,3.0,4.0,2.4,1.0,2.4,4.0,0.0,0.0,1.8,&
0.0,4.0 /),(/4,4/) )
CM2 = DCMPLX(M2)
write(NRES,*) 'BEFORE'
CALL ZEIGSUB(KSTAR,CM2,DOFMECH2,FREQ,MLF,VR)
write(NRES,*) 'AFTER'
DO 1010, I = 1, DOFMECH2
IF (FREQ(I) .GT. 10 .AND. MLF(I) .GT. 0) THEN
TTT = I
GO TO 1020
END IF
1010 CONTINUE
1020 OMM1 = OM
OM = FREQ(TTT + MOI - 1)
TESTOM = ABS(OM - OMM1)/OM
WRITE(NRES,*) 'MOI = ',MOI
flush(NRES)
WRITE(NRES,*) 'TEST = ',TESTOM,'FREQ = ',FREQ(TTT+MOI-1),'MLF = ',MLF(TTT+MOI-1)
flush(NRES)
close(NRES)
Fortran:
SUBROUTINE ZEIGSUB(A,B,N,FREQ,MLF,VR)
! ZGGEV Example Program Text
! NAG Copyright 2005.
! .. Parameters ..
INTEGER NIN, NOUT,NOUT2
PARAMETER (NIN=5,NOUT=6,NOUT2=7)
INTEGER NB, NMAX
PARAMETER (NB=64,NMAX=4)
INTEGER LDA, LDB, LDVR, LWORK
PARAMETER (LDA=NMAX,LDB=NMAX,LDVR=NMAX,LWORK=NMAX+NMAX*NB)
REAL(8), PARAMETER :: PI = 4.*datan(1.d+0)
! .. Local Scalars ..
DOUBLE PRECISION SMALL
INTEGER I, INFO, J, LWKOPT
integer, intent(in) :: N
! .. Local Arrays ..
COMPLEX *16,intent(in) :: A(LDA,NMAX)
COMPLEX *16,intent(out) :: VR(LDVR,NMAX)
COMPLEX *16 ALPHA(NMAX), B(LDB,NMAX),&
BETA(NMAX), DUMMY(1,1),&
WORK(LWORK), OM
DOUBLE PRECISION RWORK(8*NMAX)
double precision, intent(out) :: FREQ(NMAX),MLF(NMAX)
! .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
! .. External Subroutines ..
EXTERNAL ZGGEV
! .. Intrinsic Functions ..
INTRINSIC ABS
! .. Executable Statements ..
! OPEN(UNIT=NIN,FILE="/Users/ee/zggev-ex.d",STATUS="unknown")
open(unit=NOUT,file="/Users/ndilokelwaluis/Documents/Fortran/Vibration/MODEL_V1-0/fileout.dat",status='unknown')
open(unit=NOUT2,file="/Users/ndilokelwaluis/Documents/Fortran/Vibration/MODEL_V1-0/fileout2.dat",status='unknown')
WRITE (NOUT,*) 'ZGGEV Program Results'
! Skip heading in data file
! READ (NIN,*)
! READ (NIN,*) N
IF (N.LE.NMAX) THEN
!
! Read in the matrices A and B
!
! READ (NIN,*) ((A(I,J),J=1,N),I=1,N)
! READ (NIN,*) ((B(I,J),J=1,N),I=1,N)
!
! Solve the generalized eigenvalue problem
!
CALL ZGGEV('No left vectors','Vectors (right)',N,A,LDA,B,LDB,&
ALPHA,BETA,DUMMY,1,VR,LDVR,WORK,LWORK,RWORK,INFO)
!
IF (INFO.GT.0) THEN
WRITE (NOUT,*)
WRITE (NOUT,99999) 'Failure in ZGGEV. INFO =', INFO
ELSE
SMALL = DLAMCH('Sfmin')
DO 20 J = 1, N
WRITE (NOUT,*)
IF ((ABS(ALPHA(J)))*SMALL.GE.ABS(BETA(J))) THEN
WRITE (NOUT,99998) 'Eigenvalue(', J, ')',&
' is numerically infinite or undetermined',&
'ALPHA(', J, ') = ', ALPHA(J), ', BETA(', J, ') = ',&
BETA(J)
FREQ(J) = 0
MLF(J) = 0
ELSE
OM = ALPHA(J)/BETA(J)
WRITE (NOUT,99997) 'Eigenvalue(', J, ') = ',&
ALPHA(J)/BETA(J)
IF (DBLE(OM) .LE. 0) THEN
FREQ(J) = 0
MLF(J) = 0
ELSE
FREQ(J) = SQRT(DBLE(OM))/(2*PI)
MLF(J) = AIMAG(OM)/DBLE(OM)
END IF
END IF
WRITE (NOUT,*)
WRITE (NOUT,99996) 'Eigenvector(', J, ')', &
(VR(I,J),I=1,N)
20 CONTINUE
!
! call DSORT (FREQ, MLF, N, 2)
write(NOUT2,*) 'Ordered EIGENVALUES'
do 30 J = 1, N
write(NOUT2,*) 'EIG = ',J,', OM = ',FREQ(J),', MLF = ',MLF(J)
30 continue
LWKOPT = WORK(1)
IF (LWORK.LT.LWKOPT) THEN
WRITE (NOUT,*)
WRITE (NOUT,99995) 'Optimum workspace required = ', &
LWKOPT, 'Workspace provided = ', LWORK
END IF
END IF
ELSE
WRITE (NOUT,*)
WRITE (NOUT,*) 'NMAX too small'
END IF
STOP
!
99999 FORMAT (1X,A,I4)
99998 FORMAT (1X,A,I2,2A,/1X,2(A,I2,A,'(',1P,E11.4,',',1P,E11.4,')'))
99997 FORMAT (1X,A,I2,A,'(',1P,E11.4,',',1P,E11.4,')')
99996 FORMAT (1X,A,I2,A,/3(1X,'(',1P,E11.4,',',1P,E11.4,')',:))
99995 FORMAT (1X,A,I5,/1X,A,I5)
! close(NIN)
flush(NOUT)
flush(NOUT2)
close(NOUT)
close(NOUT2)
END SUBROUTINE ZEIGSUB
END PROGRAM TEST