DEV Community: chrisarg

Coding a Perl interface to a foreign function with AI coding assistance

chrisarg — Thu, 04 Sep 2025 09:10:39 +0000

My takes as an amateur coder after a "vibecoding" session with #AI in #GithubCopilot:
👉 Agentic bots are unhelpful
🤏“Ask” mode could generate useful code, riddled with mistakes.
👍As auto-complete for documentation helping me overcome writer's block

Full prompt and analysis of Perl generated code 👇

https://chrisarg.github.io/Killing-It-with-PERL/2025/09/02/AI-assisted-coding-of-FFI.html

Taking VelociPerl for a ride.

chrisarg — Thu, 04 Sep 2025 09:08:13 +0000

VelociPerl is a closed source fork of Perl that claims performance gains of 45% over the stock ("your dad's Perl" in their parlance) based on some public benchmarks.
I will not go into how they achieved this performance boost, or why they released it as closed source, or even "but why the heck did you release it as closed source?", as you can follow the Reddit discussion.
However, even a modest speed gain may be useful in some applications, so I decided to dive in a a little bit deeper.

Some of the benchmarks are numerical e.g. a linear system solver, generating random numbers but others are more relevant to garden variety Perl tasks, e.g. generating objects as blessed hashes.
These tasks may appear in the context of some applications, so it is nice to know there are some benefits, but why not come up with a composite task and benchmark it? My usage of Perl involves random number generation, operations on tasks, creation of objects, string concatenation and function calling, so I figured there should be a way to combine all of them together and take VelociPerl for a ride.

Consider the following code that executes two benchmarks:

Generating a Hash (in which the keys are random strings and the values are obtained as a blessed reference to an anonymous scalar)
Accessing the said hash

use v5.38;
use Benchmark qw(timethis);    # for benchmarking

my $Length_of_accesstest = 1_000_000;
my $Length_of_gentest=100_000;
my $class = 'HashingAround';    
sub generate_random_string(@alphabet) {
    my $len = scalar @alphabet;
    my $length = int( rand(100) ) + 1;    # random length between 1 and 100
    return join '', @alphabet[ map { rand  $len } 1 .. $length ];
}

my %hash_of_seqs;

for ( 1 .. $Length_of_accesstest ) {
    my $key = generate_random_string(('A' .. 'Z'));
    my $val = bless \do { my $anon_scalar = $key }, $class;
    $hash_of_seqs{$key} = $val;
}

say "=" x 80;
say "\tTiming hash access using timethis";
timethis(
    20,
    sub {
        while ( my ( $seq_id, $seq ) = each %hash_of_seqs ) {
            ## worthless access to seq to force the sequence to be copied in memory
            my $sequence = $seq;
        }
    }
);

say "=" x 80;
say "\tTiming hash generation using timethis";
timethis(
    10,
    sub {
        my %hash_of_seqs;
        for ( 1 .. $Length_of_gentest ) {
    my $key = generate_random_string(('A' .. 'Z'));
    my $val = bless \do { my $anon_scalar = $key }, $class;
            $hash_of_seqs{$key} = $val;
        }
    }
);

To switch the Perl interpreter one simply changes the shebang line. In my oldish dual Xeon E5-2697 v4, I obtained the following results

Stock Perl :

(base) chrisarg@chrisarg-HP-Z840-Workstation:~/software-dev/velociperlHash$ ./testHash_speed.pl 
================================================================================
        Timing hash access using timethis
timethis 20: 12 wallclock secs (11.61 usr +  0.00 sys = 11.61 CPU) @  1.72/s (n=20)
================================================================================
        Timing hash generation using timethis
timethis 10: 12 wallclock secs (12.58 usr +  0.01 sys = 12.59 CPU) @  0.79/s (n=10)

VelociPerl :

(base) chrisarg@chrisarg-HP-Z840-Workstation:~/software-dev/velociperlHash$ ./testHash_speed_vperl.pl 
================================================================================
        Timing hash access using timethis
timethis 20: 11 wallclock secs (11.04 usr +  0.00 sys = 11.04 CPU) @  1.81/s (n=20)
================================================================================
        Timing hash generation using timethis
timethis 10: 10 wallclock secs ( 9.80 usr +  0.01 sys =  9.81 CPU) @  1.02/s (n=10)

Walking the hash was not materially different between the two interpreters, but the more compute intense task that involved random numbers, string operation and blessing of objects was ~30% faster.
This gain may or may not justify using a closed source version of Perl to you. But it is worth noting that one may be able to tweak the compilation of the Perl source (which appears to be a major source of the claimed gains) to generate faster executing Perl code. Perhaps an approach that one can try with an open sourced fork?

Vibe coding a Perl interface to a C library - Part 1

chrisarg — Tue, 01 Jul 2025 03:10:07 +0000

Introduction

In this multipart series we will explore the benefits (and pitfalls) of vibe coding a Perl interface to an external (or foreign) library through a large language model.
Those of you who follow me on X/Twitter (as @ChristosArgyrop and @ArgyropChristos),
Bluesky , mast.hpc, mstdn.science,
mstdn.social know that I have been very critical of the hype behind AI and the hallucinations of both models and the
human AI influencers (informally known as botlickers in some corners of the web).

However, there are application areas of vibe coding with AI, e.g. semi-automating the task of creating API from one one language to another, in which the chatbots may actually deliver well and act as productivity boosters.
In this application area, we will be leveraging the AI tools as more or less enhanced auto-complete tools that can help a developer navigate less familiar, more technical and
possibly more boring aspects of the target language's 'guts'. If AI were to deliver in this area, then meaningless language wars can be averted, and wrong (at least in my opinion) reasons
to prefer one language, i.e. the availability of some exotic library, may be avoided.

For my foray in this area, I chose to interface to the library Bit that I wrote to support text fingerprinting for some research
applications. The library based on David Hanson's Bit_T library discussed in
Chapter 13 of "C Interfaces and Implementations"
has been extended to incorporate additional operations on bitsets (such as counts on unions/differences/intersections of sets) and fast population counts in both CPU and GPU
(with TPU implementations coming down the road). Hence, this is a test case that can illustrate the utility of Perl in using code that executes transparently in various hardware assets.
Similar to Hanson's original implementation (after all, my work is based on his!) the library interface is implemented to
an Abstract Data Type in C; a crucial aspect of the implementation is to manage memory and avoid leaks without looking (at all, or as little as possible!) under the hood.
For our experiments we will use Claude 3.7 Thinking Sonnet through the Github Copilot vscode interface. This is going to be a multipart series post that will be published throughout the summer.
Our focus will be on interactions between me and the bot, and in particular critiquing the responses it has given me for both high level (e.g. choice of approach) and technical aspects of the project.

The Prompt

For the prompt I provided the "bit.c", "bit.h", makefile and the readme file of the Bit github repository as context to Claude and then I issued the following:

Look at the bit.h file that defines two abstract data types Bit_T and Bit_T_DB
and their interfaces in C. I would like to wrap these interfaces in Perl.
The build process will  consist of using the Alien package in Perl to make bit available to Perl.
Explore possible options for the creation of the interface: 
A. Creation of the bit  C library :
1) compiling bit.c and bit.h as dynamic libraries
2) compiling bit.c nd bit.h into a static library
B. Interfacing with the bit C library from Perl:
1) using SWIG i files
2) using Platypus::FFI
3) using Inline::C
4) using XS

Since you are an expert with deep knowledge of the Perl ecosystem, consider pros and cons of
each approach. and then recommend the optimal Consider lightweight options for encapsulating
the Bit_T and Bit_T_DB abstract data types as Perl objects.
Provide wrapping recommendations targetting a beginner, intermediate and experienced Perl
programmer since those may differ.  When providing recommendations , consider the dependencies
that are required for each implementation and highlight them in your analysis.
Once you have provided the recommendation, provide the implementation of the wrapper for
these 3 programmer levels. To reduce the token output (and your risk of hallucinating),
provide a full implementation of the wrapper for the beginner level and only the
changes (if any) that are required for the intermediate and advanced programmer levels.

In my prompt specification, I clearly specified that

linkage to the foreign library could be either as a static or a dynamic object,
the desired way to package the dependency for use in Perl is an Alien package
all four possible ways to interface to a foreign library from Perl should be considered.
asked it to consider the problem of wrapping from the perspective of Perl developers with different levels of expertise.

Note that I gave the bot considerable freedom over the "low" level aspects of wrapping (type of library, choice of the cross-language interface),
but restricted it in terms of how I wanted the library served.
Furthermore, I guided it to consider alternative viewpoints and level of expertise, providing a unique opportunity to test Claude's
representation¹ of conversations about wrappers. Finally, I alerted the chatbot to consider the token-sequence length dependent risk of hallucinations and
tried to limit the length of the output of the alternative implementations.

Claude's Response

The output based on this prompt is shown below (inside the horizontal lines); the formatting was also part of the output

Perl Wrapper for Bit_T and Bit_T_DB C Library

Let me analyze the options for creating a Perl interface to the Bit_T and Bit_T_DB C library and provide implementation recommendations for programmers at different experience levels.

A. Creation of the Bit C Library

1. Dynamic Library Approach

I see the Makefile already creates a dynamic library (libbit.so) in the build directory with this rule:

$(TARGET): $(OBJ)
    $(CC) $(CFLAGS) -shared -o $@ $^ $(LDFLAGS)

Pros:

Runtime binding, reducing memory footprint
Updates to the library don't require recompiling Perl code
Already implemented in your build system
Better for distributing the library separately from the wrapper

Cons:

Deployment requires ensuring the library is in the system path
Potential version compatibility issues

2. Static Library Approach

To compile as a static library, you would modify the Makefile:

TARGET = $(BUILD_DIR)/libbit.a

$(TARGET): $(OBJ)
    ar rcs $@ $^

Pros:

Self-contained distribution
No runtime dependency resolution
Potentially faster function calls
No path issues during deployment

Cons:

Larger binary size
Need to recompile wrapper if library changes
Less memory efficient if multiple programs use the library

B. Perl Interfacing Options

1. SWIG

Pros:

Automated wrapper generation
Supports multiple languages (future portability)
Good for complex interfaces

Cons:

Learning curve for SWIG interface files
Additional build dependency
Complex for beginners

2. FFI::Platypus

Pros:

No compilation required at installation time
Modern, actively maintained
Cleaner code than XS
Suitable for your library structure

Cons:

Requires FFI::Platypus and related modules
Less efficient for very complex interfaces

3. Inline::C

Pros:

Simple integration of C code directly in Perl
No separate interface files
Good for prototyping

Cons:

C compiler required at runtime
Less suitable for large libraries
Performance overhead for initialization

4. XS

Pros:

Standard Perl way to interface with C
Excellent performance
Tight integration with Perl's memory management

Cons:

Steep learning curve
Verbose, low-level coding
More complex build process

Recommendation for Each Level

Beginner Level: FFI::Platypus with Dynamic Library

FFI::Platypus provides the easiest approach for beginners with minimal C knowledge. It works with your existing dynamic library build.

Intermediate Level: Inline::C with Dynamic Library

Inline::C offers more control for intermediate users who understand some C code but want to avoid complex build processes.

Advanced Level: XS with Static Library

XS provides the highest performance and greatest control, suitable for experienced Perl programmers.

Implementation for Beginner Level: FFI::Platypus

Let's create the Alien::Bit module first to handle library discovery, then implement Bit::Set and Bit::Set::DB modules using FFI::Platypus.

Step 1: Create Alien::Bit module

# /media/chrisarg/Software-Dev/Bit/perl-wrapper/lib/Alien/Bit.pm
package Alien::Bit;

use strict;
use warnings;
use parent 'Alien::Base';
use File::Spec;

sub dynamic_libs {
    my $class = shift;

    # Try to find the library in common locations
    my @potential_paths = (
        # From the build directory
        '/media/chrisarg/Software-Dev/Bit/bit/build/libbit.so',

        # System paths
        '/usr/local/lib/libbit.so',
        '/usr/lib/libbit.so',

        # Environment variable if set
        $ENV{BIT_LIB_PATH} ? File::Spec->catfile($ENV{BIT_LIB_PATH}, 'libbit.so') : ()
    );

    foreach my $path (@potential_paths) {
        return $path if -f $path && -r $path;
    }

    die "Could not find libbit.so library. Please set BIT_LIB_PATH environment variable to its location.\n";
}

sub include_dirs {
    return (
        '/media/chrisarg/Software-Dev/Bit/bit/include',
        '/usr/local/include',
        '/usr/include',
        $ENV{BIT_INCLUDE_PATH} || ()
    );
}

1;

Step 2: Create Bit::Set module for individual bitsets

# /media/chrisarg/Software-Dev/Bit/perl-wrapper/lib/Bit/Set.pm
package Bit::Set;

use strict;
use warnings;
use FFI::Platypus;
use FFI::Platypus::Memory qw(malloc free);
use FFI::Platypus::Buffer qw(scalar_to_buffer buffer_to_scalar);
use Alien::Bit;

# Create FFI::Platypus object
my $ffi = FFI::Platypus->new(api => 1);

# Add path to our dynamic library
$ffi->lib(Alien::Bit->dynamic_libs);

# Define opaque types for our bitset pointers
$ffi->type('opaque' => 'Bit_T');
$ffi->type('opaque*' => 'Bit_T_Ptr');

# Define our Bit_T functions
$ffi->attach(Bit_new => ['int'] => 'Bit_T' => sub {
    my ($xsub, $self, $length) = @_;
    die "Length must be a positive integer" unless defined $length && $length > 0;
    my $bit_set = $xsub->($length);
    die "Failed to create bit set" unless $bit_set;
    return bless { _handle => $bit_set }, $self;
});

$ffi->attach(Bit_free => ['Bit_T_Ptr'] => 'opaque' => sub {
    my ($xsub, $self) = @_;
    my $ptr = \$self->{_handle};
    return $xsub->($ptr);
});

$ffi->attach(Bit_length => ['Bit_T'] => 'int' => sub {
    my ($xsub, $self) = @_;
    return $xsub->($self->{_handle});
});

$ffi->attach(Bit_count => ['Bit_T'] => 'int' => sub {
    my ($xsub, $self) = @_;
    return $xsub->($self->{_handle});
});

$ffi->attach(Bit_buffer_size => ['int'] => 'int');

$ffi->attach(Bit_bset => ['Bit_T', 'int'] => 'void' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    $xsub->($self->{_handle}, $index);
});

$ffi->attach(Bit_bclear => ['Bit_T', 'int'] => 'void' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    $xsub->($self->{_handle}, $index);
});

$ffi->attach(Bit_get => ['Bit_T', 'int'] => 'int' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    return $xsub->($self->{_handle}, $index);
});

$ffi->attach(Bit_set => ['Bit_T', 'int', 'int'] => 'void' => sub {
    my ($xsub, $self, $lo, $hi) = @_;
    die "Low index must be non-negative" unless defined $lo && $lo >= 0;
    die "High index must be greater than or equal to low index" unless defined $hi && $hi >= $lo;
    $xsub->($self->{_handle}, $lo, $hi);
});

$ffi->attach(Bit_clear => ['Bit_T', 'int', 'int'] => 'void' => sub {
    my ($xsub, $self, $lo, $hi) = @_;
    die "Low index must be non-negative" unless defined $lo && $lo >= 0;
    die "High index must be greater than or equal to low index" unless defined $hi && $hi >= $lo;
    $xsub->($self->{_handle}, $lo, $hi);
});

# Comparison operations
$ffi->attach(Bit_eq => ['Bit_T', 'Bit_T'] => 'int' => sub {
    my ($xsub, $self, $other) = @_;
    die "Other bitset must be a Bit::Set object" unless ref $other eq ref $self;
    return $xsub->($self->{_handle}, $other->{_handle});
});

# Set operations
$ffi->attach(Bit_union => ['Bit_T', 'Bit_T'] => 'Bit_T' => sub {
    my ($xsub, $self, $other) = @_;
    die "Other bitset must be a Bit::Set object" unless ref $other eq ref $self;
    my $result_handle = $xsub->($self->{_handle}, $other->{_handle});
    return bless { _handle => $result_handle }, ref $self;
});

$ffi->attach(Bit_inter => ['Bit_T', 'Bit_T'] => 'Bit_T' => sub {
    my ($xsub, $self, $other) = @_;
    die "Other bitset must be a Bit::Set object" unless ref $other eq ref $self;
    my $result_handle = $xsub->($self->{_handle}, $other->{_handle});
    return bless { _handle => $result_handle }, ref $self;
});

# Count operations
$ffi->attach(Bit_inter_count => ['Bit_T', 'Bit_T'] => 'int' => sub {
    my ($xsub, $self, $other) = @_;
    die "Other bitset must be a Bit::Set object" unless ref $other eq ref $self;
    return $xsub->($self->{_handle}, $other->{_handle});
});

# Constructor and destructor
sub new {
    my ($class, $length) = @_;
    return $class->Bit_new($length);
}

sub DESTROY {
    my ($self) = @_;
    $self->Bit_free() if defined $self->{_handle};
}

# Convenient accessor methods
sub length {
    my ($self) = @_;
    return $self->Bit_length();
}

sub count {
    my ($self) = @_;
    return $self->Bit_count();
}

sub get {
    my ($self, $index) = @_;
    return $self->Bit_get($index);
}

sub set {
    my ($self, $index) = @_;
    $self->Bit_bset($index);
    return $self;
}

sub clear {
    my ($self, $index) = @_;
    $self->Bit_bclear($index);
    return $self;
}

sub set_range {
    my ($self, $lo, $hi) = @_;
    $self->Bit_set($lo, $hi);
    return $self;
}

sub clear_range {
    my ($self, $lo, $hi) = @_;
    $self->Bit_clear($lo, $hi);
    return $self;
}

sub equals {
    my ($self, $other) = @_;
    return $self->Bit_eq($other) ? 1 : 0;
}

sub union {
    my ($self, $other) = @_;
    return $self->Bit_union($other);
}

sub intersection {
    my ($self, $other) = @_;
    return $self->Bit_inter($other);
}

sub intersection_count {
    my ($self, $other) = @_;
    return $self->Bit_inter_count($other);
}

1;

Step 3: Create Bit::Set::DB module for packed containers

# /media/chrisarg/Software-Dev/Bit/perl-wrapper/lib/Bit/Set/DB.pm
package Bit::Set::DB;

use strict;
use warnings;
use FFI::Platypus;
use FFI::Platypus::Memory qw(malloc free);
use FFI::Platypus::Buffer qw(scalar_to_buffer buffer_to_scalar);
use Bit::Set;
use Alien::Bit;

# Create FFI::Platypus object
my $ffi = FFI::Platypus->new(api => 1);

# Add path to our dynamic library
$ffi->lib(Alien::Bit->dynamic_libs);

# Define opaque types for our bitset pointers
$ffi->type('opaque' => 'Bit_T');
$ffi->type('opaque' => 'Bit_T_DB');
$ffi->type('opaque*' => 'Bit_T_DB_Ptr');

# Define struct for SETOP_COUNT_OPTS
$ffi->type('record(SETOP_COUNT_OPTS)' => [
    num_cpu_threads => 'int',
    device_id => 'int',
    upd_1st_operand => 'bool',
    upd_2nd_operand => 'bool',
    release_1st_operand => 'bool',
    release_2nd_operand => 'bool',
    release_counts => 'bool'
]);

# Define Bit_T_DB functions
$ffi->attach(BitDB_new => ['int', 'int'] => 'Bit_T_DB' => sub {
    my ($xsub, $self, $length, $num_of_bitsets) = @_;
    die "Length must be a positive integer" unless defined $length && $length > 0;
    die "Number of bitsets must be a positive integer" unless defined $num_of_bitsets && $num_of_bitsets > 0;
    my $db = $xsub->($length, $num_of_bitsets);
    die "Failed to create bitset DB" unless $db;
    return bless { _handle => $db }, $self;
});

$ffi->attach(BitDB_free => ['Bit_T_DB_Ptr'] => 'opaque' => sub {
    my ($xsub, $self) = @_;
    my $ptr = \$self->{_handle};
    return $xsub->($ptr);
});

$ffi->attach(BitDB_length => ['Bit_T_DB'] => 'int' => sub {
    my ($xsub, $self) = @_;
    return $xsub->($self->{_handle});
});

$ffi->attach(BitDB_nelem => ['Bit_T_DB'] => 'int' => sub {
    my ($xsub, $self) = @_;
    return $xsub->($self->{_handle});
});

$ffi->attach(BitDB_count_at => ['Bit_T_DB', 'int'] => 'int' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    return $xsub->($self->{_handle}, $index);
});

$ffi->attach(BitDB_clear => ['Bit_T_DB'] => 'void' => sub {
    my ($xsub, $self) = @_;
    $xsub->($self->{_handle});
});

$ffi->attach(BitDB_clear_at => ['Bit_T_DB', 'int'] => 'void' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    $xsub->($self->{_handle}, $index);
});

$ffi->attach(BitDB_get_from => ['Bit_T_DB', 'int'] => 'Bit_T' => sub {
    my ($xsub, $self, $index) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    my $bit_handle = $xsub->($self->{_handle}, $index);
    return bless { _handle => $bit_handle }, 'Bit::Set';
});

$ffi->attach(BitDB_put_at => ['Bit_T_DB', 'int', 'Bit_T'] => 'void' => sub {
    my ($xsub, $self, $index, $bitset) = @_;
    die "Index must be non-negative" unless defined $index && $index >= 0;
    die "Bitset must be a Bit::Set object" unless ref $bitset eq 'Bit::Set';
    $xsub->($self->{_handle}, $index, $bitset->{_handle});
});

# CPU-specific intersection count function
$ffi->attach(BitDB_inter_count_cpu => ['Bit_T_DB', 'Bit_T_DB', 'SETOP_COUNT_OPTS'] => 'int*' => sub {
    my ($xsub, $self, $other, $opts) = @_;
    die "Other must be a Bit::Set::DB object" unless ref $other eq ref $self;
    $opts ||= { 
        num_cpu_threads => 1, 
        device_id => -1, 
        upd_1st_operand => 0, 
        upd_2nd_operand => 0,
        release_1st_operand => 0,
        release_2nd_operand => 0,
        release_counts => 0
    };
    return $xsub->($self->{_handle}, $other->{_handle}, $opts);
});

# Constructor and destructor
sub new {
    my ($class, $length, $num_of_bitsets) = @_;
    return $class->BitDB_new($length, $num_of_bitsets);
}

sub DESTROY {
    my ($self) = @_;
    $self->BitDB_free() if defined $self->{_handle};
}

# Convenient accessor methods
sub length {
    my ($self) = @_;
    return $self->BitDB_length();
}

sub num_of_bitsets {
    my ($self) = @_;
    return $self->BitDB_nelem();
}

sub count_at {
    my ($self, $index) = @_;
    return $self->BitDB_count_at($index);
}

sub clear {
    my ($self) = @_;
    $self->BitDB_clear();
    return $self;
}

sub clear_at {
    my ($self, $index) = @_;
    $self->BitDB_clear_at($index);
    return $self;
}

sub get_bitset {
    my ($self, $index) = @_;
    return $self->BitDB_get_from($index);
}

sub put_bitset {
    my ($self, $index, $bitset) = @_;
    $self->BitDB_put_at($index, $bitset);
    return $self;
}

sub intersection_count_cpu {
    my ($self, $other, $opts) = @_;
    return $self->BitDB_inter_count_cpu($other, $opts);
}

1;

Step 4: Example usage script

#!/usr/bin/perl

use strict;
use warnings;
use Bit::Set;
use Bit::Set::DB;

# Create a new bitset with 1024 bits
my $bitset1 = Bit::Set->new(1024);
my $bitset2 = Bit::Set->new(1024);

# Set some bits
$bitset1->set(42);
$bitset1->set(100);
$bitset2->set(42);
$bitset2->set(200);

# Check if a bit is set
print "Bit 42 in bitset1: ", $bitset1->get(42), "\n";
print "Bit 100 in bitset1: ", $bitset1->get(100), "\n";
print "Bit 200 in bitset1: ", $bitset1->get(200), "\n";

# Count the number of bits set
print "Number of bits set in bitset1: ", $bitset1->count(), "\n";
print "Number of bits set in bitset2: ", $bitset2->count(), "\n";

# Calculate intersection count
my $intersection_count = $bitset1->intersection_count($bitset2);
print "Intersection count: $intersection_count\n";

# Create a new bitset as the union of the two bitsets
my $union = $bitset1->union($bitset2);
print "Number of bits set in union: ", $union->count(), "\n";

# Create a BitDB with 3 bitsets of length 1024
my $db = Bit::Set::DB->new(1024, 3);

# Put our bitsets into the DB
$db->put_bitset(0, $bitset1);
$db->put_bitset(1, $bitset2);
$db->put_bitset(2, $union);

# Count bits in each bitset in the DB
print "Bits set in DB at index 0: ", $db->count_at(0), "\n";
print "Bits set in DB at index 1: ", $db->count_at(1), "\n";
print "Bits set in DB at index 2: ", $db->count_at(2), "\n";

Modifications for Intermediate Level (Inline::C)

For intermediate-level Perl programmers, here are the changes required to implement the wrapper using Inline::C:

# Bit::Set using Inline::C
package Bit::Set;

use strict;
use warnings;
use Inline C => Config =>
    BUILD_NOISY => 1,
    CLEAN_AFTER_BUILD => 0,
    LIBS => '-L/media/chrisarg/Software-Dev/Bit/bit/build -lbit',
    INC => '-I/media/chrisarg/Software-Dev/Bit/bit/include';

use Inline C => <<'END_C';
#include "bit.h"

typedef struct Bit_T* Bit_T;

// Wrapper functions
Bit_T create_bitset(int length) {
    return Bit_new(length);
}

void free_bitset(Bit_T set) {
    Bit_free(&set);
}

int get_length(Bit_T set) {
    return Bit_length(set);
}

int get_count(Bit_T set) {
    return Bit_count(set);
}

void set_bit(Bit_T set, int index) {
    Bit_bset(set, index);
}

void clear_bit(Bit_T set, int index) {
    Bit_bclear(set, index);
}

int get_bit(Bit_T set, int index) {
    return Bit_get(set, index);
}

// Add more wrapper functions as needed
END_C

# Perl OO interface
sub new {
    my ($class, $length) = @_;
    my $handle = create_bitset($length);
    return bless { _handle => $handle }, $class;
}

# Add more methods similar to the FFI version

Modifications for Advanced Level (XS)

For advanced Perl programmers, here are the key changes required to implement using XS:

# Bit::Set.xs
#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "bit.h"

typedef struct Bit_T* Bit_T;
typedef struct Bit_T_DB* Bit_T_DB;

MODULE = Bit::Set    PACKAGE = Bit::Set

PROTOTYPES: DISABLE

Bit_T
Bit_new(length)
    int length

void
Bit_free(set)
    Bit_T set
    CODE:
        Bit_free(&set);

int
Bit_length(set)
    Bit_T set

int
Bit_count(set)
    Bit_T set

void
Bit_bset(set, index)
    Bit_T set
    int index

# Add more functions as needed

And update the Makefile.PL:

use ExtUtils::MakeMaker;

WriteMakefile(
    NAME              => 'Bit::Set',
    VERSION_FROM      => 'lib/Bit/Set.pm',
    LIBS              => ['-L/media/chrisarg/Software-Dev/Bit/bit/build -lbit'],
    INC               => '-I/media/chrisarg/Software-Dev/Bit/bit/include',
    OBJECT            => '$(O_FILES)',
);

Conclusion

For most users, especially beginners, the FFI::Platypus approach offers the best balance of simplicity,
performance, and maintainability. It works well with your existing dynamic library build process and
provides a clean, object-oriented interface to the Bit_T and Bit_T_DB types. The implementation provided
covers all the core functionality while being extensible for additional features as needed.

My critique of Claude's choices

In this first post, we will not consider the implementation(s) which will be covered in separate posts, but will restrict ourselves to Claude's output with respect to the other questions. I give Claude bonus points for clearly identifying the niches for each potential choice:

Highighting that SWIG can support wrapping for other languages)
Proposing Foreign Function Interface for linking to the dynamic library for beginners.
A experienced programmer would opt for XS is not a surprise at all.

However:

the choice of Inline for the intermediate user is head-scratching: it seems that the chatbot closed on the intermediate level of programming experience in the prompt, and the selection of the approach was driven entirely by the fact that the user could (presumably) do more stuff in C.
SWIG was not considered as suitable (perhaps because few people in the training databases use SWIG) for implementing at any level. Without going into the specifics of the implementation though, I'd feel comfortable opting for FFI as an initial step for largely the reasons identified by Claude. We will have more things to say about the FFI implementation in the subsequent post in this series.

Note, I did not use the word understanding, as I do not think that LLMs can understant: they are merely noisy statistical pattern generators that can be tasked to create rough solutions for refining.
I alerted the bot to the (substantial) risk of hallucinations and decreased

[Boost]

chrisarg — Sat, 07 Jun 2025 19:55:10 +0000

LNATION for LNATION

Jun 6 '25

Learning Perl – Introduction

#perl #beginners #programming #tutorial

Comments 1

2 min read

Using Perl to Monitor Peak DRAM use in R

chrisarg — Sun, 19 Jan 2025 18:24:27 +0000

A story in two parts:

buff.ly/4hmFepy goes over the subtleties of monitoring DRAM use in #rstats
buff.ly/4hlMstF shows the #Perl solution and how one can make it play nice from within R

Code is released under the MIT license

Taking VelociPerl for a ride.

chrisarg — Sat, 21 Sep 2024 16:44:06 +0000

Consider the following code that executes two benchmarks:

Generating a Hash (in which the keys are random strings and the values are obtained as a blessed reference to an anonymous scalar)
Accessing the said hash

use v5.38;
use Benchmark qw(timethis);    # for benchmarking

my $Length_of_accesstest = 1_000_000;
my $Length_of_gentest=100_000;
my $class = 'HashingAround';    
sub generate_random_string(@alphabet) {
    my $len = scalar @alphabet;
    my $length = int( rand(100) ) + 1;    # random length between 1 and 100
    return join '', @alphabet[ map { rand  $len } 1 .. $length ];
}

my %hash_of_seqs;

for ( 1 .. $Length_of_accesstest ) {
    my $key = generate_random_string(('A' .. 'Z'));
    my $val = bless \do { my $anon_scalar = $key }, $class;
    $hash_of_seqs{$key} = $val;
}

say "=" x 80;
say "\tTiming hash access using timethis";
timethis(
    20,
    sub {
        while ( my ( $seq_id, $seq ) = each %hash_of_seqs ) {
            ## worthless access to seq to force the sequence to be copied in memory
            my $sequence = $seq;
        }
    }
);

say "=" x 80;
say "\tTiming hash generation using timethis";
timethis(
    10,
    sub {
        my %hash_of_seqs;
        for ( 1 .. $Length_of_gentest ) {
    my $key = generate_random_string(('A' .. 'Z'));
    my $val = bless \do { my $anon_scalar = $key }, $class;
            $hash_of_seqs{$key} = $val;
        }
    }
);

To switch the Perl interpreter one simply changes the shebang line. In my oldish dual Xeon E5-2697 v4, I obtained the following results

Stock Perl :

(base) chrisarg@chrisarg-HP-Z840-Workstation:~/software-dev/velociperlHash$ ./testHash_speed.pl 
================================================================================
        Timing hash access using timethis
timethis 20: 12 wallclock secs (11.61 usr +  0.00 sys = 11.61 CPU) @  1.72/s (n=20)
================================================================================
        Timing hash generation using timethis
timethis 10: 12 wallclock secs (12.58 usr +  0.01 sys = 12.59 CPU) @  0.79/s (n=10)

VelociPerl :

(base) chrisarg@chrisarg-HP-Z840-Workstation:~/software-dev/velociperlHash$ ./testHash_speed_vperl.pl 
================================================================================
        Timing hash access using timethis
timethis 20: 11 wallclock secs (11.04 usr +  0.00 sys = 11.04 CPU) @  1.81/s (n=20)
================================================================================
        Timing hash generation using timethis
timethis 10: 10 wallclock secs ( 9.80 usr +  0.01 sys =  9.81 CPU) @  1.02/s (n=10)

The Day Perl Stood Still: Unveiling A Hidden Power Over C

chrisarg — Fri, 16 Aug 2024 05:11:38 +0000

Sometimes the unexpected happens and must be shared with the world ... this one is such a case.

Recently, I've started experimenting with Perl for workflow management and high-level supervision of low level code for data science applications. A role I'd reserve for Perl in this context is that of lifecycle management of memory buffers, using the Perl application to "allocate" memory buffers and shuttle it between computing components written in C, Assembly, Fortran and the best hidden gem of the Perl world, the Perl Data Language.
There at least 3 ways that Perl can be used to allocate memory buffers:

Generate a list of bytes and use the pack function to convert them to a string.
Use the repetition operator (x) to generate a string of length equal to the size of the buffer one wants minus one (strings are terminated via null bytes in Perl, thus the null byte compensates for the mine one).
Access an external memory allocator library either through Inline or FFI::Platypus to allocate the buffer.

The following Perl code implements these three methods (pack , string and malloc in C) and allows one to experiment with different buffer sizes, initial values and precision of the results (by averaging over many iterations of the allocation routines)

#!/home/chrisarg/perl5/perlbrew/perls/current/bin/perl
use v5.38;
use Inline (
    C         => 'DATA',
    cc        => 'g++',
    ld        => 'g++',
    inc       => q{},      # replace q{} with anything else you need
    ccflagsex => q{},      # replace q{} with anything else you need
    lddlflags => join(
        q{ },
        $Config::Config{lddlflags},
        q{ },              # replace q{ } with anything else you need
    ),
    libs => join(
        q{ },
        $Config::Config{libs},
        q{ },              # replace q{ } with anything else you need
    ),
    myextlib => ''
);

use Benchmark qw(cmpthese);
use Getopt::Long;
my ($buffer_size, $init_value, $iterations);
GetOptions(
    'buffer_size=i' => \$buffer_size,
    'init_value=s'  => \$init_value,
    'iterations=i'  => \$iterations,
) or die "Usage: $0 --buffer_size <size> --init_value <value> --iterations <count>\n";
my $init_value_byte = ord($init_value);
my %code_snippets   = (
    'string' => sub {
        $init_value x ( $buffer_size - 1 );
    },
    'pack' => sub {
        pack "C*", ( ($init_value_byte) x $buffer_size );
    },
    'C' => sub {
        allocate_and_initialize_array( $buffer_size, $init_value_byte );
    },
);

cmpthese( $iterations, \%code_snippets );

__DATA__

__C__

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

SV* allocate_and_initialize_array(size_t length, short initial_value) {
    // Allocate memory for the array
    char* array = (char*)malloc(length * sizeof(char));
    char initial_value_byte = (char)initial_value;
    if (array == NULL) {
        fprintf(stderr, "Memory allocation failed\n");
        exit(1);
    }

    // Initialize each element with the initial_value
    memset(array, initial_value_byte, length);
    return newSVuv(PTR2UV(array));
}

calling the script as:

./time_mem_alloc.pl -buffer_size=1000000 -init_value=A -iterations=20000

yielded the surprising result :

          Rate   pack      C string
pack     322/s     --   -92%   -99%
C       4008/s  1144%     --   -92%
string 50000/s 15417%  1147%     --

with the Perl string method outperforming C by 10 fold.

Not believing the massive performance gain, and thinking I am dealing with a bug in Inline::C, I recoded the allocation in pure C (adding the usual embellishments for commandline processing/timing etc) :

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

char* allocate_and_initialize_array(size_t length, char initial_value) {
    // Allocate memory for the array
    char* array = (char*)malloc(length * sizeof(char));
    if (array == NULL) {
        fprintf(stderr, "Memory allocation failed\n");
        exit(1);
    }

    // Initialize each element with the initial_value
    memset(array, initial_value, length);

    return array;
}

double time_allocation_and_initialization(size_t length, char initial_value) {
    clock_t start, end;
    double cpu_time_used;

    start = clock();
    char* array = allocate_and_initialize_array(length, initial_value);
    end = clock();

    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
    /* This rudimentary loop prevents the compiler from optimizing out the 
     * allocation/initialization with the de-allocation
    */
    for(size_t i = 1; i < length; i++) {
       array[i]++;
       if(i % 100000 == 0) {
           printf("array[%zu] = %c\n", i, array[i]);
       }
    }
    free(array); // Free the allocated memory

    return cpu_time_used;
}

int main(int argc, char *argv[]) {
    if (argc != 3) {
        fprintf(stderr, "Usage: %s <length> <initial_value>\n", argv[0]);
        return 1;
    }

    size_t length = strtoull(argv[1], NULL, 10);
    char initial_value = argv[2][0];

    double time_taken = time_allocation_and_initialization(length, initial_value);
    printf("Time taken to allocate and initialize array: %f seconds\n", time_taken);
    printf("Initializes per second: %f\n", 1/time_taken);

    return 0;
}

/*
Compilation command:
gcc -O2 -o time_array_allocation time_array_allocation.c -std=c99

Example invocation:
./time_array_allocation 10000000 A
*/

Invoking the C program as the comment in the C code says to do,
I obtained the following result:

Time taken to allocate and initialize array: 0.000203 seconds
Initializes per second: 4926.108374

which is practically executes in the same order of magnitude as the equivalent allocation as the Inline::C malloc/C approach.
After researching the issue further, I discovered that the malloc I grew to admire trades speed in memory allocation for generality, and there is a plethora of faster memory allocators out there. It seems that Perl is using one such allocator for its strings, and kicks C's butt in this task of allocating buffers.

The Quest for Performance Part IV : May the SIMD Force be with you

chrisarg — Fri, 02 Aug 2024 05:01:28 +0000

At this point one may wonder how numba, the Python compiler around numpy Python code, delivers a performance premium over numpy. To do so, let's inspect timings individually for all trigonometric functions (and yes, the exponential and the logarithm are trigonometric functions if you recall your complex analysis lessons from high school!). But the test is not relevant only for those who want to do high school trigonometry: physics engines e.g. in games will use these function, and machine learning and statistical calculations heavily use log and exp. So getting the trigonometric functions right is one small, but important step, towards implementing a variety of applications. The table below shows the timings for 50M in place transformations:

Function	Library	Execution Time
Sqrt	numpy	1.02e-01 seconds
Log	numpy	2.82e-01 seconds
Exp	numpy	3.00e-01 seconds
Cos	numpy	3.55e-01 seconds
Sin	numpy	4.83e-01 seconds
Sin	numba	1.05e-01 seconds
Sqrt	numba	1.05e-01 seconds
Exp	numba	1.27e-01 seconds
Log	numba	1.47e-01 seconds
Cos	numba	1.82e-01 seconds

The table holds the first hold to the performance benefits: the square root, a function that has a dedicated SIMD instruction for vectorization takes exactly the same time to execute in numba and numpy, while all the other functions are speed up by 2-2.5 time, indicating either that the code auto-vectorizes using SIMD or auto-threads. A second clue is provided by examining the difference in results between the numba and numpy, using the ULP (Unit in the Last Place). ULP is a measure of accuracy in numerical calculations and can easily be computed for numpy arrays using the following Python function:

def compute_ulp_error(array1, array2):
## maxulp set up to a very high number to avoid throwing an exception in the code
    return np.testing.assert_array_max_ulp(array1, array2, maxulp=100000)

These numerical benchmarks indicate that the square root function utilizes pretty much equivalent code in numpy and numba, while for all the other trigonometric functions the mean, median, 99.9th percentile and maximum ULP value over all 50M numbers differ. This is a subtle hint that SIMD is at play: vectorization changes slightly the semantics of floating point code to make use of associative math, and floating point numerical operations are not associative.

Function	Mean ULP	Median ULP	99.9th ULP	Max ULP
Sqrt	0.00e+00	0.00e+00	0.00e+00	0.00e+00
Sin	1.56e-03	0.00e+00	1.00e+00	1.00e+00
Cos	1.43e-03	0.00e+00	1.00e+00	1.00e+00
Exp	5.47e-03	0.00e+00	1.00e+00	2.00e+00
Log	1.09e-02	0.00e+00	2.00e+00	3.00e+00

Finally, we can inspect the code of the numba generated functions for vectorized assembly instructions as detailed here, using the code below:

@njit(nogil=True, fastmath=False, cache=True)
def compute_sqrt_with_numba(array):
    np.sqrt(array, array)


@njit(nogil=True, fastmath=False, cache=True)
def compute_sin_with_numba(array):
    np.sin(array, array)


@njit(nogil=True, fastmath=False, cache=True)
def compute_cos_with_numba(array):
    np.cos(array, array)


@njit(nogil=True, fastmath=False, cache=True)
def compute_exp_with_numba(array):
    np.exp(array, array)


@njit(nogil=True, fastmath=False, cache=True)
def compute_log_with_numba(array):
    np.log(array, array)

## check for vectorization
## code lifted from https://tbetcke.github.io/hpc_lecture_notes/simd.html
def find_instr(func, keyword, sig, limit=5):
    count = 0
    for l in func.inspect_asm(func.signatures[sig]).split("\n"):
        if keyword in l:
            count += 1
            print(l)
            if count >= limit:
                break
    if count == 0:
        print("No instructions found")

# Compile the function to avoid the overhead of the first call
compute_sqrt_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_sin_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_exp_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_cos_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))
compute_exp_with_numba(np.array([np.random.rand() for _ in range(1, 6)]))

And this is how we probe for the presence f the vmovups instruction indicating that the YMM AVX2 registers are being used in calculations

print("sqrt")
find_instr(compute_sqrt_with_numba, keyword="vsqrtsd", sig=0)
print("\n\n")
print("sin")
find_instr(compute_sin_with_numba, keyword="vmovups", sig=0)

As we can see in the output below the square root function uses the x86 SIMD vectorized square root instruction vsqrtsd , and the sin (but also all trigonometric functions) use SIMD instructions to access memory.

sqrt
        vsqrtsd %xmm0, %xmm0, %xmm0
        vsqrtsd %xmm0, %xmm0, %xmm0


sin
        vmovups (%r12,%rsi,8), %ymm0
        vmovups 32(%r12,%rsi,8), %ymm8
        vmovups 64(%r12,%rsi,8), %ymm9
        vmovups 96(%r12,%rsi,8), %ymm10
        vmovups %ymm11, (%r12,%rsi,8)

Let's shift attention to Perl and C now (after all this is a Perl blog!). In Part I
we saw that PDL and C gave similar performance when evaluating the nested function cos(sin(sqrt(x))), and in Part II that the single-threaded Perl code was as fast as numba. But what about the individual trigonometric functions in PDL and C without the Inline module? The C code block that will be evaluated in this case is:

#pragma omp for simd
  for (int i = 0; i < array_size; i++) {
    double x = array[i];
    array[i] = foo(x);
  }

where foo is one of sqrt, sin, cos, log, exp. We will use the simd omp pragma in conjunction with the following compilation flags and the gcc compiler to see if we can get the code to pick up the hint and auto-vectorize the machine code generated using the SIMD instructions.

CC = gcc
CFLAGS = -O3 -ftree-vectorize  -march=native -mtune=native -Wall -std=gnu11 -fopenmp -fstrict-aliasing -fopt-info-vec-optimized -fopt-info-vec-missed
LDFLAGS = -fPIE -fopenmp 
LIBS =  -lm

During compilation, gcc informs us about all the wonderful missed opportunities to optimize the loop. The performance table below also demonstrates this; note that the standard C implementation and PDL are equivalent and equal in performance to numba. Perl through PDL can deliver performance in our data science world.

Function	Library	Execution Time
Sqrt	PDL	1.11e-01 seconds
Log	PDL	2.73e-01 seconds
Exp	PDL	3.10e-01 seconds
Cos	PDL	3.54e-01 seconds
Sin	PDL	4.75e-01 seconds
Sqrt	C	1.23e-01 seconds
Log	C	2.87e-01 seconds
Exp	C	3.19e-01 seconds
Cos	C	3.57e-01 seconds
Sin	C	4.96e-01 seconds

To get the compiler to use SIMD, we replace the flag -O3 by -Ofast and all these wonderful opportunities for performance are no longer missed, and the C code now delivers (but with the usual caveats that apply to the -Ofast flag).

Function	Library	Execution Time
Sqrt	C - Ofast	1.00e-01 seconds
Sin	C - Ofast	9.89e-02 seconds
Cos	C - Ofast	1.05e-01 seconds
Exp	C - Ofast	8.40e-02 seconds
Log	C - Ofast	1.04e-01 seconds

With these benchmarks, let's return to our initial Perl benchmarks and contrast the timings obtained with the non-SIMD aware invokation of the Inline C code:

use Inline (
    C    => 'DATA',
    build_noisy => 1,
    with => qw/Alien::OpenMP/,
    optimize => '-O3 -march=native -mtune=native',
    libs => '-lm'
);

and the SIMD aware one (in the code below, one has to incluce the vectorized version of the mathematics library for the code to compile):

use Inline (
    C    => 'DATA',
    build_noisy => 1,
    with => qw/Alien::OpenMP/,
    optimize => '-Ofast -march=native -mtune=native',
    libs => '-lmvec'
);

The non-vectorized version of the code yield the following table

Inplace  in  base Python took 11.9 seconds
Inplace  in PythonJoblib took 4.42 seconds
Inplace  in         Perl took 2.88 seconds
Inplace  in Perl/mapCseq took 1.60 seconds
Inplace  in    Perl/mapC took 1.50 seconds
C array  in            C took 1.42 seconds
Vector   in       Base R took 1.30 seconds
C array  in   Perl/C/seq took 1.17 seconds
Inplace  in     PDL - ST took 0.94 seconds
Inplace in  Python Numpy took 0.93 seconds
Inplace  in Python Numba took 0.49 seconds
Inplace  in   Perl/C/OMP took 0.24 seconds
C array  in   C with OMP took 0.22 seconds
C array  in    C/OMP/seq took 0.18 seconds
Inplace  in     PDL - MT took 0.16 seconds

while the vectorized ones this one:

Inplace  in  base Python took 11.9 seconds
Inplace  in PythonJoblib took 4.42 seconds
Inplace  in         Perl took 2.94 seconds
Inplace  in Perl/mapCseq took 1.59 seconds
Inplace  in    Perl/mapC took 1.48 seconds
Vector   in       Base R took 1.30 seconds
Inplace  in     PDL - ST took 0.96 seconds
Inplace in  Python Numpy took 0.93 seconds
Inplace  in Python Numba took 0.49 seconds
C array  in   Perl/C/seq took 0.30 seconds
C array  in            C took 0.26 seconds
Inplace  in   Perl/C/OMP took 0.24 seconds
C array  in   C with OMP took 0.23 seconds
C array  in    C/OMP/seq took 0.19 seconds
Inplace  in     PDL - MT took 0.17 seconds

To facilitate comparisons against the various flavors of Python and R, we inserted the results we presented previously in these 2 tables.

The take home points (some of which may be somewhat surprising are):

Performance of base Python was horrendous - nearly 4 times slower than base Perl. Even under parallel processing (Joblib) the 8 threaded Python code was slower than Perl
R performed the best out of the 3 dynamically typed languages considered : nearly 9 times faster than base Python and 2.3 times faster than base Perl.
Inplace modification of Perl arrays by accessing the array containers in C (Perl/mapC) narrowed the gap between Perl and R
Considering variability in timings, the three codes, base R, Perl/mapC and an equivalent inplace modification of a C array are roughly equivalent
The uni-threaded PDL code delivered the same performance as numpy
SIMD aware code, e.g. through numba or (for the case of native C code compiler pragmas, e.g. cntract C array in C timings between the vectorized and the non-vectorized versions of the table) delivered the best single thread performance.
OpenMP pragmas can really speed up operations (such as map) of Perl containers.
Adding threads (OMP code) or the multi-threaded versions of the PDL delivered the best performance.

These observations generate the following big picture questions:

Observations 1 and 2 make it quite surprising that base Python earned a reputation for a data science language. In fact, with these performance characteristics of a rather simple numerical code, one should approach the replacement of R (or Perl for those of us who use it) by base Python in data analysis codes.
Since hybrid implementations that involve C and Perl can deliver real performance benefits using SIMD aware code, even if a single thread is used, should we be upgrading our Inline::C codes to use SIMD friendly compiler flags? Or (for those who are are afraid of -Ofast, perhaps a carefully prepared mixology of intrinsics, and even assembly?
Should we be looking into upgrading various C/XS modules for Perl, so that they use OpenMP?
Why are not more people using the jewel of software that is PDL? The auto-threading property alone should make people think about using it for demanding, data intensive tasks.
Is it possible to create hybrid applications that rely on both PDL and OpenMP/Inline::C for different computations?
Should we look into special builds for PDL that leverage the SIMD properties and compiler pragmas to allow users to have their cake (SIMD vectorization) and eat it (autothreading) too?

The Quest for Performance Part III : C Force

chrisarg — Fri, 02 Aug 2024 04:58:10 +0000

In the two prior installments of this series, we considered the performance of floating operations in Perl,
Python and R in a toy example that computed the function cos(sin(sqrt(x))), where x was a very large array of 50M double precision floating numbers.
Hybrid implementations that delegated the arithmetic intensive part to C were among the most performant implementations. In this installment, we will digress slightly and look at the performance of a pure C code implementation of the toy example.
The C code will provide further insights about the importance of memory locality for performance (by default elements in a C array are stored in sequential addresses in memory, and numerical APIs such as PDL or numpy interface with such containers) vis-a-vis containers,
e.g. Perl arrays which do not store their values in sequential addresses in memory. Last, but certainly not least, the C code implementations will allow us to assess whether flags related to floating point operations for the low level compiler (in this case gcc) can affect performance.
This point is worth emphasizing: common mortals are entirely dependent on the choice of compiler flags when "piping" their "install" or building their Inline file. If one does not touch these flags, then one will be blissfully unaware of what they may missing, or pitfalls they may be avoiding.
The humble C file makefile allows one to make such performance evaluations explicitly.

The C code for our toy example is listed in its entirety below. The code is rather self-explanatory, so will not spend time explaining other than pointing out that it contains four functions for

Non-sequential calculation of the expensive function : all three floating pointing operations take place inside a single loop using one thread
Sequential calculations of the expensive function : each of the 3 floating point function evaluations takes inside a separate loop using one thread
Non-sequential OpenMP code : threaded version of the non-sequential code
Sequential OpenMP code: threaded of the sequential code

In this case, one may hope that the compiler is smart enough to recognize that the square root maps to packed (vectorized) floating pointing operations in assembly, so that one function can be vectorized using the appropriate SIMD instructions (note we did not use the simd program for the OpenMP codes).
Perhaps the speedup from the vectorization may offset the loss of performance from repeatedly accessing the same memory locations (or not).


#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <omp.h>

// simulates a large array of random numbers
double*  simulate_array(int num_of_elements,int seed);
// OMP environment functions
void _set_openmp_schedule_from_env();
void _set_num_threads_from_env();



// functions to modify C arrays 
void map_c_array(double* array, int len);
void map_c_array_sequential(double* array, int len);
void map_C_array_using_OMP(double* array, int len);
void map_C_array_sequential_using_OMP(double* array, int len);

int main(int argc, char *argv[]) {
    if (argc != 2) {
        printf("Usage: %s <array_size>\n", argv[0]);
        return 1;
    }

    int array_size = atoi(argv[1]);
    // printf the array size
    printf("Array size: %d\n", array_size);
    double *array = simulate_array(array_size, 1234);

    // Set OMP environment
    _set_openmp_schedule_from_env();
    _set_num_threads_from_env();

    // Perform calculations and collect timing data
    double start_time, end_time, elapsed_time;
    // Non-Sequential calculation
    start_time = omp_get_wtime();
    map_c_array(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Non-sequential calculation time: %f seconds\n", elapsed_time);
    free(array);

    // Sequential calculation
    array = simulate_array(array_size, 1234);
    start_time = omp_get_wtime();
    map_c_array_sequential(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Sequential calculation time: %f seconds\n", elapsed_time);
    free(array);

    array = simulate_array(array_size, 1234);
    // Parallel calculation using OMP
    start_time = omp_get_wtime();
    map_C_array_using_OMP(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Parallel calculation using OMP time: %f seconds\n", elapsed_time);
    free(array);

    // Sequential calculation using OMP
    array = simulate_array(array_size, 1234);
    start_time = omp_get_wtime();
    map_C_array_sequential_using_OMP(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Sequential calculation using OMP time: %f seconds\n", elapsed_time);

    free(array);
    return 0;
}



/*
*******************************************************************************
* OMP environment functions
*******************************************************************************
*/
void _set_openmp_schedule_from_env() {
  char *schedule_env = getenv("OMP_SCHEDULE");
  printf("Schedule from env %s\n", getenv("OMP_SCHEDULE"));
  if (schedule_env != NULL) {
    char *kind_str = strtok(schedule_env, ",");
    char *chunk_size_str = strtok(NULL, ",");

    omp_sched_t kind;
    if (strcmp(kind_str, "static") == 0) {
      kind = omp_sched_static;
    } else if (strcmp(kind_str, "dynamic") == 0) {
      kind = omp_sched_dynamic;
    } else if (strcmp(kind_str, "guided") == 0) {
      kind = omp_sched_guided;
    } else {
      kind = omp_sched_auto;
    }
    int chunk_size = atoi(chunk_size_str);
    omp_set_schedule(kind, chunk_size);
  }
}

void _set_num_threads_from_env() {
  char *num = getenv("OMP_NUM_THREADS");
  printf("Number of threads = %s from within C\n", num);
  omp_set_num_threads(atoi(num));
}
/*
*******************************************************************************
* Functions that modify C arrays whose address is passed from Perl in C
*******************************************************************************
*/

double*  simulate_array(int num_of_elements, int seed) {
  srand(seed); // Seed the random number generator
  double *array = (double *)malloc(num_of_elements * sizeof(double));
  for (int i = 0; i < num_of_elements; i++) {
    array[i] =
        (double)rand() / RAND_MAX; // Generate a random double between 0 and 1
  }
  return array;
}

void map_c_array(double *array, int len) {
  for (int i = 0; i < len; i++) {
    array[i] = cos(sin(sqrt(array[i])));
  }
}

void map_c_array_sequential(double* array, int len) {
  for (int i = 0; i < len; i++) {
    array[i] = sqrt(array[i]);
  }
  for (int i = 0; i < len; i++) {
    array[i] = sin(array[i]);
  }
  for (int i = 0; i < len; i++) {
    array[i] = cos(array[i]);
  }
}

void map_C_array_using_OMP(double* array, int len) {
#pragma omp parallel
  {
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = cos(sin(sqrt(array[i])));
    }
  }
}

void map_C_array_sequential_using_OMP(double* array, int len) {
#pragma omp parallel
  {
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = sqrt(array[i]);
    }
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = sin(array[i]);
    }
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = cos(array[i]);
    }
  }
}

A critical question is whether the use of fast floating compiler flags, a trick that trades speed for accuracy of the code, can affect performance.
Here is the makefile withut this compiler flag

CC = gcc
CFLAGS = -O3 -ftree-vectorize  -march=native  -Wall -std=gnu11 -fopenmp -fstrict-aliasing 
LDFLAGS = -fPIE -fopenmp
LIBS =  -lm

SOURCES = inplace_array_mod_with_OpenMP.c
OBJECTS = $(SOURCES:.c=_noffmath_gcc.o)
EXECUTABLE = inplace_array_mod_with_OpenMP_noffmath_gcc

all: $(SOURCES) $(EXECUTABLE)

clean:
    rm -f $(OBJECTS) $(EXECUTABLE)

$(EXECUTABLE): $(OBJECTS)
    $(CC) $(LDFLAGS) $(OBJECTS) $(LIBS) -o $@

%_noffmath_gcc.o : %.c 
    $(CC) $(CFLAGS) -c $< -o $@

and here is the one with this flag:

CC = gcc
CFLAGS = -O3 -ftree-vectorize  -march=native -Wall -std=gnu11 -fopenmp -fstrict-aliasing -ffast-math
LDFLAGS = -fPIE -fopenmp
LIBS =  -lm

SOURCES = inplace_array_mod_with_OpenMP.c
OBJECTS = $(SOURCES:.c=_gcc.o)
EXECUTABLE = inplace_array_mod_with_OpenMP_gcc

all: $(SOURCES) $(EXECUTABLE)

clean:
    rm -f $(OBJECTS) $(EXECUTABLE)

$(EXECUTABLE): $(OBJECTS)
    $(CC) $(LDFLAGS) $(OBJECTS) $(LIBS) -o $@

%_gcc.o : %.c 
    $(CC) $(CFLAGS) -c $< -o $@

And here are the results of running these two programs

Without -ffast-math

OMP_SCHEDULE=guided,1 OMP_NUM_THREADS=8 ./inplace_array_mod_with_OpenMP_noffmath_gcc 50000000
Array size: 50000000
Schedule from env guided,1
Number of threads = 8 from within C
Non-sequential calculation time: 1.12 seconds
Sequential calculation time: 0.95 seconds
Parallel calculation using OMP time: 0.17 seconds
Sequential calculation using OMP time: 0.15 seconds

With -ffast-math

OMP_SCHEDULE=guided,1 OMP_NUM_THREADS=8 ./inplace_array_mod_with_OpenMP_gcc 50000000
Array size: 50000000
Schedule from env guided,1
Number of threads = 8 from within C
Non-sequential calculation time: 0.27 seconds
Sequential calculation time: 0.28 seconds
Parallel calculation using OMP time: 0.05 seconds
Sequential calculation using OMP time: 0.06 seconds

Note that one can use the fastmath in Numba code as follows (the default is fastmath=False):

@njit(nogil=True,fastmath=True)
def compute_inplace_with_numba(array):
    np.sqrt(array,array)
    np.sin(array,array)
    np.cos(array,array)

A few points that are worth noting:

The -ffast-math gives major boost in performance (about 300% for both the single threaded and the multi-threaded code), but it can generate erroneous results
Fastmath also works in Numba, but should be avoided for the same reasons it should be avoided in any application that strives for accuracy
The sequential C single threaded code gives performance similar to the single threaded PDL and Numpy
Somewhat surprisingly, the sequential code is about 20% faster than the non-sequential code when the correct (non-fast) math is used.
Unsurprisingly, multi-threaded code is faster than single threaded code :)
I still cannot explain how numbas delivers a 50% performance premium over the C code of this rather simple function.

title: " The Quest for Performance Part III : C Force "

date: 2024-07-07

The C code for our toy example is listed in its entirety below. The code is rather self-explanatory, so will not spend time explaining other than pointing out that it contains four functions for

Non-sequential calculation of the expensive function : all three floating pointing operations take place inside a single loop using one thread
Sequential calculations of the expensive function : each of the 3 floating point function evaluations takes inside a separate loop using one thread
Non-sequential OpenMP code : threaded version of the non-sequential code
Sequential OpenMP code: threaded of the sequential code


#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <omp.h>

// simulates a large array of random numbers
double*  simulate_array(int num_of_elements,int seed);
// OMP environment functions
void _set_openmp_schedule_from_env();
void _set_num_threads_from_env();



// functions to modify C arrays 
void map_c_array(double* array, int len);
void map_c_array_sequential(double* array, int len);
void map_C_array_using_OMP(double* array, int len);
void map_C_array_sequential_using_OMP(double* array, int len);

int main(int argc, char *argv[]) {
    if (argc != 2) {
        printf("Usage: %s <array_size>\n", argv[0]);
        return 1;
    }

    int array_size = atoi(argv[1]);
    // printf the array size
    printf("Array size: %d\n", array_size);
    double *array = simulate_array(array_size, 1234);

    // Set OMP environment
    _set_openmp_schedule_from_env();
    _set_num_threads_from_env();

    // Perform calculations and collect timing data
    double start_time, end_time, elapsed_time;
    // Non-Sequential calculation
    start_time = omp_get_wtime();
    map_c_array(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Non-sequential calculation time: %f seconds\n", elapsed_time);
    free(array);

    // Sequential calculation
    array = simulate_array(array_size, 1234);
    start_time = omp_get_wtime();
    map_c_array_sequential(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Sequential calculation time: %f seconds\n", elapsed_time);
    free(array);

    array = simulate_array(array_size, 1234);
    // Parallel calculation using OMP
    start_time = omp_get_wtime();
    map_C_array_using_OMP(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Parallel calculation using OMP time: %f seconds\n", elapsed_time);
    free(array);

    // Sequential calculation using OMP
    array = simulate_array(array_size, 1234);
    start_time = omp_get_wtime();
    map_C_array_sequential_using_OMP(array, array_size);
    end_time = omp_get_wtime();
    elapsed_time = end_time - start_time;
    printf("Sequential calculation using OMP time: %f seconds\n", elapsed_time);

    free(array);
    return 0;
}



/*
*******************************************************************************
* OMP environment functions
*******************************************************************************
*/
void _set_openmp_schedule_from_env() {
  char *schedule_env = getenv("OMP_SCHEDULE");
  printf("Schedule from env %s\n", getenv("OMP_SCHEDULE"));
  if (schedule_env != NULL) {
    char *kind_str = strtok(schedule_env, ",");
    char *chunk_size_str = strtok(NULL, ",");

    omp_sched_t kind;
    if (strcmp(kind_str, "static") == 0) {
      kind = omp_sched_static;
    } else if (strcmp(kind_str, "dynamic") == 0) {
      kind = omp_sched_dynamic;
    } else if (strcmp(kind_str, "guided") == 0) {
      kind = omp_sched_guided;
    } else {
      kind = omp_sched_auto;
    }
    int chunk_size = atoi(chunk_size_str);
    omp_set_schedule(kind, chunk_size);
  }
}

void _set_num_threads_from_env() {
  char *num = getenv("OMP_NUM_THREADS");
  printf("Number of threads = %s from within C\n", num);
  omp_set_num_threads(atoi(num));
}
/*
*******************************************************************************
* Functions that modify C arrays whose address is passed from Perl in C
*******************************************************************************
*/

double*  simulate_array(int num_of_elements, int seed) {
  srand(seed); // Seed the random number generator
  double *array = (double *)malloc(num_of_elements * sizeof(double));
  for (int i = 0; i < num_of_elements; i++) {
    array[i] =
        (double)rand() / RAND_MAX; // Generate a random double between 0 and 1
  }
  return array;
}

void map_c_array(double *array, int len) {
  for (int i = 0; i < len; i++) {
    array[i] = cos(sin(sqrt(array[i])));
  }
}

void map_c_array_sequential(double* array, int len) {
  for (int i = 0; i < len; i++) {
    array[i] = sqrt(array[i]);
  }
  for (int i = 0; i < len; i++) {
    array[i] = sin(array[i]);
  }
  for (int i = 0; i < len; i++) {
    array[i] = cos(array[i]);
  }
}

void map_C_array_using_OMP(double* array, int len) {
#pragma omp parallel
  {
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = cos(sin(sqrt(array[i])));
    }
  }
}

void map_C_array_sequential_using_OMP(double* array, int len) {
#pragma omp parallel
  {
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = sqrt(array[i]);
    }
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = sin(array[i]);
    }
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      array[i] = cos(array[i]);
    }
  }
}

A critical question is whether the use of fast floating compiler flags, a trick that trades speed for accuracy of the code, can affect performance.
Here is the makefile withut this compiler flag

CC = gcc
CFLAGS = -O3 -ftree-vectorize  -march=native  -Wall -std=gnu11 -fopenmp -fstrict-aliasing 
LDFLAGS = -fPIE -fopenmp
LIBS =  -lm

SOURCES = inplace_array_mod_with_OpenMP.c
OBJECTS = $(SOURCES:.c=_noffmath_gcc.o)
EXECUTABLE = inplace_array_mod_with_OpenMP_noffmath_gcc

all: $(SOURCES) $(EXECUTABLE)

clean:
    rm -f $(OBJECTS) $(EXECUTABLE)

$(EXECUTABLE): $(OBJECTS)
    $(CC) $(LDFLAGS) $(OBJECTS) $(LIBS) -o $@

%_noffmath_gcc.o : %.c 
    $(CC) $(CFLAGS) -c $< -o $@

and here is the one with this flag:

CC = gcc
CFLAGS = -O3 -ftree-vectorize  -march=native -Wall -std=gnu11 -fopenmp -fstrict-aliasing -ffast-math
LDFLAGS = -fPIE -fopenmp
LIBS =  -lm

SOURCES = inplace_array_mod_with_OpenMP.c
OBJECTS = $(SOURCES:.c=_gcc.o)
EXECUTABLE = inplace_array_mod_with_OpenMP_gcc

all: $(SOURCES) $(EXECUTABLE)

clean:
    rm -f $(OBJECTS) $(EXECUTABLE)

$(EXECUTABLE): $(OBJECTS)
    $(CC) $(LDFLAGS) $(OBJECTS) $(LIBS) -o $@

%_gcc.o : %.c 
    $(CC) $(CFLAGS) -c $< -o $@

And here are the results of running these two programs

Without -ffast-math

OMP_SCHEDULE=guided,1 OMP_NUM_THREADS=8 ./inplace_array_mod_with_OpenMP_noffmath_gcc 50000000
Array size: 50000000
Schedule from env guided,1
Number of threads = 8 from within C
Non-sequential calculation time: 1.12 seconds
Sequential calculation time: 0.95 seconds
Parallel calculation using OMP time: 0.17 seconds
Sequential calculation using OMP time: 0.15 seconds

With -ffast-math

OMP_SCHEDULE=guided,1 OMP_NUM_THREADS=8 ./inplace_array_mod_with_OpenMP_gcc 50000000
Array size: 50000000
Schedule from env guided,1
Number of threads = 8 from within C
Non-sequential calculation time: 0.27 seconds
Sequential calculation time: 0.28 seconds
Parallel calculation using OMP time: 0.05 seconds
Sequential calculation using OMP time: 0.06 seconds

Note that one can use the fastmath in Numba code as follows (the default is fastmath=False):

@njit(nogil=True,fastmath=True)
def compute_inplace_with_numba(array):
    np.sqrt(array,array)
    np.sin(array,array)
    np.cos(array,array)

A few points that are worth noting:

The -ffast-math gives major boost in performance (about 300% for both the single threaded and the multi-threaded code), but it can generate erroneous results
Fastmath also works in Numba, but should be avoided for the same reasons it should be avoided in any application that strives for accuracy
The sequential C single threaded code gives performance similar to the single threaded PDL and Numpy
Somewhat surprisingly, the sequential code is about 20% faster than the non-sequential code when the correct (non-fast) math is used.
Unsurprisingly, multi-threaded code is faster than single threaded code :)
I still cannot explain how numbas delivers a 50% performance premium over the C code of this rather simple function.

The Quest for Performance Part II : Perl vs Python

chrisarg — Wed, 31 Jul 2024 03:35:28 +0000

Having run a toy performance example, we will now digress somewhat and contrast the performance against
a few Python implementations. First let's set up the stage for the calculations, and provide commandline
capabilities to the Python script.

import argparse
import time
import math
import numpy as np
import os
from numba import njit
from joblib import Parallel, delayed

parser = argparse.ArgumentParser()
parser.add_argument("--workers", type=int, default=8)
parser.add_argument("--arraysize", type=int, default=100_000_000)
args = parser.parse_args()
# Set the number of threads to 1 for different libraries
print("=" * 80)
print(
    f"\nStarting the benchmark for {args.arraysize} elements "
    f"using {args.workers} threads/workers\n"
)

# Generate the data structures for the benchmark
array0 = [np.random.rand() for _ in range(args.arraysize)]
array1 = array0.copy()
array2 = array0.copy()
array_in_np = np.array(array1)
array_in_np_copy = array_in_np.copy()

And here are our contestants:

Base Python

  for i in range(len(array0)):
    array0[i] = math.cos(math.sin(math.sqrt(array0[i])))

Numpy (Single threaded)

np.sqrt(array_in_np, out=array_in_np)
np.sin(array_in_np, out=array_in_np)
np.cos(array_in_np, out=array_in_np)

Joblib (note that this example is not a true in-place one, but I have not been able to make it run using the out arguments)


def compute_inplace_with_joblib(chunk):
    return np.cos(np.sin(np.sqrt(chunk))) #parallel function for joblib

chunks = np.array_split(array1, args.workers)  # Split the array into chunks
numresults = Parallel(n_jobs=args.workers)(
        delayed(compute_inplace_with_joblib)(chunk) for chunk in chunks
    )# Process each chunk in a separate thread
array1 = np.concatenate(numresults)  # Concatenate the results

Numba

@njit
def compute_inplace_with_numba(array):
    np.sqrt(array,array)
    np.sin(array,array)
    np.cos(array,array)
    ## njit will compile this function to machine code
compute_inplace_with_numba(array_in_np_copy)

And here are the timing results:

In place in (  base Python): 11.42 seconds
In place in (Python Joblib): 4.59 seconds
In place in ( Python Numba): 2.62 seconds
In place in ( Python Numpy): 0.92 seconds

The numba is surprisingly slower!? Could it be due to the overhead of compilation as pointed out by mohawk2 in an IRC exchange about this issue?
To test this, we should call compute_inplace_with_numba once before we execute the benchmark. Doing so, shows that Numba is now faster than Numpy.

In place in (  base Python): 11.89 seconds
In place in (Python Joblib): 4.42 seconds
In place in ( Python Numpy): 0.93 seconds
In place in ( Python Numba): 0.49 seconds

Finally, I decided to take base R for ride in the same example:

n<-50000000
x<-runif(n)
start_time <- Sys.time()
result <- cos(sin(sqrt(x)))
end_time <- Sys.time()

# Calculate the time taken
time_taken <- end_time - start_time

# Print the time taken
print(sprintf("Time in base R: %.2f seconds", time_taken))

which yielded the following timing result:

Time in base R: 1.30 seconds

Compared to the Perl results we note the following about this example:

Inplace operations in base Python were ~ 3.5 slower than Perl
Single threaded PDL and numpy gave nearly identical results, followed closely by base R
Failure to account for the compilation overhead of Numba yields the false impression that it is slower than Numpy. When accounting for the compilation overhead, Numba is x2 faster than Numpy
Parallelization with Joblib did improve upon base Python, but was still inferior to the single thread Perl implementation
Multi-threaded PDL (and OpenMP) crushed (not crashed!) every other implementation in all lanugages). Hopefully this post provides some food for thought about the language to use for your next data/compute intensive operation. The next part in this series will look into the same example using arrays in C. This final installment will (hopefully) provide some insights about the impact of memory locality and the overhead incurred by using dynamically typed languages.

The Quest for Performance Part I : Inline C, OpenMP and the Perl Data Language (PDL)

chrisarg — Fri, 26 Jul 2024 18:58:39 +0000

Sometimes, one's code must simply perform and principles, such as aesthetics, "cleverness" or commitment to a single language solution simply go out of the window.
At the TPRC I gave a talk (here are the slides) about how this
can be done for bioinformatics applications, but I feel that a simpler example is warranted to illustrate the potential venues to maximize performance that a Perl
programmer has at their disposal when working in data intensive applications.

So here is a toy problem to illustrate these options. Given a very large array of double precision floats transform them in place with the following function : cos(sin(sqrt(x))).
The function has 3 nested floating point operations. This is an expensive function to evaluate, especially if one has to calculate for a large number of values. We can generate reasonably
quickly the array values in Perl (and some copies for the solutions we will be examining) using the following code:

my $num_of_elements = 50_000_000;
my @array0 = map { rand } 1 .. $num_of_elements;    ## generate random numbers
my @array1 = @array0;                               ## copy the array
my @array2 = @array0;                               ## another copy
my @array3 = @array0;                               ## yet another copy
my @rray4  = @array0;                               ## the last? copy
my $array_in_PDL      = pdl(@array0);    ## convert the array to a PDL ndarray
my $array_in_PDL_copy = $array_in_PDL->copy;    ## copy the PDL ndarray

The posssible solutions include the following:

Inplace modification using a for-loop in Perl.

for my $elem (@array0) {
    $elem = cos( sin( sqrt($elem) ) );
}

Using Inline C code to walk the array and transform in place in C. . Effectively one does a inplace map using C. Accessing elements of Perl arrays (AV* in C) in C is particularly
performant if one is using perl 5.36 and above because of an optimized fetch function introduced in that version of Perl.

void map_in_C(AV *array) {
  int len = av_len(array) + 1;
  for (int i = 0; i < len; i++) {
    SV **elem = av_fetch_simple(array, i, 0); // perl 5.36 and above
    if (elem != NULL) {
      double value = SvNV(*elem);
      value = cos(sin(sqrt(value))); // Modify the value
      sv_setnv(*elem, value);
    }
  }
}

Using Inline C code to transform the array, but break the transformation in 3 sequential C for-loops. This is an experiment really about tradeoffs: modern x86 processors have a specialized,
vectorized square root instruction, so perhaps the compiler can figure how to use it to speed up at least one part of the calculation. On the other hand, we will be reducing the arithmetic intensity
of each loop and accessing the same data value twice, so there will likely be a price to pay for these repeated data accesses.

void map_in_C_sequential(AV *array) {
  int len = av_len(array) + 1;
  for (int i = 0; i < len; i++) {
    SV **elem = av_fetch_simple(array, i, 0); // perl 5.36 and above
    if (elem != NULL) {
      double value = SvNV(*elem);
      value = sqrt(value); // Modify the value
      sv_setnv(*elem, value);
    }
  }
  for (int i = 0; i < len; i++) {
    SV **elem = av_fetch_simple(array, i, 0); // perl 5.36 and above
    double value = SvNV(*elem);
    value = sin(value); // Modify the value
    sv_setnv(*elem, value);
  }
  for (int i = 0; i < len; i++) {
    SV **elem = av_fetch_simple(array, i, 0); // perl 5.36 and above
    double value = SvNV(*elem);
    value = cos(value); // Modify the value
    sv_setnv(*elem, value);
  }
}

Parallelize the C function loop using OpenMP. In a previous entry we discussed how to control the OpenMP environment from within Perl and compile OpenMP aware Inline::C code for
use by Perl, so let's put this knowledge into action! On the Perl side of the program we will do this:

use v5.38;
use Alien::OpenMP;
use OpenMP::Environment;
use Inline (
    C    => 'DATA',
    with => qw/Alien::OpenMP/,
);
my $env = OpenMP::Environment->new();
my $threads_or_workers = 8; ## or any other value
## modify number of threads and make C aware of the change
$env->omp_num_threads($threads_or_workers);
_set_num_threads_from_env();

## modify runtime schedule and make C aware of the change
$env->omp_schedule("guided,1");    ## modify runtime schedule
_set_openmp_schedule_from_env();

On the C part of the program, we will then do this (the helper functions for the OpenMP environment have been discussed
previously, and thus not repeated here).

#include <omp.h>
void map_in_C_using_OMP(AV *array) {
  int len = av_len(array) + 1;
#pragma omp parallel
  {
#pragma omp for schedule(runtime) nowait
    for (int i = 0; i < len; i++) {
      SV **elem = av_fetch_simple(array, i, 0); // perl 5.36 and above
      if (elem != NULL) {
        double value = SvNV(*elem);
        value = cos(sin(sqrt(value))); // Modify the value
        sv_setnv(*elem, value);
      }
    }
  }
}

Perl Data Language (PDL) to the rescue. The PDL set of modules is yet another way to speed up operations and can save the programmer from C. It also autoparallelizes given the right directives so why not use it?

use PDL;
## set the minimum size problem for autothreading in PDL
set_autopthread_size(0);
my $threads_or_workers = 8; ## or any other value

## PDL
## use PDL to modify the array - multi threaded
set_autopthread_targ($threads_or_workers);
$array_in_PDL->inplace->sqrt;
$array_in_PDL->inplace->sin;
$array_in_PDL->inplace->cos;


## use PDL to modify the array - single thread
set_autopthread_targ(0);

$array_in_PDL_copy->inplace->sqrt;
$array_in_PDL_copy->inplace->sin;
$array_in_PDL_copy->inplace->cos;

Using 8 threads we get something like this

Inplace benchmarks
Inplace  in         Perl took 2.85 seconds
Inplace  in Perl/mapCseq took 1.62 seconds
Inplace  in    Perl/mapC took 1.54 seconds
Inplace  in   Perl/C/OMP took 0.24 seconds

PDL benchmarks
Inplace  in     PDL - ST took 0.94 seconds
Inplace  in     PDL - MT took 0.17 seconds

Using 16 threads we get this!

Starting the benchmark for 50000000 elements using 16 threads/workers

Inplace benchmarks
Inplace  in         Perl took 3.00 seconds
Inplace  in Perl/mapCseq took 1.72 seconds
Inplace  in    Perl/mapC took 1.62 seconds
Inplace  in   Perl/C/OMP took 0.13 seconds

PDL benchmarks
Inplace  in     PDL - ST took 0.99 seconds
Inplace  in     PDL - MT took 0.10 seconds

A few observations:

The OpenMP and the multi-threaded (MT) of the PDL are responsive to the number of workers, while the solutions are not. Hence, the timings of the pure Perl and the inline non-OpenMP solution timings in these benchmarks give an idea of the natural variability in performance
Writing the map version of the code in C improved performance by about 180% (contrast Perl and Perl/mapC).
Using PDL in a single thread improved performance by 285-300% (contrast PDL - ST and Perl timings).
There was a price to pay for repeated memory access (contrast Perl/mapC to Perl/mapCseq)
OpenMP and multi-threaded PDL operations gave similar performance (though PDL appeared faster in these examples). The code run 23-30 times faster.

In summary, there are both native (PDL modules) and foreign (C/OpenMP) solutions to speed up data intensive operations in Perl, so why not use them widely and wisely to make Perl programs performant?

Marrying Perl to Assembly

chrisarg — Thu, 04 Jul 2024 18:34:13 +0000

This is probably one of the things that should never be allowed to exist, but why not use Perl and its capabilities to inline foreign code, to FAFO with assembly without a build system? Everything in a single file! In the process one may find ways to use Perl to enhance NASM and vice versa. But for now, I make no such claims : I am just using the perlAssembly git repo to illustrate how one can use Perl to drive (and learn to code!) assembly programs from a single file.
(Source code may be found in the perlAssembly repo )

x86-64 examples

Adding Two Integers

Simple integer addition in Perl - this is the Hello World version of the perlAssembly repo
But if we can add two numbers, why not add many, many more?

The sum of an array of integers

Explore multiple equivalent ways to add large arrays of short integers (e.g. between -100 to 100) in Perl. The Perl and the C source files contain the code for:

ASM_blank : tests the speed of calling ASM from Perl (no computations are done)
ASM : passes the integers as bytes and then uses conversion operations and scalar floating point addition
ASM_doubles : passes the array as a packed string of doubles and do scalar double floating addition in assembly
ASM_doubles_AVX: passes the array as a packed string of doubles and do packed floating point addition in assembly
ForLoop : standard for loop in Perl
ListUtil: sum function from list utilities
PDL : uses summation in PDL

Scenarios w_alloc : allocate memory for each iteration to test the speed of pack, those marked
as wo_alloc, use a pre-computed data structure to pass the array to the underlying code.
Benchmarks of the first scenario give the true cost of offloading summation to of a Perl array to a given
function when the source data are in Perl. Timing the second scenario benchmarks speed of the
underlying implementation.

This example illustrates

an important (but not the only one!) strategy to create a data structure that is suitable for Assembly to work with, i.e. a standard array of the appropriate type, in which one element is laid adjacent to the previous one in memory
the emulation of declaring a pointer as constant in the interface of a C function. In the AVX code, we don't FAFO with the pointer (RSI in the calling convention) to the array directly, but first load its address to another register that we manipulate at will.

Results

Here are the timings!

	mean	median	stddev
ASM_blank	2.3e-06	2.0e-06	1.1e-06
ASM_doubles_AVX_w_alloc	3.6e-03	3.5e-03	4.2e-04
ASM_doubles_AVX_wo_alloc	3.0e-04	2.9e-04	2.7e-05
ASM_doubles_w_alloc	4.3e-03	4.1e-03	4.5e-04
ASM_doubles_wo_alloc	8.9e-04	8.7e-04	3.0e-05
ASM_w_alloc	4.3e-03	4.2e-03	4.5e-04
ASM_wo_alloc	9.2e-04	9.1e-04	4.1e-05
ForLoop	1.9e-02	1.9e-02	2.6e-04
ListUtil	4.5e-03	4.5e-03	1.4e-04
PDL_w_alloc	2.1e-02	2.1e-02	6.7e-04
PDL_wo_alloc	9.2e-04	9.0e-04	3.9e-05

Let's say we wanted to do this toy experiment in pure C (using Inline::C of course!)
This code obtains the integers as a packed "string" of doubles and forms the sum in C

double sum_array_C(char *array_in, size_t length) {
    double sum = 0.0;
    double * array = (double *) array_in;
    for (size_t i = 0; i < length; i++) {
        sum += array[i];
    }
    return sum;
}

Here are the timing results:

	mean	median	stddev
C_doubles_w_alloc	4.1e-03	4.1e-03	2.3e-04
C_doubles_wo_alloc	9.0e-04	8.7e-04	4.6e-05

What if we used SIMD directives and parallel loop constructs in OpenMP? All three combinations were tested, i.e. SIMD directives
alone (the C equivalent of the AVX code), OpenMP parallel loop threads and SIMD+OpenMP.
Here are the timings!

	mean	median	stddev
C_OMP_w_alloc	4.0e-03	3.7e-03	1.4e-03
C_OMP_wo_alloc	3.1e-04	2.3e-04	9.5e-04
C_SIMD_OMP_w_alloc	4.0e-03	3.8e-03	8.6e-04
C_SIMD_OMP_wo_alloc	3.1e-04	2.5e-04	8.5e-04
C_SIMD_w_alloc	4.1e-03	4.0e-03	2.4e-04
C_SIMD_wo_alloc	5.0e-04	5.0e-04	8.9e-05

Discussion of the sum of an array of integers example

For calculations such as this, the price that must be paid is all in memory currency: it takes time to generate these large arrays, and for code with low arithmetic intensity this time dominates the numeric calculation time.
Look how insanely effective sum in List::Util is : even though it has to walk the Perl array whose elements (the doubles, not the AV*) are not stored in a contiguous area in memory, it is no more than 3x slower than the equivalent C code C_doubles_wo_alloc.
Look how optimized PDL is compared to the C code in the scenario without memory allocation.
Manual SIMD coded in assembly is 40% faster than the equivalent SIMD code in OpenMP (but it is much more painful to write)
The threaded OpenMP version achieved equivalent performance to the single thread AVX assembly programs, with no obvious improvement from combining SIMD+parallel loop for pragmas in OpenMP.
For the example considered here, it thus makes ZERO senso to offload a calculation as simple as a summation because ListUtil is already within 15% of the assembly solution (at a latter iteration we will also test AVX2 and AVX512 packed addition to see if we can improve the results).
If however, one was managing the array, not as a Perl array, but as an area in memory through a Perl object, then one COULD consider offloading. It may be fun to consider an example in which one adds the output of a function that has an efficient PDL and assembly implementation to see how the calculus changes (in the to-do list for now).

Disclaimer

The code here is NOT meant to be portable. I code in Linux and in x86-64, so if you are looking into Window's ABI or ARM, you will be disappointed. But as my knowledge of ARM assembly grows, I intend to rewrite some examples in Arm assembly!