Multicore matrix multiplication using PDL::Parallel::threads

matmult_pdl_thr.pl

#!/usr/bin/env perl

##
## Usage:
##    perl matmult_pdl_thr.pl   1024  ## Default size is 512:  $c = $a * $b
##
## by David Mertens
## based on code by Mario Roy

#################
# Preliminaries #
#################

use strict;
use warnings;
 
my $prog_name = $0; $prog_name =~ s{^.*[\\/]}{}g;
 
use Time::HiRes qw(time);
 
use PDL;
use PDL::Parallel::threads qw(retrieve_pdls);
use PDL::Parallel::threads::SIMD qw(parallel_id parallelize);


# Get the matrix size and croak on bad input
my $tam = @ARGV ? shift : 512;
die "Error: $tam must be an integer greater than 1.\n"
  if $tam !~ /^\d+$/ or $tam < 2;

my $cols = $tam;
my $rows = $tam;
 
###########################
# Create some shared data #
###########################

sequence($cols, $rows)->share_as('left_input');
sequence($rows, $cols)->share_as('right_input');
my $output = zeroes($rows, $rows)->share_as('output');
my $N_threads = 4;

###################################
# Run the calculation in parallel #
###################################

my $start = time();
parallelize {
	my ($l, $r, $o) = retrieve_pdls('left_input', 'right_input', 'output');
	my $pid = parallel_id;
	
	# chop up the input matrix based on the number of rows and the number
	# of threads.
	my $step = int($rows / $N_threads);
	my $start = $pid * $step;
	my $stop = ($pid + 1) * $step - 1;
	$stop = $rows - 1 if $stop >= $rows;
	
	use PDL::NiceSlice;
	$o(:, $start:$stop) .= $l(:,$start:$stop) x $r;
	no PDL::NiceSlice;
} $N_threads;
my $end = time();

#########################
# Print out the results #
#########################

printf "## $prog_name $tam: compute time: %0.03f secs\n\n", $end - $start;

for my $pair ([0, 0], [324, 5], [42, 172], [$tam-1, $tam-1]) {
	my ($row, $col) = @$pair;
	$row %= $rows;
	$col %= $cols;
	print "## ($row, $col): ", $output->at($row, $col), "\n";
}