Automatic differentiation of Fortran code, opinions?

FortranCalculus (FC) compiler is available and free. Download it at my website, goal-driven.net . Its roots start at NASA’s Apollo Space program that got us to the moon and back. The first language, developed around 1960 by TRW, was named Slang. It has 10+ solvers stored in a library that a user can call by name. These solvers used Automatic Differentiation (AD) to calculate the Jacobian and hessian matrices. For more, please read my profile.

I’m 77 years old and ready to retire. Anyone interested to taking on my website? My website has some 5+ Apps showing the power of FC and the FC compiler & manuals. The average download customer downloads ~4 Apps per day and the number of customers over the last year (2022) ranges from 50 to 90 downloads per customer! China is the number one country of interest. They are from 5 to 30 times the number of USA customers. They must be involved in Optimization.

2 Likes

There have been 3 AD based compilers that I’m a where of. NASA’s Apollo Space program had need for faster compilers as their engineers kept changing their minds (see the movie hidden figures?). So NASA hired TRW to write a new type compiler where solvers were stored in a library and required no development. Around 1960, TRW developed the first AD based compiler, Slang, that dropped solving problems to days instead of years. The first commercial AD compiler, PROSE, was made available in 1974 on CDC computers. Today there is a version named FortranCalculus that became available around 1990. They all 3 compilers handle Optimization problem nicely.

1 Like

Slang, PROSE, and FortranCalculus are AD based compilers.

1 Like

Hi Phil,

welcome to Discourse! I couldn’t find anything about AD-supporting compilers on the NASA Technical Reports Server (https://ntrs.nasa.gov/). Who do the initials TRW stand for?

Have you compared FortranCalculus with any AD frameworks in languages that are popular today, like Python, MATLAB or Julia? Do you still operate a company based around your (extended) Fortran compiler?

1 Like

The only option for Backward mode that I am aware of is Tapenade which I can recommend highly.

I see from the article (at arxiv: https://arxiv.org/pdf/1502.05767.pdf), page 24 of 43, the ADIFOR and NAGWare also work in reverse mode. In addition, a complete list of AD implementation in different languages is given.

1 Like

TRW = Thomson Ramo Wooldridge. Wikipedia will tell you more about it.

2 Likes

Thank you, I’ve found a paper from them:

The Q approach to problem solving (1969) | ACM

According to Google Scholar this paper only has 3 citations! Quite a shame. Too many new ideas at once.

1 Like

Hi @ivanpribec, TRW Systems, Inc. I believe its the full name, but sometimes I see it as just TRW, Inc. I believe the main developer was Joseph Thames who lived in the L.A., CA area, if that helps. Joe gave a SIAM conference presentation on PROSE in 1974 in Palo Alto (?), CA. The conference Heads couldn’t believe what he was saying and wouldn’t publish his presentation. Joe sent me a copy, but not sure where it is.

Have you compared FortranCalculus with any AD frameworks in languages that are popular today, like Python, MATLAB or Julia?

No, we are speaking about different worlds. FC / PROSE solve problems that the others have never though of. For example, nesting of solvers, implicit problems, etc. see my casebook for more. Most problems are solved in less than 50 lines of code. Try it, you’ll love it. Here is my Computer-Aided Designs presentation to High School students that may give you an idea about the FC language. It’s about 30 minutes long with 8 minutes of questions.

Do you still operate a company based around your (extended) Fortran compiler?

My company, Optimal Designs Enterprise (ODE), has just me. It took me some 22 years just to write my casebook, solve each problem, document code & solution. So, today I volunteer my time with middle & high school students via zoom who are thinking about a STEM career in industry … visit nepris.com and become a volunteer!

I do have one contract job I’m still thinking about. It’s regarding Oil Refiners and process control. Know anyone in such arenas? The problem is huge and may require building larger computers … about 10 times today’s computers; faster and more memory. It involves solving around 400 PDEs at the same time … nesting solvers. The problem may require about 20 hours per run, but must be done under 6 hours. Any interest? It will increase an Oil Refinery’s productivity by 10+%.

I think automatic differentiation should be one of the fundamentals of a heavy numeric language as Fortran. Could this be implemented in stdlib? @gardhor project seems like one of the more complete and usable projects for what I have seen (also some other dual numbers modules). But I think autodiff deserves a full representation in stdlib including both forward and backward modes

3 Likes

Thanks @fedebenelli for your nice comment !
Porting an existing AD library into the Fortran stdlib should be relatively easy, although it will be a lot of works. Personally, I don’t known how to start.
Right now, my library (AD_dnSVM) can be used only with forward mode.

1 Like

I have made a fork of Joshua Hodson’s DNAD where I have refactored it into a set of fypp macros/templates. I have also extended it to support hyper-dual numbers for computing Hessians, currently only for a small subset of the functions/oparators that the normal dual number type supports. The approach here is forward mode autodiff via operator overloading.

Using meta-programming for this has some advantages:

  • You can create several dual and/or hyper-dual types and give them names reflecting the design variables that their derivative is taken with respect to! Then the compiler will arrest you if you try to use them later in incompatible ways.
  • You can “inject” the function overloads in the same module where they are used to have them inlined without interprocedural optimization.
  • You can turn on/off “stringent” evaluation at the code generation level.

I am currently generating the code directly into the example, src and test folders and incuded them in the repo, so the examples and tests can be built and run directly (for instance fpm run --example main_example_hdual or fpm test) without having to run the code generator (fypp) first.

test/test_hdual_chain_mod.f90 (generated from test/test_hdual_chain_mod.fypp) shows an example of using two different dual number types together, one with design vector x and another with design vector y. A variable fy of type hdual_y_t can be converted into fx of type hdual_x_t by calling fx = chain_duals(fy, yx) where yx has type hdual_x_t and represents the y-design-vector.

Currently, the size of the design vector(s) need to be known at compile time. It should’nt be too difficult to lift this restriction and allow allocatable gradients and hessians. This is probably something that should be controlled at the code-generation level as well!

3 Likes

My compile, FortranCalculus (FC), uses Automatic Differentiation (AD) coupled with “operator overload” in order to compute the user’s Jacobian or Hessian matrix that is used by the chosen Solver. The matrix requires the system to compute the partial derivative of the entire user’s model routine and thus yields the Solver’s matrix. A true dynamic process! Every operator (+, -, X, /) requires another derivative. Your response " I think autodiff deserves a full representation in stdlib including both forward and backward modes" is impossible in a standard library do to the continual need for derivatives. Hope this gives you some understanding of FC.

NASA’s Apollo Space project hired TRW (System ?), Inc. to write/create a new computer compiler that would be easier & faster to get them to the moon and back. The result as the 1st Calculus (level) language, Slang (1950s). Slang require some 100 man-years of time to develop. Next was PROSE (1974), 1st commercial version, it required another 100 man-years of development time. FC is today’s version for Windows.

The Apollo program was in the 60s. The article by Wengert was published in 1964. The paper by McCully from TRW is dated to 1969, and his references mention what looks like an internal report - “Introduction to SLANG” - from TRW by Adamson, dated 1968.

Still insane to me that they wrote effectively a Fortran interpreter for auto-differentiation back then. At the moment the JAX and Julia programming language communities appear to be carrying the auto-diff torch.

2 Likes

Yes, autodiff with Julia is just another level thanks to it not being statically typed. It makes it easy to obtain high-order derivatives of any external function. This would be impossible in Fortran since it is statically typed.

Another thing is the recursive possibility it also brings (doing autodiff over an already autodiffed function). In my work I highly depend on obtaining Jacobian matrices of functions that are obtained through derivatives of other functions, I did draft an implementation of recursive dual numbers in Fortran for this but the performance was really bad. This made me consider moving to Julia at some point, but I still prefer the general Fortran “dev experience”, waiting for Julia’s jit to compile complex stuff or trying to do simple things such as debugging is just painful. Thanks to all the work with fpm. fortls and the vscode extension I can have a much better time on Fortran :slight_smile:

1 Like

The Apollo program was in the 60s.

My source on publications was Joseph Thames who is now deceased. Whether 1950s or 1960 is not clear, but have you seen the movie Hidden Figures, if you haven’t seen it, you should. It may help understand the communication problem between different types of people; e.g. NASA’s Engineers, Scientists, and Mathematicians. The movie shows engineers trying to solve some equations and getting no where fast! A mathematician comes along and solves the problem. But numbers say little to engineers. Near the movie’s end, the mathematician draws a graph showing a solution. The engineers finally get the ‘picture’ of what the equations are trying to say to them. Remember, a picture is worth a 1,000 words, right?

Still insane to me that they wrote effectively a Fortran interpreter for auto-differentiation back then. At the moment the JAX and Julia programming language communities appear to be carrying the auto-diff torch.

Good to hear! TRW’s work was done using Compass (assembly language) and Fortran (?). Switching to PROSE & TSPROSE still was in the same languages. Thames was the first to believe there was a way to use Fortran completely, no assembly!

Slang, PROSE, nor FortranCalculus are Fortran interpreter languages; All used Operator Overload techniques (some call it “symbolic differentiation at a point”) to take derivatives on the fly.

The following link is for match-n-freq that is an example with source code (*.fc files). This App has nested solved (2 or 3 at times). Install & Run a few Demos to get the idea what’s going on.

Here is a chart mentioning what key issues that are going on:
MetaScience Modeling

Another just showing the key Find stmt.

Hope these give you an idea of what “symbolic differentiation at a point” can do for you.

This comment is more on Fortran than AD. As languages increase in accuracy, all constants need more accuracy. PI is a great example! Try pi = 4 * atan( 1.) that yields to most possible digits.

I once got docked 10 points on a math exam where the answer was PI because I tried to write out PI to as many decimal points as I could remember (being an engineer and not a mathematician). The teacher explained I didn’t give the correct answer because the symbol for PI represents the number to an infinite number of digits and I only gave the answer to 6 digits. He wanted the symbol for PI and not numbers.

pi = 4 * atan(1d0) gives more correct digits than pi = 4 * atan( 1.), and most compilers offer at least one real type with even higher precision.

Good to know, thanks!

In expressions like this where the intention is to write portable code, the programmer should never use literal constants to specify the kind values (or the precision). A better way would be something like

   integer, parameter :: wp = kind(pi)
   ...
   pi = 4.0_wp * atan( 1.0_wp )

or one of the other more or less equivalent ways to write the expression. This way, when the programmer wants to change precision (either more or less), the expression automatically adapts accordingly.