I was playing with an old Fortran IV program - PCA from Wahlstedt and Davis - Kansas Geological Survey, they have an equivalence statement at the beginning, but the program does not work with this statement in.
has the meaning of equivalence changed. I have checked the results with new data.
This is the code as I currently have it:
! PCA.f90
!
! FUNCTIONS:
! PCA - Entry point of console application.
!
!****************************************************************************
!
! PROGRAM: PCA
!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************
program PCA
use Base
implicit none
Integer M,N
LOGICAL TEST, OPDAT, THERE
CHARACTER*10 LABEL(mn,4)
Character*10 Title(20),FMT(mn_12)
Integer manual
Integer astat
Integer L
Integer IVX
Integer i
Integer J
Integer k
Integer EN
REAL (KIND=dp) :: X(ml,nn), SSD(mn,mn)
REAL (KIND=dp) :: SX(mn), SD(mn),PERCON(mn)
REAL (KIND=dp) :: SS(mn,mn), R(mn,mn),COVMAT(mn,mn)
COMMON / WJ/ M,N, X,SSD,Title
!EQUIVALENCE (SSD(1,1), SS(1,1),SX(1), PERCON(1))
inquire(file = "data.in", EXIST = THERE)
if(there) then
open(8,file="data.in",status = "UNKNOWN", IOSTAT = astat)
else
stop 'NO file'
end if
500 READ (8, *, ERR = 5000, END = 5000) M, MANUAL,OPDAT
write(*,5010) M, MANUAL, OPDAT
5010 Format(" Data Reflection ::",/," M :: ",I4,/," Manual :: ",I4,/," OP DATA :: ",L4)
IF(M .le. 0) then
CALL EXITM()
end if
do 103 L = 1,M
read(8,*, ERR = 5000, END = 5000)(LABEL(L,IVX),IVX = 1,4)
103 end do
do 100 i = 1,M
SX(i) = ZERO
do 201 j = 1,M
SS(i,j) = ZERO
201 end do
100 end do
READ (8, 10) TITLE
10 Format(20A4)
READ (8,10)FMT
Write(*,1000)TITLE
1000 Format(9x,20A4,/)
SX = ZERO
SS = ZERO
N = 0
do 16 n=1,300
read(8,*)(x(n,j),j=1,M),TEST
write(*,310)n,(x(n,j),j=1,M)
310 Format(" ",I4,5(" ",F8.2))
do 20 i=1,M
SX(i) = SX(i)+X(N,i)
do 25 j = 1,M
SS(I,j) = SS(i,j)+X(n,i)*X(n,j)
25 end do
20 end do
write(*,310)N,(SX(i),i=1,M)
if(TEST) then
goto 26
end if
16 end do
write(*,*)"To many data points"
call exitM()
26 write(*,27)N
27 Format(//" Number of data points :: ", I4)
If(OPDAT) THEN
go to 1643
end if
write(*,1002)
1002 Format(//"------------------------------------------------------------",//," Input Data",/)
do 101 i = 1,N
write(*,1003) I
1003 Format(" Sample Number :: ",i4)
101 End do
1643 EN = N
do 30 i = 1,M
do 30 j = 1,M
SSD(I,J) = (SS(i,j)-(sx(i)*SX(J))/float(EN))/(float(EN) - 1.0D0)
COVMAT(I,J) = SSD(I,J)
30 end do
Do 35 I = 1,M
SD(I) = SQRT(SSD(I,I))
SX(I) = SX(I)/FLOAT(EN)
35 end do
Write(*,1000)Title
Write(*,1005)M,N
1005 Format(" Number of variables :: ",I4,/," Number of data points :: ", I4,///," Variable Mean Standard Deviation")
do 102 I =1,M
write(*,1001)i,SX(I),SD(I)
102 end do
1001 Format(10x,I4,2(5x,f15.6))
do 60 i=1,M
do 61 j = 1,M
R(I,J) = SSD(I,J)/(SD(I)*SD(J))
61 end do
60 end do
if(MANUAL .eq. 2) then
write(*,*)" You have selected the Manual = 2 option"
goto 72
end if
72 if(MANUAL .eq. 1) then
write(*,*)" You have selected the Manual = 1 option"
goto 500
end if
write(*,73)
73 Format(//" You have selected the Manual = 3 option",/)
stop 'End of Program'
5000 Stop 'File input error'
end program PCA
subroutine exitM()
implicit none
STOP " EXIT SUBROUTINE"
return
end subroutine
Module Base
INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
INTEGER, PARAMETER :: sw = 2, srA = 1, st = 11, sa = 3, sm = 4, sc = 12, sRT = 25
Integer, parameter :: slog = 9, sOUT = 10, sCAD = 21, sINI = 13, sWULF = 18
Integer, parameter :: array_size = 3000
integer, parameter :: nt = 2000
integer, parameter :: mt = 2000 ! Number of members
integer, parameter :: mn = 30
integer, parameter :: ml = 300
integer, parameter :: mn_12 = 12
integer, parameter :: mn_4 = 4
integer, parameter :: count1 = 20000
integer, parameter :: nn = 5
integer :: readFlag = 0
REAL (KIND=dp) :: g = 9.806, pi = 3.14159265D0
REAL (KIND=dp) :: delta = 0.000000000001d0
REAL (KIND=dp) :: h = 0.0005d0
REAL (KIND=dp) :: alpha = 0.16666d0
REAL (KIND=dp) :: beta = 0.25d0
REAL (KIND=dp) :: gamma = 0.5d0
REAL (KIND=dp) :: ONE = 1.0d0
REAL (KIND=dp) :: SMALLMASS = 0.001d0
REAL (KIND=dp) :: TWO = 2.0d0
REAL (KIND=dp) :: ZERO = 0.0d0
REAL (KIND=dp) :: DAMPRATIO = 100.0d0
REAL (KIND=dp) :: RATIO1 = 0.45d0
REAL (KIND=dp) :: MULTIPLER = 0.4d0
REAL (KIND=dp) :: RATIO2 = 0.25d0
REAL (KIND=dp) :: SuggestG = 0.0d0
contains
end module Base
