PLA Status Updates

On Friday night we dropped the kid off with his grandparents for a sleepover. With the apartment to ourselves, my wife and I did what we've been wanting to do for months: she went to bed early and I stayed up to program. I started making some real progress in Parrot-Linear-Algebra, and also uncovered some interesting bugs which needed to be fixed.

I added a short file for adding PLA support to programs written in NQP. The file, pla.nqp, can be included into an NQP program like this once it's installed:

INIT { pir::load_bytecode("pla.pbc"); }

That file loads the PLA binary library, stores a global reference to the library, and registers the type names with P6metaclass. That done, you can start writing more natural-looking programs in NQP. I updated the Gaussian Elimination example I wrote a few weeks ago to use this new bootstrap file (and some new methods I added), which means the example can now run without Kakapo or any other dependencies.

Next up, I started prototyping bindings for the library for Rakudo. These are still preliminary, but I also have a working example of using numerical matrices in Rakudo. Within the next few days I will probably create a complete module for Rakudo (Math::LinearAlgebra, or something like that). I'll post more details as I do more work on it.

I added a new item_at method which returns and optionally sets a value at given coordinates in a matrix. This method is common to all my core matrix types. The list of common methods for these matrix types is:

• resize() : pre-allocate size for the matrix, growing (never shrinking) it to hold a certain size
• fill() : Fill a matrix, or a region of a matrix, with a constant value. Automatically resize if necessary
• transpose() : Transpose (swap rows with columns) the matrix lazily
• mem_transpose() : Eagerly transpose the actual memory contents
• iterate_function_inplace() : Execute a function for every element of the matrix, replacing that element with the function result
• iterate_function_external() : Execute a function for every element of the matrix, creating a new matrix with the results. This is similar to Perl's map function.
• initialize_from_array() : Insert values into the matrix from an array
• initialize_from_args() : Similar to the _from_array variant, but initializes the matrix using elements from a slurpy argument list
• get_block() : Return a block, or "submatrix" from the matrix
• set_block() : Set a block in the matrix
• item_at() : New this weekend, gets or sets a value in the matrix at the specified coordinates

There are maybe a handful of other methods I would like to add, but in terms of the core types, this is the standard API (plus a series of VTABLE calls, which are standard). For mathematics types I also have GEMM (the BLAS-based matrix multiply operation), and elementary matrix operations (add rows, swap rows, scale rows). This is not a shabby interface at all, and can start to be used for some real-world applications.

This weekend I also started seeing some very weird errors in the PLA test suite. Tests were running fine, but I was seeing exceptions (and, in one case, a segfault) occur after all tests had passed. This sounded very similar to another problem I've seen in the past. Here's the test that set it off. Can you spot the problem?

method test_METHOD_iterate_function_inplace_TRANSPOSE() {
    my $m := self.fancymatrix2x2();
    $m.transpose();
    my $n := self.matrix2x2(self.fancyvalue(0) * 2, self.fancyvalue(2) * 2,
                            self.fancyvalue(1) * 2, self.fancyvalue(3) * 2);
    my $sub := -> $matrix, $value, $x, $y {
        return ($value * 2);
    };
    $m.iterate_function_inplace($sub);
    assert_equal($m, $n, "external iteration does not respect transpose");
}

What's maddening is that this test has been a problem for months, but never caused a failure. It was silently wrong, probably since the day I wrote it.

See it yet?

The problem is that return statement. What should happen is this: The pointy block creates an anonymous subroutine, which I pass to the iterate_function_inplace method. The method should loop over ever element in our matrix and call that pointy block for each one. The result should be the same matrix with every element multiplied by two.

What actually happens is this: the iterate_function_inplace method, in order to call the callback, must recurse on the C stack into a new runloop function. The pointy block executes in this inferior runloop. However, where things get weird is that return statement. In NQP, returns are performed by constructing and throwing control exceptions. In the case of the pointy block (and I'm not sure whether or not this is a bug), the return statement jumps to the return handler for the test_METHOD_iterate_function_inplace_TRANSPOSE() function, not the anonymous sub.

The sub executes after having only called the callback once, it never hits the final assertion (which, it turns out, would have failed). Control flow continues inside the inner runloop until the end of the test program, then the C runloop function returns, tries to continue executing from that point, and things go to hell.

The solution is really quite simple. Change this:

my $sub := -> $matrix, $value, $x, $y {
    return ($value * 2);
};

into this:

my $sub := sub ($matrix, $value, $x, $y) {
    return ($value * 2);
}; 

Problem solved, and now more tests are legitimately passing.

It's been a productive couple of days in PLA, and I'm hoping I can keep this momentum up in the days ahead. I need to finish implementing my GEMM wrapper for ComplexMatrix2D, and then get started on the Linear Algebra module for Rakudo.