Old Fortran Program

[JohnNichols]

Steve:

I came across an old Fortran program from UCLA.  It is a very useful program and the site provides a copy of the program complied in Compaq Fortran in 32 bit.

I tried unsuccessfully yesterday to compile the program.  I broke it down into reasonable units and slowly worked to fix some things.

But when I get to the FEM module - it starts and never stops on the compiler.  I stepped through it by adding back a line at time, but could not spot the error.  If I did not have a working EXE I would wonder, but there is a working exe in Windows.

Any ideas, I enclose the original code and the modifications I made. It is an interesting program that has the potential to save lives in earthquakes so I am keen to get it working.  If I cannot I will rewrite the algorithms in C#. 

JMN

[TimP]

F90 introduced the first reasonable alternative to statement functions.  Up to about 10 years ago, statement functions often had better in-line optimization than alternatives.

About the time compilers began fixing limitations of internal subprograms, they also started to drop the limitations on how the old statement functions could be combined with f90 syntax (according to restrictions placed on them in the standard).

I put in the caution about using statement functions in accordance with f77 since much usage of them prior to f77 involved inadvertent use of extensions of various compilers (some of those reflected in the internal procedures of f90).

[JohnNichols]

Tim:

I spent an enjoyable morning inside the program.  Once I had been taught that it would compile as release it was not that hard, just had to go back to putting write statements at the error points

This reminded me of Microsoft Fortran 3.3.1 and the old Compaq portable. It was nice in those days, it took 40 minutes to compile a program so you had lots of free time to do things, Now with the flash drive this thing runs in no time, so you only leave your seat to get another coke.

Thanks for the note on inline functions, I have programmed since 1978, but there is no way I would use such a beast.  it appears to be working.

The problems this morning related to DATA statements not being read properly, my fault. 

I am now at the stage of sorting out the Stiffness matrix and turning it from a vector to an array.

JMN

[mecej4]

John Nichols: Attached is a zip archive containing two files: the source file from http://peterbird.name/oldFTP/neotec/SHELLS/SHELLS.f90.txt, modified just enough to allow us to replace the calls to IMSL routines by equivalent calls to MKL routines, and an input file that I created.

I tested the IMSL and MKL versions using Step-7 of the instructions at http://peterbird.name/oldFTP/neotec/SHELLS/Earth5/Earth5-049/READ_ME.txt. To compile the source files, do ifort /Qmkl /fpscomp:filesfromcmd shellsm.f OR (ifort /Qimsl /fpscomp:filesfromcmd shells.f, for the unmodified version). Download the input data files named in the instructions to the working directory on your machine. Then, at the command line, do shellsm < shells.inp. The program works with either the 32-bit or 64-bit versions of IFort.

The IMSL and MKL versions ran to completion and gave similar results, but I have no feel for what to expect for output. Using the multi-threaded MKL routines, the program runs about ten times faster with eight threads (4 core i7-2720 CPU with hyper-threading) as compared to the single-threaded IMSL version.

[JohnNichols]

Dear MECEJ4

Thank you for your kind work on this project. I have not used the Fortran compiler directly since the 1990's so it did not cross my mind. It is a brilliant suggestion.

I have never played with the MKL libraries before this program. I spent yesterday afternoon looking at the documentation for the MKL libraries and had thought I had found the packed solver in the MKL library would replace the ISML library. I went through the parameters and had sorted out how to make the call.

Your kind email landed about that time and looking inside your code, I noticed that you had found the error in the packing code that one variable was ignored in the transform equation. I had printed out results and the variable always had a value of zero, so I had ignored that issue, but thanks for the fix. 

As you kindly showed it works nicely with ifort.

I would also like to learn why I am having so much trouble with MKL. I have set the 32 bit environment variables, turn on the release mode and set the box inside the program properties to INCLUDE MKL libraries.

The program compiles, but does not find the MKL.

I am running mkvars in the VS command prompt window. But when I check the environment variables there is no sign of the added MKL variables. 

Am I supposed to make all the link connections inside the program properties instead of taking the shortcut step. 

I prefer not to go fooling with the system variables.

There is no problem using the IFORT, but not getting the other way going knaws at my mind. ( I realize knaw is obsolete, but so am I in this modern world.)

Professor Bird from UCLA has developed an excellent idea. He has split the world's surface into essentially triangular shapes, he joins them with appropriate boundary conditions and then solves the FEM mesh he has created to determine the stress and strain tensors. It is a major undertaking and he is to be congratulated on doing such an excellent job on such a massive task. This code is 12000 lines and it is well commented.

He also has a program to generate synthetic earthquake catalogs. ie estimates of future seismicity (SEC). This program also uses the packed solver, so I will go straight to the IFORT solution. SEC takes as one of its input files an output file from shells. I am looking to slowly vary the variables within the SHELLS program to generate a set of solutions and then use them in the SEC to get a huge range of catalogs. From that I will look to estimates of the death toll.

So thank you for all the help.

While I have been working on this project I stumbled across a plan of a pipe crusher from IOWA State Uni in a 1913 book. YOu can see the PDF of the original as an attachment. This was state of the art in 1913 at the best AG Eng school in the world. We have come a long way in a century.

JMN

JMN

[dboggs]

John, regarding your comments on storing a square matrix in a vector:

This is indeed common, especially if the matrix is tightly banded or sparse. A variety of techniques has been developed to take advantage of the extreme memory storage savings that is possible. One that I learned many years ago is called "skyline" storage, but even this is old technology and better ones probably exist. In general, they all aim to extract only the populated elements of the matrix and then "line them up" in a vector with appropriate pointers to maintain their original matrix locations.

Is your problem, by any chance, related to finite element analysis (FEA)? If so the above is very likely.

[mecej4]

JMN: I did not use your modified version of the PeterBird source file, since I wanted to start with the reference source file, which the authors have documented quite well. In substituting the IMSL calls by MKL/Lapack calls, all that I need to do in addition was to replace the INDEXK statement function declarations to 

       indexk(irow, jcolum) = (jcolum-1)*lda + nlca + nuca + irow - jcolum + 1

To make this work, I added NLCA to the blank common, and set LDA = NLB+NLB+1+NUB in Subroutine KSIZE.

Creating a project for IFort using MKL is easy. Start VS, create a new Intel Fortran console application, and add the single Fortran file to the project sources. Next, in the project properties, go to Configuration Properties-Fortran-Libraries and select one of the options for "Use Intel Math Kernel Library". If you have a working IFort/MKL installation, you should be able to build the project with no special Fortran options.

The (subjective) impression that I have gathered about the SHELLS program is that it is high quality. One measure of the quality is that I was able to make the conversion from IMSL to MKL with little effort. I have built and run the program with three different compilers, 32-bit and 64-bit, and they all yielded similar results. We have to be careful in judging such codes that we do not apply today's standards to code from fifteen years earlier.

[ASIDE:]By "knaw" do you mean "gnaw" but by a toothless antagonist?

[JohnNichols]

Steve:

----------------------------------------------------------------------------------------------------------------------------------------------

Do I have to do something like this to get the Shells program working nicely with the MKL routine

ie use an OBJCOMMENT?

!call dgbsv(n,nlca,nuca,nrhs,A,ldab,ipiv,B,ldb,info)

Because I can not get them to run nicely at the moment.

The LAPACK program works nicely once I add in the statement below. 

------------------------------------------------------------------------------------

John - add the line:

!DEC$ OBJCOMMENT LIB:"mkl_lapack95.lib"

to your source (just after the PROGRAM statement is fine.) The LAPACK routines are not linked in by default.

[Steven_L_Intel1]

John, the OBJCOMMENT directive is used to solve linking problems. I have no idea what you mean by "run nicely".

[JohnNichols]

Dear MECEJ4:

Professor Bird has a vast set of programs all in Fortran.  I will slowly work my way through them working out how to use them to solve my problems.  I was looking at the MKL problem with the VS version of shells, so I took some time and read the ISML, MKL and Numerical Recipes in Fortran stuff on the packed routines.  (I would not be surprised if they all share a common ancestry, although it would be nice if they used the same symbols. So I added the F90 Numerical Recipes in Fortran routines into shell and failed. Going back one step I got the F77 recipes code to work for the first stage of the dgbs program.

The shells program iterates about 50 times on an array that in real terms is 32000 by 32000.  In the MKL version an iteration is a few minutes, in the ISML it is a lot longer, but in the F77 version the following code :

  do 18 k=1,n
        write(*,*)k,n
        dum=a(k,1)
        i=k
        if(l.lt.n)l=l+1
        do 14 j=k+1,l
          if(abs(a(j,1)).gt.abs(dum))then
            dum=a(j,1)
            i=j
          endif
14      continue
        indx(k)=i
        if(dum.eq.0.) a(k,1)=TINY
        if(i.ne.k)then
          d=-d
          do 15 j=1,mm
            dum=a(k,j)
            a(k,j)=a(i,j)
            a(i,j)=dum
15        continue
        endif
        do 17 i=k+1,l
          dum=a(i,1)/a(k,1)
          al(k,i-k)=dum
          do 16 j=2,mm
            a(i,j-1)=a(i,j)-dum*a(k,j)
16        continue
          a(i,mm)=0.
17      continue
18    continue
       

takes 4.5 hours for a single iteration.

MKL is way faster, but not nice to link in.

JMN

[JohnNichols]

Steve:

I cannot get the shells program to execute after compiling within VS 2013 using the MKL  libraries. it executes after being compiled with IFORT and /QMKL etc from a command prompt.

This is a brand new installation of Windows 7.1 on a Dell Intel Core 7 processor, with lots of ram and a 1 TB flash drive.

VS 2013 Professional was installed two weeks ago, with all the updates needed, the Latest Fortran was installed along with MKL after VS. There are no other VS studios or Intel Fortran

I have  used the VARBATS command and when I develop the program properties, they are in release mode with the MKL Link option selected.

if they compile, then when the program runs it crashed VS or gives me a MKL___ .lib.dll not found error.

How do I take an old Fortran program - create a new VS solution, add a call to a MKL routine mainly GBSV for example so that it will run when it is compiled. 

Shells is a pretty straightforward program, the implicit none problem aside, it does run in IFORT and gives answers that look reasonable.

MKL is significantly faster than ISML for this routine and much much faster than the Numerical Recipes in Fortran.  

Clearly I am missing something simple.

can you please try SHELLSM.f. It compiles with IFORT but if you can get it to compile on VS and even start I will be very happy.

You do not need input files, it will not even start.

You might try in debug mode as it will not even compile or tell you the offending code.

Am I missing some link or DEC commands?

That is what I mean by not run nicely.

JMN

[Steven_L_Intel1]

So the problem is the run-time error about the missing DLL? That's because the path to the MKL DLLs aren't added to PATH by default.

Right click on "Computer" on your desktop and select Properties. Click on Advanced System Settings, then Environment Variables on the Advanced tab. Under System Variables, scroll down to Path. Select that and click Edit.

Insert at the end:

;%IFORT_COMPILER14%redist\intel64\mkl;%IFORT_COMPILER14%\redist\ia32\mkl

Click OK, OK, OK.

That should remove the DLL error. You may need to log out and in again for that to take effect.

I can see the problem with the large Shell source not compiling and will send that on to the developers.

[jimdempseyatthecove]

Consider this:

  do 918 k=1,n
        write(*,*)k,n
        smallestMag(k)=abs(a(k,1)) ! smallestMag(n)
        i=k
        if(l.lt.n)l=l+1
        do 14 j=k+1,l
          if(abs(a(j,1)).gt.smallestMag(k))then
            smallestMag(k)=abs(a(j,1))
            i=j
          endif
14      continue
        indx(k)=i
918   continue
  do 18 k=1,n
        if(smallestMag(k).eq.0.) a(k,1)=TINY
        if(i.ne.k)then
          d=-d
          do 15 j=1,mm
            dum=a(k,j)
            a(k,j)=a(i,j)
            a(i,j)=dum
15        continue
        endif
        do 17 i=k+1,l
          dum=a(i,1)/a(k,1)
          al(k,i-k)=dum
          do 16 j=2,mm
            a(i,j-1)=a(i,j)-dum*a(k,j)
16        continue
          a(i,mm)=0.
17      continue
18    continue
       

Jim Dempsey

[JohnNichols]

Jim:

I will run it tonight and let you know. 

The problem with the packed array, is that it is almost all still zeros. The bound limit is 1400 for 32000 nodes, but on most lines there are only three triangles, ie triangle 1 can touch say 2,3 and 4 but not the other 10,000 triangles. (I printed out a few sets of the rows).

it would be a lot more efficient to pack a duple with the (I,J) and a vector with the value from A. Should be able to do it with a lot less entries. Why play with zeros at all?

In the matrix manipulation you never play with more than one value or two values at a time, so a lookup of an ordered I,J has got to be a lot faster.

Perhaps I am missing a point.

JMN

[mecej4]

John Nichols wrote:

Why play with zeros at all?

In the matrix manipulation you never play with more than one value or two values at a time, so a lookup of an ordered I,J has got to be a lot faster.

Perhaps I am missing a point.

Matrix entries within the band that are zero may get filled in when the L-U (or L-D-U) factors are computed. In other words, after factorization the bandwidth is preserved, but zero entries within the band do not necessarily yield zero entries in the factors.

[JohnNichols]

True:

but I was thinking that if we used a second working vector we could invert without needing a large matrix. I suppose it depends on the density of points in the matrix, before and after the manipulation.

Just a thought.

JMN

[jimdempseyatthecove]

The point of using two do k loops is to increase cache locality (at least L2), .AND. should the 32000x32000 plus other data not fit entirely in RAM, then decrease then number of page faults.

Note, with cache line of 64 bytes, 16 floats fit in cache line. The first loop is reading down a column with a stride of 32000. Assuming no data in cache to begin with the first column read will encounter 32000 (cache lines) requiring 32000 x 64 bytes of cache, which is 2MB, which generally fits in L3 cache.

Running the code will verify if it helps or not. It is a simple enough change to test.

Jim Dempsey

[mecej4]

John Nichols, #30 wrote:
(emphasis added)

The shells program iterates about 50 times on an array that in real terms is 32000 by 32000.  In the MKL version an iteration is a few minutes, in the ISML it is a lot longer, but in the F77 version the following code :

<<--CUT-->>       

takes 4.5 hours for a single iteration.

MKL is way faster, but not nice to link in.

JMN

I don't know which problem you ran and obtained results, but when I ran the MKL version SHELLSM.f on a laptop with an i7-2720QM CPU, with 4 cores X 2 hyperthreads/core= 8 threads, on the problem data listed in #24, the matrix is 32016 X 32016, with upper and lower bandwidth = 1451, and the program took only 8.5 seconds per iteration.

The same program when run on a tablet/laptop with an i5-4200U CPU (2 cores, 4 threads) takes 9.5 seconds/iteration, so you should be able to get decent throughput on a properly outfitted engineering workstation with a multi-core CPU.

[JohnNichols]

I willrun it.

[JohnNichols]

I got part of the changes running and it was way faster.

As it finished the code gave me an access violation issue and I was trying to sort it out. 7 year old got in the road.

I ran the Shells program from the IFORT build and got 45 seconds per iteration. If you are getting 9 seconds per run, then I am not running on separate processors, but you are.

MKL and VS are a dud together and it is time to forget this combination. Pity.

IFORT and MKL are the way to go. (Like going back to the 80s)

I have two machines, a Dell Precision which I like because I can take it anywhere and a special build using an ASUS motherboard and a lot of processing power, but it is not portable so it plays music most of the day. The IT guys say it is the fastest in the building. They like it, because when I upgrade it they got the parts, that are better than what they normally have. I believe I have the only special build in the building.

it is time to move onto the next program in the shell set.

Thanks

JMN

[Steven_L_Intel1]

MKL works fine from VS. You probably need to add the MKL DLLs to PATH as I suggested earlier. Note that if you build from the command prompt, the programs will run ok as long as you start them from that command prompt session, but not if you run them independently (without updating PATH.)

[mecej4]

John Nichols wrote:

MKL and VS are a dud together and it is time to forget this combination. Pity.

IFORT and MKL are the way to go. (Like going back to the 80s)

I think that you are laboring under a couple of misconceptions.

VS does not include a Fortran compiler, libraries or runtime DLLs. It gives you a language-sensitive and customizable editor and a GUI front-end for using the same compiler and libraries as you use when building at the command line. Compiler options that you specify in the compiler invocation from the command line have equivalent configuration / menu choices. When you say "Build" from the VS IDE, what happens is the same sequence of events that happens when you compile from the command line. Think of VS as an agent which collects all your desired options, and issues the correct command on your behalf, in the correct command line syntax (which, you say, is "like going back to the 80s"). In a sense, the IDE lets you issue commands without being caught in the act of committing decadent deeds.

As Steve Lionel has stated, it is quite easy to use MKL from VS or the command line for a project that consists entirely of Fortran sources and the simplest way of using MKL is sufficient.