I'll add to this answer if I find anything more detailed, but I have a place to start looking for the arguments to GJR
.
This function is part of the Sperry UNIVAC MATH-PACK library - a full list of functions in the library can be found in http://www.dtic.mil/dtic/tr/fulltext/u2/a170611.pdf GJR
is described as "determinant; inverse; solution of simultaneous equations". Marginally helpful.
A better description comes from http://nvlpubs.nist.gov/nistpubs/jres/74B/jresv74Bn4p251_A1b.pdf
A FORTRAN subroutine, one of the Univac 1108 Math Pack programs,
available on the library tapes at the University of Maryland computing
center. It solves simultaneous equations, computes a determinant, or
inverts a matrix or any combination of the three above by using a
Gauss-Jordan elimination technique with column pivoting.
This is slightly more useful, but what we really want is "MATH-PACK, Programmer Reference", UP-7542 Rev. 1 from Sperry-UNIVAC (Unisys) I find a lot of references to this document but no full-text PDF of the document itself.
I'd take a look at the arguments in the function call, how they are set up and how the results are used, then look for equivalent routines in LAPACK or BLAS. See http://www.netlib.org/lapack/
I have a few books on piping networks including "Analysis of Flow in Pipe Networks" by Jeppson (same author as in the original PDF hosted by USU) https://books.google.com/books/about/Analysis_of_flow_in_pipe_networks.html?id=peZSAAAAMAAJ - I'll see if I can dig that up. The book may have a more portable matrix solver than the proprietary Sperry-UNIVAC library.
Update:
From p. 41 of http://ngds.egi.utah.edu/files/GL04099/GL04099_1.pdf I found documentation for the CGJR
function, the complex version of GJR
from the same library. It is likely the only difference in the arguments is variable type (COMPLEX
instead of REAL
):
CGJR
is a subroutine which solves simultaneous equations, computes a determinant, inverts a matrix, or does any combination of these three operations, by using a Gauss-Jordan elimination technique with column pivoting.
The procedure for using CGJR
is as follows:
Calling statement: CALL CGJR(A,NC,NR,N,MC,$K,JC,V)
where
A
is the matrix whose inverse or determinant is to be determined. If simultaneous equations are solved, the last MC-N
columns of the matrix are the constant vectors of the equations to be solved. On output, if the inverse is computed, it is stored in the first N columns of A
. If simultaneous equations are solved, the last MC-N
columns contain the solution vectors. A
is a complex array.
NC
is an integer representing the maximum number of columns of the array A
.
NR
is an integer representing the maximum number of rows of the array A
.
N
is an integer representing the number of rows of the array A
to be operated on.
MC
is the number of columns of the array A
, representing the coefficient matrix if simultaneous equations are being solved; otherwise it is a dummy variable.
K
is a statement number in the calling program to which control is returned if an overflow or singularity is detected.
1) If an overflow is detected, JC(1)
is set to the negative of the last correctly completed row of the reduction and control is then returned to statement number K
in the calling program.
2) If a singularity is detected, JC(1)
is set to the number of the last correctly completed row, and V
is set to (0.,0.)
if the determinant was to be computed. Control is then returned to statement number K
in the calling program.
JC
is a one dimensional permutation array of N
elements which is used for permuting the rows and columns of A
if an inverse is being computed .. If an inverse is not computed, this array must have at least one cell for the error return identification. On output, JC(1)
is N
if control is returned normally.
V
is a complex variable. On input REAL(V)
is the option indicator, set as follows:
- invert matrix
- compute determinant
- do 1. and 2.
- solve system of equations
- do 1. and 4.
- do 2. and 4.
- do 1., 2. and 4.
Notes on usage of row dimension arguments N
and NR
:
The arguments N
and NR
refer to the row dimensions of the A
matrix.
N
gives the number of rows operated on by the subroutine, while NR
refers to the total number of rows in the matrix as dimensioned by the
calling program. NR
is used only in the dimension statement of the
subroutine. Through proper use of these parameters, the user may specify that only a submatrix, instead of the entire matrix, be operated on by the subroutine.
In your application (pipe flow), look at how matrix A
and vector V
are populated before the call to GJR
and how they are used after the call.
You may be able to replace the call to GJR
with a call to LAPACK's SGESV
or DGESV
without much difficulty.
Aside: The Fortran community really needs a drop-in 'Rosetta library' that wraps LAPACK, etc. for replacing legacy/proprietary IBM, UNIVAC, and Numerical Recipes math functions. The perfect case would be that maintainers would replace legacy functions with de facto standard math functions but in the real world, many of these older programs are un(der)maintained and there simply isn't the will (or, as in this case, the ability) to update them.
Update 2:
I started work on a compatibility library for the Sperry MATH-PACK and STAT-PACK routines as well as a few other legacy libraries, posted at https://bitbucket.org/apthorpe/alfc
Further, I located my copy of Jeppson's Analysis of Flow in Pipe Networks which is a slightly more legible version of the PDF of Steady Flow Analysis of Pipe Networks: An Instructional Manual and modernized the codes listed in the text. I have posted those at https://bitbucket.org/apthorpe/jeppson_pipeflow
Note that I found a number of errors in both the code listings and in the example problems given for many of the codes. If you're trying to learn how to write a pipe flow solver based on Jeppson's paper or text, I'd strongly suggest reviewing my updated codes and test cases because they will save you hours of effort trying to understand why the code doesn't work and why you can't replicate the example cases. This took a fair amount of forensic computing to sort out.
Update 3:
The source to CGJR
and DGJR
can be found in http://www.dtic.mil/dtic/tr/fulltext/u2/a110089.pdf. DGJR
is the closest to what you want, though it references more routines that aren't available (proprietary UNIVAC error-handling routines). It should be easy to convert `DGJR' to single precision and skip the proprietary calls. Otherwise, use the compatibility library mentioned above.
$98
is strange.. <non standard> alternate return?