如何从Julia调用Numerical Recipes的svdcmp？遇Fortran size函数异常

Got it, let's figure out why your Fortran svdcmp subroutine is returning wonky size() results when called from Julia. I've run into this exact kind of interop issue before, so here's what's likely going on and how to fix it:

Common Culprits & Fixes #

1. Array Layout Mismatch (Row-Major vs Column-Major) #

First off: Fortran uses column-major (columns stored first) memory layout, while Julia defaults to row-major (rows stored first). If you're passing a Julia array straight to Fortran without adjusting, the dimensions Fortran sees will be transposed—so size(A,1) will give you Julia's column count instead of row count, and vice versa. That's a super common gotcha.

How to fix it:

• Either use the FortranArrays.jl package to create native column-major arrays in Julia, or explicitly transpose your array before passing it. For example:

  # Convert Julia's row-major array to Fortran's column-major
  fortran_A = permutedims(your_julia_array)
  # Now pass fortran_A to your svdcmp library
  

• Alternatively, when setting up your ccall in Julia, make sure you pass the dimensions in the reverse order that you'd use for Julia (since Fortran sees the transposed array).

2. Missing Explicit Dimension Arguments #

Older Numerical Recipes code (especially the Fortran 77/90 versions) often relies on the caller to pass explicit m and n arguments (the rows and columns of the matrix) instead of letting Fortran infer them with size(). If you're not passing these correctly from Julia, Fortran's size() will read garbage values from memory.

How to fix it:

• Check the original svdcmp subroutine signature—you'll almost certainly see it expects m and n as input arguments. Make sure you're passing the correct values from Julia, adjusted for the column-major layout if needed.
• Example ccall setup:

  # Get Julia's array dimensions (rows, columns)
  julia_m, julia_n = size(your_julia_array)
  # Pass transposed dimensions if you're not permuting the array first
  fortran_m = julia_n
  fortran_n = julia_m
  # Call the Fortran library with explicit dimensions
  ccall((:svdcmp, "libsvdcmp.so"), 
        Cvoid, 
        (Ref{Int}, Ref{Int}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
        fortran_m, fortran_n, your_julia_array, s, u, v)
  

3. Assumed-Shape Array Issues in Fortran #

If your svdcmp uses assumed-shape arrays (like real :: A(:,:)), Fortran needs hidden metadata about the array's dimensions to make size() work. If you compiled the shared library with outdated flags, or Julia isn't passing the array in a way that preserves this metadata, size() will fail.

How to fix it:

• When compiling your Fortran code to a shared library, use flags that support modern Fortran features. For gfortran, that looks like:

  gfortran -shared -fPIC -O2 svdcmp.f90 -o libsvdcmp.so
  

• If the original code is fixed-form Fortran (super common in old Numerical Recipes), consider updating it to use explicit shape declarations or add an interface block to clarify dimensions for the compiler.

4. Wrong Argument Types in Julia's ccall #

If you're not using the right pointer/reference types in your ccall, Fortran will misinterpret the data it receives—including array metadata. For example, passing a raw integer instead of a Ref{Int} for dimensions can cause Fortran to read incorrect values.

How to fix it:

• Match every argument in the Fortran subroutine to the correct Julia type:
  • For integer scalars (like m and n), use Ref{Int} (since Fortran expects pointers for intent(in) scalars in most cases).
  • For arrays, use Ptr{Float64} (or the appropriate type, e.g., Ptr{Float32} if your Fortran code uses single-precision).
  • Make sure output arrays are pre-allocated in Julia before passing them to Fortran.

Quick Debug Hack #

To confirm what's happening, temporarily modify your Fortran svdcmp subroutine to print the results of size(A,1) and size(A,2) as soon as it's called. Compare those numbers to the dimensions you're passing from Julia—this will instantly tell you if there's a layout mismatch or argument-passing error.

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