-
-
Notifications
You must be signed in to change notification settings - Fork 1.3k
GSoC 2026 ‐ Kaustubh Patange
My self Kaustubh Patange, an undergraduate in Computer Science & Engineering at the Walchand Institute of Technology, Solapur. I am passionate about algebra, programming and backend systems. My prior work spans full-stack and backend engineering beyond this I like exploring machine learning and scaling systems under load.
Extending Level 2 and Level 3 BLAS routines for linear algebra #2039 is stdlib's long-standing tracking issue for a complete, native BLAS surface in JavaScript. It enumerates every routine in the reference BLAS specification across Levels 1, 2, and 3.The major goals associated with the project are:
add JavaScript implementations for the remaining Level 2 and Level 3 routines, add C implementations, re-implement BLAS subroutines using the existing lapack-netlib Fortran implementation, write Node.js bindings to allow calling BLAS routines in compiled C/Fortran from JavaScript, complete the blas/base/ndarray/* layer that sits between the low-level strided kernels. The approach I followed was the one set out in the proposal. In each level, the priority was first real, followed by complex datatypes; for real datatypes, double precision first, then single precision, then the generic variations. Starting from reviewing the Level 1 routines during community bonding, I worked up through real Level 2 and Level 3 before the midterm, lifting each routine into the ndarray layer once its underlying kernels were merged and tested.
During the community bonding period, I began coding by focusing on the BLAS Level 1 routines, as outlined in tracking issue #2039, and the ndarray layer; it let me get used to stdlib's conventions on routines whose mathematics I understood first. The complex routines were harder to grasp at first, but after a few packages it became routine. From there I moved up to the real Level 2 routines, then the packed and triangular variants, and the complex Level 2 routines. Alongside this I worked on migrating the option arguments from strings to enumerated constants, work I had underestimated at the start, but which turned out to be the prerequisite for the ndarray layer.
Native Implementation Level 1 Signature:
dzasum( N, x, strideX )Ndarray Implementation Level 1 Signature:
dzasum( N, x, strideX, offsetX )The difference between the two is a single parameter, but it changes what the routine can address. The native form always begins reading at index zero, so operating on a slice of a longer buffer means creating a view first. With offsetX the caller names the starting position directly, and the routine walks from there which means an arbitrary window of an existing array can be processed in place, with no copy and no wrapper object. dzasum was also where I first met stdlib's handling of complex arrays: a Complex128Array is reinterpreted as a Float64Array of twice the length, with strides and offsets doubled to match, so the loop moves pairs of floats and never allocates a complex number. I reviewed both the JavaScript and the C/Fortran implementations of this package, and the Fortran review was the first time I had to read reference code closely enough to judge whether it was faithful rather than merely whether it was tidy.
Moving from Level 1 to Level 2 meant moving from vectors to matrices, and that one change is what makes these routines interesting. A vector needs only a stride; a matrix needs a rule for turning a pair of indices into a single position in a flat buffer, and BLAS carries that rule in the leading dimension. Because the reference implementation is Fortran, it assumes a column-major rule and offers no alternative but JavaScript has no native two-dimensional array either, so a matrix here is still one contiguous block of memory and we are free to choose. Rather than inherit Fortran's assumption, stdlib exposes an order parameter and lets the caller declare which way their data is already laid out, which avoids forcing a transpose on anyone whose data arrives in the other convention.
Taking a small matrix and writing out both conventions:
A = [ 1, 2 ]
[ 3, 4 ] (3X2)
[ 5, 6 ]
A = [ 1, 2, 3, 4, 5, 6 ] (row-major, LDA = 2)
A = [ 1, 3, 5, 2, 4, 6 ] (column-major, LDA = 3)
Level 2 also brings in the option flags that decide which part of the matrix is touched at all, and each of them multiplies the number of paths through the implementation: uplo selects which triangle is referenced, trans selects whether the matrix, its transpose, or its conjugate transpose enters the product, and diag says whether the diagonal is stored or implicitly one. On top of that, the same mathematical matrix can arrive under several different storage conventions — full, packed into a single contiguous triangle, or banded — each of which needs its own index arithmetic even though the operation itself is unchanged. But the more interesting step comes after all of this.
Native Implementation Level 2 Signature:
dsymv( order, uplo, N, α, A, LDA, x, strideX, β, y, strideY )The vectors in this signature already carry an independent stride each, so there is no reason the matrix should not do the same. Replacing order and LDA with a stride per dimension gives:
sa1: stride along the first dimension of matrix A. sa2: stride along the second dimension of matrix A.
and once those exist, order becomes redundant — row-major is simply the case where sa2 is 1, and column-major the case where sa1 is 1. That is what I mean by the more interesting step: the two layouts stop being separate code paths and become two values of the same pair of parameters, so a single loop nest serves both. It also makes the ndarray form strictly more expressive than the one it replaces. Swapping sa1 and sa2 gives a transposed view without moving any data; shifting the offset selects a submatrix of a larger matrix; and a negative stride walks a dimension backwards, none of which order and LDA can describe between them.
Ndarray implementation signature:
zgemv( trans, M, N, α, A, sa1, sa2, offsetA, x, strideX, offsetX, β, y, strideY, offsetY )A short example of how the matrix stride parameters behave in practice. Suppose the values we care about are scattered through a larger 4X4 buffer, stored row-major, with the rest of the positions holding data we must not touch:
A = [ 999.0, 999.0, 999.0, 999.0 ],
[ 999.0, 1.0, 999.0, 2.0 ],
[ 999.0, 999.0, 999.0, 999.0 ],
[ 999.0, 3.0, 999.0, 4.0 ]
A = [ 999.0, 999.0, 999.0, 999.0, 999.0, 1.0, 999.0, 2.0, 999.0, 999.0, 999.0, 999.0, 999.0, 3.0, 999.0, 4.0 ]
Now let's say our required sub-matrix for operation is:
A = [ 1, 2 ],
[ 3, 4 ]
The offset takes us to the first element, and the two strides describe how far to travel to reach the next row and the next column:
sa1 = 8
sa2 = 2
offsetA = 5
Not every package was already in this shape when I started. dsymv still carried order and LDA in its ndarray implementation, and I reviewed the refactor that moved it onto the stride-based form, a useful way to watch the migration happen from the reviewing side rather than the implementing side.
The other Level 2 routines follow the same pattern, but what changes between them is how the matrix is stored and how much of it is read:
DGEMV - where A is a general MXN matrix. DSYMV - where A is an NXN symmetric matrix and only one triangle is referenced. DSPMV - where A is an NXN symmetric matrix supplied in packed form, as a single contiguous triangle. DTRMV - where A is an NXN triangular matrix whose diagonal may be implicitly one and never read.
The packed and triangular variants took the longest. Packing destroys the regular stride pattern, since consecutive rows of the stored triangle have different lengths, so the index mapping is arithmetic rather than a pair of multiplications. The triangular routines brought a compensation, though: a row-major problem over the upper triangle turns out to be exactly a column-major problem over the lower triangle, so recognising that identity lets one loop nest cover four combinations of uplo and layout instead of one.
Running alongside the implementation work was the migration of the option arguments from strings to enumerated constants. These had been validated with boolean predicates such as isMatrixTriangle, which report whether a value is acceptable but hand back nothing, so every caller downstream keeps passing the original string around. The resolver packages accept either a string or the matching enumerated integer and return a canonical value, which folds validation and normalisation into a single call:
u = resolveStr( uplo );
if ( u === null ) {
throw new TypeError( ... );
}I had filed this under cleanup until I started the same routines one layer higher, in blas/base/ndarray/*.
That namespace is the third form each routine takes. Where the two signatures above accept raw typed arrays and describe their geometry through separate stride, offset and dimension parameters, this one accepts stdlib ndarray objects and reads shape, strides, offset and dtype straight off them, so a single argument replaces the whole group:
Ndarray-object implementation signature:
dsymv( [ A, x, y, uplo, alpha, beta ] )
Everything in that list is an ndarray, including alpha, beta and uplo, which arrive as zero-dimensional ones. That uniformity is the point: an argument list of nothing but ndarrays is what stdlib's dispatch machinery can consume, cache, and eventually hand to a compiled kernel. It is also where the enum work stopped looking like housekeeping a flag can only travel as a zero-dimensional int8 ndarray if it has a canonical integer form, and the resolvers are the only thing that produces one. A JavaScript string cannot be passed to a C or WebAssembly routine; an int8 can. Because the geometry now comes from the ndarray itself, an entire class of caller error simply stops being expressible: there is no LDA to get wrong, and a transposed or strided view is described by the object rather than by parameters the caller has to compute.
In the later weeks I built out a family of triangular routines in the extended BLAS namespace routines that are BLAS-shaped but have no counterpart in the reference specification. triu and tril copy the upper or lower triangular part of one matrix into another, while triu2tril and tril2triu reflect one triangle into the opposite triangle of the destination. Each exists in five dtype variants and, like the Level 2 routines, in all three forms:
dtriu( order, M, N, k, A, LDA, B, LDB )
dtriu( M, N, k, A, sa1, sa2, offsetA, B, sb1, sb2, offsetB )
dtriu( [ A, B, k ] )The parameter worth dwelling on is k, the diagonal the copy is measured from. At zero it is the main diagonal; a positive value moves the boundary upwards so only part of the triangle is copied; a negative value moves it downwards so the triangle plus one or more sub-diagonals is copied. That single parameter is what turns four fixed routines into a general banded-region extractor. The reflecting pair earns its place separately: given one stored triangle they reconstruct the full symmetric matrix in a single pass, which is exactly what is needed before handing a matrix to a routine that will not accept packed or triangle-restricted input. Because A and B carry independent strides in the second form, the reflection is not a second pass over the output, writing with sb1 and sb2 swapped is what performs it.
Writing the tests was rarely the shorter half of the work. Because the ndarray form accepts arbitrary strides and offsets, the interesting failures are the ones that only appear away from the defaults, so the test matrix has to sweep both layouts, every combination of the option flags, positive, negative and non-unit strides on each vector, non-zero offsets, submatrix views, and the early-return paths where a dimension or a scalar is zero. Laying that grid out on paper before writing the loop nest repeatedly changed the implementation, because seeing all the cases at once makes the symmetries between them obvious, and a symmetry is what lets one branch do the work of four. Where the multidimensional behaviour was ambiguous I generated fixtures against the reference implementations and cross-checked them against NumPy and SciPy.
That is the arc of the work as I experienced it. A fair amount of it was not new implementation at all but going back over packages that already existed, lifting them onto shared base implementations, replacing the string arguments, filling in the test cases nobody had written, and those were often the changes that unblocked the most downstream.
TODO: include a list of links to all relevant PRs along with a short description for each. For small bug fix PRs, you can group them together and describe them as, e.g., "various maintenance work".
TODO: add a summary of the current state of the project.
TODO: add a summary of what remains left to do. If there is a tracking issue, include a link to that issue.
TODO: add a summary of any unexpected challenges that you faced, along with any lessons learned during the course of your project.
TODO: add a report summary and include any acknowledgments (e.g., shout outs to contributors/maintainers who helped out along the way, etc).