老式 FORTRAN 中的 Univac Math pack 子例程(77 之前)

Univac Math pack subroutines in old-school FORTRAN (pre-77)

我一直在看一篇工程论文 here,它描述了一个用于求解管道流方程的旧 FORTRAN 代码(它的日期是 1974 年,在 FORTRAN 被标准化为 Fortran 77 之前)。在本文档的第 42 页,旧代码调用以下子例程:

C SYSTEM SUBROUTINE FROM UNIVAC MATH-PACK TO
C SOLVE LINEAR SYSTEM OF EQ. 
CALL GJR(A,51,50,NP,NPP,,JC,V)

有点扯远了,有老手或者老代码爱好者还记得这个系统子程序和它的输入参数吗?我无法找到有关它的任何信息。

如果我可以将旧代码改编成我当前的应用程序,我可能会用 C++ 或 VBA 重写它,并且会在这些语言中寻找等效的函数。

如果我发现更详细的内容,我会添加到这个答案中,但我有一个地方可以开始寻找 GJR.

的参数

此函数是 Sperry UNIVAC MATH-PACK 库的一部分 - 库中函数的完整列表可在 http://www.dtic.mil/dtic/tr/fulltext/u2/a170611.pdf 中找到 GJR 描述为 "determinant; inverse; solution of simultaneous equations"。有点帮助。

更好的描述来自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.

这有点用处,但我们真正想要的是 "MATH-PACK, Programmer Reference",来自 Sperry-UNIVAC (Unisys) 的 UP-7542 Rev. 1 我找到了很多对这个文档的引用,但没有全文文档本身的 PDF。

我会查看函数调用中的参数、它们的设置方式以及结果的使用方式,然后在 LAPACK 或 BLAS 中寻找等效例程。参见 http://www.netlib.org/lapack/

我有几本关于管道网络的书,包括 Jeppson 的 "Analysis of Flow in Pipe Networks"(与 USU 托管的原始 PDF 中的同一作者)https://books.google.com/books/about/Analysis_of_flow_in_pipe_networks.html?id=peZSAAAAMAAJ - 我会看看是否可以挖掘它。这本书可能有一个比专有的 Sperry-UNIVAC 库更便携的矩阵求解器。

更新:

从页。 http://ngds.egi.utah.edu/files/GL04099/GL04099_1.pdf 中的 41 我找到了 CGJR 函数的文档,它是同一库中 GJR 的复杂版本。参数的唯一区别可能是变量类型(COMPLEX 而不是 REAL):


CGJR 是一个子程序,它通过使用带有列旋转的高斯-若尔当消元技术来求解联立方程、计算行列式、求逆矩阵或执行这三种操作的任意组合。

使用CGJR的步骤如下:

调用语句:CALL CGJR(A,NC,NR,N,MC,$K,JC,V)

哪里

  • A 是要求逆或行列式的矩阵。如果求解联立方程组,则矩阵的最后 MC-N 列是待求解方程组的常数向量。在输出时,如果计算了逆,它存储在 A 的前 N ​​列中。如果求解联立方程,则最后 MC-N 列包含解向量。 A 是复数数组。

  • NC是一个整数,表示数组A.

  • 的最大列数
  • NR是一个整数,表示数组A.
  • 的最大行数
  • N是一个整数,表示要操作的数组A的行数。
  • MC为数组的列数A,表示求解联立方程时的系数矩阵;否则它是一个虚拟变量。
  • K 是调用程序中的语句编号,如果检测到溢出或奇异性,则将控制 return 编辑到该调用程序中。 1) 如果检测到溢出,JC(1) 设置为归约的最后正确完成的行的负值,然后控制 returned 到调用程序中的语句编号 K。 2) 如果检测到奇点,JC(1) 设置为最后正确完成的行的编号,如果要计算行列式,则 V 设置为 (0.,0.)。然后控制被 returned 到调用程序中的语句编号 K
  • JCN 元素的一维置换数组,如果正在计算逆运算,则用于置换 A 的行和列。如果逆运算是未计算,此数组必须至少有一个单元格用于错误 return 标识。在输出中,如果控件 return 正常编辑,JC(1)N
  • V 是一个复杂的变量。输入REAL(V)是option indicator,设置如下:
    1. 逆矩阵
    2. 计算行列式
    3. 做 1. 和 2.
    4. 求解方程组
    5. 做 1. 和 4.
    6. 做 2. 和 4.
    7. 执行 1.、2. 和 4.

行维度参数NNR的用法说明:

参数 NNR 指的是 A 矩阵的行维度。 N 给出子程序操作的行数,而 NR 指的是矩阵中的总行数,由 调用程序。 NR 仅在维度语句中使用 子程序。通过正确使用这些参数,用户可以指定子程序仅对子矩阵进行操作,而不是对整个矩阵进行操作。


在您的应用程序(管道流)中,查看在调用 GJR 之前如何填充矩阵 A 和向量 V 以及在调用之后如何使用它们。

您可以毫不费力地将对 GJR 的调用替换为对 LAPACK 的 SGESVDGESV 的调用。

旁白:Fortran 社区确实需要一个包含 LAPACK、 的插件 'Rosetta library' 来替换 legacy/proprietary IBM、UNIVAC 和 数值食谱 数学函数。完美的情况是,维护者会用 de facto 标准数学函数替换遗留函数,但在现实世界中,许多这些较旧的程序都没有得到维护,根本就没有更新它们的意愿(或在本例中为能力)。

更新 2:

我开始为 Sperry MATH-PACK 和 STAT-PACK 例程以及其他一些遗留库开发兼容性库,发布在 https://bitbucket.org/apthorpe/alfc

此外,我找到了 Jeppson 的 管网流动分析 的副本,它是 管道稳态流动分析的 PDF 的更易读版本Networks: An Instructional Manual 并对文中列出的代码进行了现代化改造。我已将这些张贴在 https://bitbucket.org/apthorpe/jeppson_pipeflow

请注意,我在代码清单和针对许多代码给出的示例问题中发现了一些错误。如果您正在尝试学习如何根据 Jeppson 的论文或文本编写管道流量求解器,我强烈建议您查看我更新的代码和测试用例,因为它们会节省您数小时的精力来理解为什么代码没有工作以及为什么你不能复制示例案例。这需要大量的取证计算才能解决。

更新 3:

CGJRDGJR 的来源可以在 http://www.dtic.mil/dtic/tr/fulltext/u2/a110089.pdf 中找到。 DGJR 是最接近您想要的,尽管它引用了更多不可用的例程(专有的 UNIVAC 错误处理例程)。将“DGJR”转换为单精度并跳过专有调用应该很容易。否则,使用上面提到的兼容性库。