PCA - Fortran IV Program

JohnNichols · ‎10-10-2021

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

jimdempseyatthecove · ‎10-10-2021

** UNTESTED **

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

Jim Dempsey

cean · ‎10-11-2021

confirm.

mecej4 · ‎10-11-2021

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.

JohnNichols · ‎10-11-2021

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.

mecej4 · ‎10-11-2021

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.

JohnNichols · ‎10-11-2021

I agree.

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

JohnNichols · ‎10-11-2021

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.

JohnNichols · ‎10-11-2021

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.

JohnNichols · ‎10-11-2021

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.

Steve_Lionel · ‎10-11-2021

@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.

JohnNichols · ‎10-11-2021

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

JohnNichols · ‎10-11-2021

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.

cryptogram · ‎10-13-2021

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

JohnNichols · ‎10-14-2021

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.

andrew_4619 · ‎10-11-2021

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

* is any number, std fortran (2003?)

Steve_Lionel · ‎10-12-2021

Unlimited group repeat is an F2008 feature.

JohnNichols · ‎10-12-2021

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.

Steve_Lionel · ‎10-13-2021

@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.

JohnNichols · ‎10-13-2021

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

JohnNichols · ‎10-13-2021

The letter T from the font file in decimal.

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