- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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
```

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

** UNTESTED **

Missing parens

EQUIVALENCE (SSD(1,1), SS(1,1)), (SX(1), PERCON(1))

Jim Dempsey

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

confirm.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

The posted code appears to be an incomplete conversion to Fortran 9X of the original Fortran IV code. All that the converted code appears to do is to read data from a file, form the covariance and correlation matrices from the data, and print those matrices out.

There are no subroutines or external functions, so I wonder why a common block is declared. Then comes the question of why COMMON and EQUIVALENCE, both of which are obsolescent and troublesome, are used. The block WX is declared in only one place, so it may be left out altogether. Cursory inspection of the code that is posted above indicates that the EQUIVALENCE can be left out without affecting the output of the program. Unless there is a well-reasoned justification for using EQUIVALENCE, avoid it.

MKL VSL contains summary statistics routines. Consider just reading the data in your code and then calling these high performance MKL routines to obtain the desired correlation matrix.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I stumbled across this old Fortran IV code, when I was looking at some PCA stuff. It is like a hobby, I enjoy just playing with the old code and seeing what they did with limited resources. It was a pleasant Sunday of just playing whilst we also discussed 3D.

Actually they had one good idea for data entry

```
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,4(" ",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)
```

, place a true or false value at the end of each data line. A true signals you have reached the end. In terms of a punch card deck it is not a bad idea.

The IF statements are like a crossword puzzle, just more interesting. Things like this in a crossword puzzle drive me crazy:

Clue :: designed by committee...difficult to use, -- three letters. By the way if you know the answer you have read to much computing stuff.

The working program is included in the zip file. It gets to the stage of the eigenvectors.

I agree we use the correct MKL subroutines, but this was just a fun day. I tested the results against a standard library and they are pretty good.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Thanks for posting the source and data files. I can confirm that commenting out the COMMON and EQUIVALENCE statements has no effect on the output of the program, just as I speculated in my previous post in this thread.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I agree.

It was just some harmless fun, better than watching a movie.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

How did the GE625 computer on which the original compiler program was run, cope with the author's use of HDTEST as a real, without implicit none, the program would never compile on a modern compiler. There are several such instances - was there automatic progression from real to int and back?

The original program has the correct answers I checked the eigenvalues with a modern subroutine.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I also do not understand their reliance on arithmetic if's, with minor thought they could reduce the code complexity a lot.

Harrison in '73 had a similar problem with his structures program.

Anyone understand this over reliance?

Also the code as posted is not the final code, there are minor logic errors that I fixed to make it work.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

And the other thing of interest that we were talking about last week was the use of character variable FORMAT statements instead of line numbers. These programmers actually read in the format statement for the data as a character array and then used it in the read statement.

Very clever.

It took me a minute to work out and then I smiled.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

@JohnNichols wrote:

And the other thing of interest that we were talking about last week was the use of character variable FORMAT statements instead of line numbers. These programmers actually read in the format statement for the data as a character array and then used it in the read statement.

That's still a feature of the language.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Actually it is a really good idea.

I am finishing the program and I found the location where the common is used.

One subroutine uses the common, everywhere else they passed the variables. Easy fix.

The other annoying feature, they will use a lot of IPIV, JPIV variables with long names, when I,J,K,L etc would do fine.

In the subroutine with the common, someone initializes an array to zero. This is a first.

My guess written by different authors, as the numbering suddenly went to 1, 2, 3 on the do loops from 100 etc..

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I finished the program.

`310 Format(" ",4(" ",F8.2))`

I was wondering if one could place the format files in the input data as Steve suggested, but you could change the format for the number of elements in an array for instance so have

`310 Format(" ",5(" ",F8.2))`

I have not seen code that drew for the page printer, such as a daisy wheel etc. This code does - it was interesting to see how they did it.

Looks like ::

There are mistakes in the plotting routine, I need to check all the placements and the changes from CHAR to INT. But it was interesting learning experience.

Complete program is attached.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Not sure what the question is about HDTEST. Starting with H, it would have been real by default.

As for the use of arithmetic IF, my speculation is that the subroutine was written earlier and reused, and was never converted

from FORTRAN II, where the arithmetic IF was the only option. Brings back memories of working with the IBM 1620

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

The question about HDTEST is an error on my part. I have used implicit none for so long, I always thought that integers started at H. Not sure why. You learn something in your old age.

No, it is quite clear several people wrote this deck of cards. They had different styles in numbering their format statements, it is like an FFT signal for an earthquake.

These people used a GE, which was GE's attempt to stop using IBM main frames. They said a 100 variable run took about 30 seconds, now it takes a few depending on the write statements.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

310 Format(" ",*(" ",F8.2))

* is any number, std fortran (2003?)

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

Does the star * mean unlimited repeat?

```
0x18,0x20,0x20,0x5F,
0x00,0x00,0x00,0x00,0x00,
```

If I have a file with 9124 of these numbers, how do I read them and then convert to pure binary.

The file represents a FONT type in 24 by 32 blocks, I now have a program that will print to the bitmap in this block size, but I have no desire to design fonts.

If you are asking why - it is easy to play with such files and develop new graph types without having to play with a library.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

@JohnNichols wrote:

Does the star * mean unlimited repeat?

Yes. It means that the parenthesized format group will repeat as many times as is necessary to satisfy the I/O list. Think of it as specifying a VERY large repeat count.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

I thought that writing a routine for reading in the hexadecimal values and do a decimal translation would be a less than simple task. But it is actually quite easy .

```
subroutine rgbimage_fonts(this, sr)
implicit none
class(rgbimage), intent(inout) :: this
character*20 startLine(3)
character*4 header(4)
character*4 line(96)
character*1 A, B
integer sr
integer i
Integer J
integer k
Integer l
type(HEXADECIMAL) :: FONT(2)
write(*,*)"here I am"
read(slog,*)(startline(i),I=1,3)
write(*,100)(startline(i),i=1,3)
100 format((A))
read(slog,*)(header(j),j=1,4)
do 120 k = 1,2
read(slog,*)(line(j),j=1,96)
write(*,*)(line(j),j=1,96)
do 130 l = 1, 96
font(k)%dataH(L) = line(l)
font(k)%dataS(L) = font(k)%dataH(L)(3:4)
A = font(k)%dataS(L)(1:1)
B = font(k)%dataS(L)(2:2)
call numH(A,B,font(k)%num(L))
130 end do
120 end do
return
end subroutine rgbimage_fonts
subroutine numH(A,B,num)
implicit none
character*1, intent(IN) :: A,B
integer, intent(inout) :: num
integer aNUM, bNUM
integer i
aNUM = 0
bNUM = 0
call numReturn(A,anum)
call numReturn(B,bNUM)
num = anum*16 + bnum
if(num .gt. 0) then
write(*,100)A, B, num
100 Format(" A :: ",A,/," B :: ",A,/," Number :: ",i2)
end if
return
end subroutine numH
subroutine numReturn(C, num)
implicit none
character*1, intent(IN) :: C
integer, intent(inout) :: num
select case (C)
case("0")
num = 0
case("1")
num = 1
case("2")
num = 2
case("3")
num = 4
case("4")
num = 4
case("5")
num = 5
case("6")
num = 6
case("7")
num = 7
case("8")
num = 8
case("9")
num = 9
case("A")
num = 10
case("B")
num = 11
case("C")
num = 12
case("D")
num = 13
case("E")
num = 14
case("F")
num = 15
case default
num = 0
end select
return
end subroutine numReturn
```

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

```
0 0 0
15 255 240
15 255 240
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 24 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
```

The letter T from the font file in decimal.

Somewhere I read that you use little gylphs to design your fonts. -- Interesting in Fortran.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page