Univac Math Pack中GJR子routine参数查询及C++/VBA等效实现问询

I’ve got some hands-on experience with vintage mainframe FORTRAN and UNIVAC’s math libraries, so even though official documentation for GJR is extremely hard to track down these days, I can give you educated guesses about its parameters based on 1970s UNIVAC Math-Pack conventions and common linear solver designs. I’ll also cover modern alternatives to help you port this logic to C++ or VBA.

Parameter Breakdown (Inferred from 1970s Linear Solver Norms) #

First, let’s look at your original subroutine call:

CALL GJR(A,51,50,NP,NPP,$98,JC,V)

Here’s what each parameter likely does, aligned with how UNIVAC structured its math subroutines back then:

• A: 2D array (column-major, as is classic FORTRAN) holding the coefficient matrix. It’s almost certainly overwritten during computation—either with an LU decomposition or the matrix inverse, depending on the exact implementation. The 51 value here refers to the declared leading dimension (you’d have declared A(51,51) to leave a buffer for the 50x50 working matrix).
• 51: The leading dimension of array A. Early FORTRAN required this parameter to correctly index 2D arrays, as it needed to know how many elements were in each column to jump to the next row in memory.
• 50: The actual order of the linear system—meaning you’re solving a 50x50 matrix equation.
• NP: Likely the number of right-hand side vectors. If you’re solving multiple systems with the same coefficient matrix (e.g., Ax₁=V₁, Ax₂=V₂), this tells the subroutine how many vectors are stored in V.
• NPP: The leading dimension of the V array, analogous to the 51 for A—used to index the right-hand side/solution vectors correctly in memory.
• $98: An error jump label specific to UNIVAC FORTRAN. If the subroutine hits a singular matrix or other fatal error, it branches to this label in your code to handle the issue (like printing an error message or terminating gracefully).
• JC: Integer array storing column pivot indices. This tracks column swaps made during pivoting (to avoid division by zero and reduce numerical error)—it should be at least length 50 to match your 50x50 matrix.
• V: 2D array (or 1D if NP=1) holding the right-hand side vectors. After the subroutine runs, this array will be overwritten with the solution vectors for each system.

Almost certainly, GJR implements Gauss-Jordan elimination with column pivoting—a standard method in the 1970s for solving linear systems or computing matrix inverses.

Modern Equivalents for C++ & VBA #

C++ #

You’ve got solid options depending on your performance and simplicity needs:

• Eigen Library (recommended for ease of use and readability):

  #include <Eigen/Dense>
  
  // Set up your matrices (50x50 coefficient, 50xNP right-hand sides)
  Eigen::MatrixXd A(50, 50);
  Eigen::MatrixXd V(50, NP);
  
  // Populate A and V with your pipeline flow data...
  
  // Solve Ax = V using full-pivot LU decomposition (matches GJR's pivoting logic)
  Eigen::MatrixXd solution = A.fullPivLu().solve(V);
  

• BLAS/LAPACK (for low-level, high-performance code):
  Use DGESV (for solving linear systems with partial-pivot LU decomposition) or DGEJ (Gauss-Jordan elimination, though implementation names can vary across vendors).

VBA #

For Excel VBA, you can leverage built-in functions or roll your own custom routine:

• Built-in Excel Functions (simple for small-to-medium matrices like 50x50):

  Function SolveLinearSystem(A As Variant, V As Variant) As Variant
      ' Compute inverse of A, then multiply by V to get the solution
      Dim inverseA As Variant
      inverseA = WorksheetFunction.MInverse(A)
      SolveLinearSystem = WorksheetFunction.MMult(inverseA, V)
  End Function
  

  Note: If you need to handle singular matrices or want full control over pivoting logic, you’d need to implement a custom Gauss-Jordan or Gaussian elimination routine with column pivoting.

内容的提问来源于stack exchange，提问作者Petrichor