18
$\begingroup$

I am working on a project that requires repeated calls to HankelH1[0, r] for $r$ spanning the full real axis. When I use the mathematica routine, it can be as much as 100x slower than the corresponding scipy routine (scipy.special.hankel1(0, r)), depending upon the value of the argument. Using ExternalEvaluate and ExternalFunction to directly call the scipy routine, it is even worse.

Here is the code:

Timing[
  Mean[Table[HankelH1[0, RandomReal[]*10.], {j, 1, 100000}]]
][[1]]

This takes four seconds to run!

Note that reducing the argument can make it faster or slower, as HankelH1 switches from Taylor series to asymptotic expansions to g-d knows what in the middle, but I would love a 100x faster solution -- ideally one that is platform independent (though if there is an easy way to use MathLink to call a Fortran routine or something, I would love to be introduced to that).

$\endgroup$
8
  • 1
    $\begingroup$ ParallelTable speeds it up a bit. $\endgroup$
    – Bill Watts
    Commented May 15, 2023 at 5:01
  • $\begingroup$ It should not be to hard to write a LibraryLink wrapper for some FORTRAN, C or C++ library. I could show you how to do that later the day. The only difficult part is to decide which library to link. If you want it platform independent, then the Intel Math Kernel Library is out. I see that boost implements Hankel functions: boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/… . I am a bit confused that they are called "Cyclic Hankel Functions" there. Do you think that these are the correct functions? $\endgroup$ Commented May 15, 2023 at 5:47
  • $\begingroup$ Hi Henrik- thanks for your reply -- i checked and those are indeed the correct functions! Any help you could provide on this front would be greatly appreciated. I am running OSX. $\endgroup$
    – JonBean
    Commented May 15, 2023 at 5:53
  • $\begingroup$ "I am running OSX." That's great because I do so too! (That means I can give you quite good guidance with regard to the linking.) Do you run an Intel or an Apple Silicon CPU? $\endgroup$ Commented May 15, 2023 at 6:06
  • 1
    $\begingroup$ Oh yes, "arch -arm64" is super important. I think I have coded it down already. Have to test it yet, but have to wait for some other job to finish first... $\endgroup$ Commented May 15, 2023 at 7:14

1 Answer 1

23
$\begingroup$

This is LibraryLink code with compiler and linker options specifically for macos on Apple Silicon.

These are the prerequisites to compile and run this successfully:

  1. XCode installed along with the Command Line Tools.

  2. Homebrew

  3. boost. Install it with arch -arm64 brew install boost from the command line.

  4. OpenMP. Install it with arch -arm64 brew install libomp from the command line. (You can skip this step if you use cHankelH1Thread from the bottom of this post).

Then you can copy, paste, and run the following code in Mathematica:

Needs["CCompilerDriver`"]

$OpenMPIncludeDirectory = "/opt/homebrew/opt/libomp/include";
$OpenMPLibraryDirectory = "/opt/homebrew/opt/libomp/lib";

$BoostIncludeDirectory = "/opt/homebrew/opt/boost/include";
$BoostLibraryDirectory = "/opt/homebrew/opt/boost/lib";

Quiet[LibraryFunctionUnload[cHankelH1OpenMP]];
ClearAll[cHankelH1OpenMP];

cHankelH1OpenMP::usage = "";

cHankelH1OpenMP = Module[{name, code, lib},
   name = "cHankelH1OpenMP";
   
   code = StringJoin["
#include \"WolframLibrary.h\"

#include <algorithm>

#include <boost/math/special_functions/hankel.hpp>

#include <omp.h>

EXTERN_C DLLEXPORT int ", name, 
     "(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor nu_          = MArgument_getMTensor(Args[0]);
    MTensor x_           = MArgument_getMTensor(Args[1]);

    mint    thread_count = MArgument_getInteger(Args[2]);
    
    const mint n = std::min(
        libData->MTensor_getDimensions(nu_)[0],
        libData->MTensor_getDimensions(x_)[0]
    );

    MTensor y_;

    (void)libData->MTensor_new( MType_Complex, 1, &n, &y_);

    const mreal * const nu = libData->MTensor_getRealData(nu_);
    const mreal * const x  = libData->MTensor_getRealData(x_);
         
    std::complex<mreal> * const y = reinterpret_cast<std::complex<mreal>*>(libData->MTensor_getComplexData(y_));
    

    #pragma omp parallel for num_threads( thread_count ) 
    for( mint i = 0; i < n ; ++i )
    {
        y[i] = boost::math::cyl_hankel_1( nu[i], x[i] );
    }

    MArgument_setMTensor(Res, y_);

    return LIBRARY_NO_ERROR;
}"
     ];
   
   lib = CreateLibrary[code, name,
     "TargetDirectory" -> $TemporaryDirectory,
     "Language" -> "C++",
     "CompileOptions" -> {"-Wall", "-Wextra", "-Wno-unused-parameter",
        "-Wno-deprecated-declarations", "-mmacosx-version-min=13.0", 
       "-std=c++17", "-Xpreprocessor -fopenmp", "-fno-math-errno", 
       "-mcpu=apple-m1 -mtune=native", "-ffast-math", "-Ofast", 
       "-flto", "-gline-tables-only", "-gcolumn-info"}, 
     "LinkerOptions" -> {"-lm", "-ldl", "-lomp"},
     "IncludeDirectories" -> {
       $OpenMPIncludeDirectory,
       $BoostIncludeDirectory
       },
     "LibraryDirectories" -> {
       $OpenMPLibraryDirectory,
       $BoostLibraryDirectory
       },
     "ShellOutputFunction" -> Print
     (*,"ShellOutputFunction"\[Rule]Print,*)
     
     ];
   
   LibraryFunctionLoad[lib, name,
    {
     {Real, 1, "Constant"}, (*First argument is a vector of nu values.*)
     {Real, 1, "Constant"}, (*Second argument is a vector of x values.*)
     Integer (*Third argument is the number of threads to use.*)
     },
    {Complex, 1}
    ]
   ];

And here a usage example with timings from my M1 Max (8 Cores) and a check for correctness:

n = 1000000;
\[Nu] = N@RandomInteger[{0, 10}, n];
x = RandomReal[{0, 100}, n];
threadCount = "ParallelThreadNumber" /. ("ParallelOptions" /. SystemOptions[]);

result1 = HankelH1[\[Nu], x]; // RepeatedTiming // First
result2 = cHankelH1OpenMP[\[Nu], x, threadCount]; // RepeatedTiming // First

Max[Abs[1 - result2/result1]]

27.4367

0.0116933

3.10116*10^-12

So we get a 2300-fold speedup.

Btw., the developers of the special functions in boost point out that their first and foremost goal was accuracy, not performance. So it is very likely that some other library might provide even better performance. I showed you how to use the LibraryLink interface. Of course you can link any library you like!

Edit

Because I frequently run into issues with OpenMP on various systems, I designed a hopefully more portable version by replacing OpenMP with std::thread. While the compilation options are still for Apple Silicon, it should be easier now to adapt them to other systems. The library boost is used nonetheless.

Needs["CCompilerDriver`"]

$BoostIncludeDirectory = "/opt/homebrew/opt/boost/include";

Quiet[LibraryFunctionUnload[cHankelH1Thread]];
ClearAll[cHankelH1Thread];

cHankelH1Thread::usage = "";

cHankelH1Thread = Module[{name, code, lib}, name = "cHankelH1Thread";
   code = StringJoin["
#include \"WolframLibrary.h\"

#include <algorithm>

#include <thread>

#include <boost/math/special_functions/hankel.hpp>

// Computes k-th job pointer for job_count equally sized jobs distributed on thread_count threads.
inline mint JobPointer( const mint job_count, const mint thread_count, const mint k )
{
    return (job_count / thread_count) * k + (job_count % thread_count) * k/thread_count;
}

EXTERN_C DLLEXPORT int ", name, 
     "(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor nu_          = MArgument_getMTensor(Args[0]);
    MTensor x_           = MArgument_getMTensor(Args[1]);

    mint    thread_count = MArgument_getInteger(Args[2]);
    
    const mint n = std::min(
        libData->MTensor_getDimensions(nu_)[0],
        libData->MTensor_getDimensions(x_)[0]
    );

    MTensor y_;

    (void)libData->MTensor_new( MType_Complex, 1, &n, &y_);

    const mreal * const nu = libData->MTensor_getRealData(nu_);
    const mreal * const x  = libData->MTensor_getRealData(x_);
         
    std::complex<mreal> * const y = reinterpret_cast<std::complex<mreal>*>(libData->MTensor_getComplexData(y_));
    
    auto F = [y,nu,x]( const mint i_begin, const mint i_end ) -> void
    { 
        for( mint i = i_begin; i < i_end; ++i )
        {
            y[i] = boost::math::cyl_hankel_1( nu[i], x[i] );
        }
    };

    std::vector<std::thread> threads;

    threads.reserve(thread_count);

    for( mint t = 0; t < thread_count; ++t )
    {
        threads.emplace_back (
            F, 
            JobPointer(n, thread_count, t   ),  
            JobPointer(n, thread_count, t+1 )
        );
    }

    for( auto & thread : threads )
    {
        thread.join();
    }

    MArgument_setMTensor(Res, y_);

    return LIBRARY_NO_ERROR;
}"];
   lib = CreateLibrary[code, name,
     "TargetDirectory" -> $TemporaryDirectory,
     "Language" -> "C++",
     "CompileOptions" -> {
       "-Wall", "-Wextra", 
       "-Wno-unused-parameter", "-Wno-deprecated-declarations",
       "-std=c++14", 
       "-mcpu=apple-m1 -mtune=native", 
       "-Ofast", 
       "-gline-tables-only", "-gcolumn-info"
       },
     "IncludeDirectories" -> {$BoostIncludeDirectory},
     "ShellOutputFunction" -> Print
    ];
   LibraryFunctionLoad[lib, name,
    {
     {Real, 1, "Constant"},(*First argument is a vector of nu values.*)
     {Real, 1, "Constant"},(*Second argument is a vector of x values.*)
     Integer (*Third argument is the number of threads to use.*)
     },
    {Complex, 1}
    ]];

Test:

n = 1000000;
\[Nu] = N@RandomInteger[{0, 10}, n];
x = RandomReal[{0, 100}, n];
threadCount = "ParallelThreadNumber" /. ("ParallelOptions" /. SystemOptions[]);

result1 = HankelH1[\[Nu], x]; // RepeatedTiming // First

result2 = cHankelH1OpenMP[\[Nu], x, threadCount]; // RepeatedTiming // First
result3 = cHankelH1Thread[\[Nu], x, threadCount]; // RepeatedTiming // First
Max[Abs[1 - result1/result2]]
Max[Abs[1 - result1/result3]]

27.5302

0.011383

0.0124506

3.5144*10^-12

3.5144*10^-12

As we can see here, cHankelH1Thread tends to be about 10% slower than cHankelH1OpenMP for n = 1000000. For n = 100000, it is even 40% slower. But for greater n, the relative difference seems to converge to 0.

$\endgroup$
3
  • $\begingroup$ woah just circled back and so this! Thank you so much! Amazing answer!!! $\endgroup$
    – JonBean
    Commented May 16, 2023 at 4:45
  • 1
    $\begingroup$ i just tried it and indeed i see a 1200x speedup! this is so terrific! I cannot thank you enough Henrik! $\endgroup$
    – JonBean
    Commented May 16, 2023 at 5:04
  • $\begingroup$ You're welcome! =) $\endgroup$ Commented May 16, 2023 at 5:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.