Skip to content

GSoC 2026 ‐ Kaustubh Patange

Jump to bottom
Kaustubh Patange edited this page Aug 21, 2026 · 28 revisions

About me

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.

Project overview

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.

Project recap

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.

Sitting above everything under base is the top-level blas/* namespace, A top-level package takes an ndarray of whatever numeric dtype, works out which base kernel applies, and adds reduction semantics on top, so the same call works on a vector, on every row of a matrix, or across an entire multidimensional stack.

Top-level Implementation Signature:

nrm2( x[, options] )
x = [ 1.0, -2.0, 2.0 ]

nrm2( x )
// returns <ndarray>[ 3.0 ]

The options are what make it a reduction rather than a single computation: dims names which dimensions to reduce over, defaulting to all of them; dtype sets the output dtype; and keepdims decides whether the reduced dimensions survive as singletons, which is what lets the result broadcast back against the input. I implemented nrm2 and asum in this form, and this is also where the earlier design question resolved itself.

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. 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.

Completed work

Merged PRs

1. blas/base/*

The above routine computes the sum of the absolute values of the real and imaginary components of a complex double-precision vector, that is Σ |Re(x[i])| + |Im(x[i])|. The reference LAPACK implementation contains parameters:

dzasum( N, x, strideX )

Native Implementation:

dzasum( N, x, strideX )

Ndarray Implementation:

dzasum.ndarray( N, x, strideX, offsetX )

The complex array is reinterpreted as a real one of twice the length with the stride and offset doubled, so the loop reads pairs of floats rather than constructing complex values.

The above matrix-matrix multiplication routines short-circuit on α = 0, β = 0 and β = 1, so a suite exercising only α = 1, β = 1 never reaches those branches. These added the missing cases and fixed an argument being forwarded incorrectly to the base implementation.

Open PRs

The above routine performs the symmetric rank 2 operation A = α*x*yᵀ + α*y*xᵀ + A, where α is a scalar, x and y are N element vectors, and A is an N by N symmetric matrix, for any array-like input rather than a fixed dtype:

gsyr2( order, uplo, N, α, x, sx, y, sy, A, LDA )

Ndarray Implementation:

gsyr2.ndarray( uplo, N, α, x, sx, ox, y, sy, oy, A, sa1, sa2, oa )

The above two routines are the packed variants, performing the symmetric rank 1 operation A = α*x*xᵀ + A and the matrix-vector operation y = α*A*x + β*y respectively, where A is supplied as AP, a single contiguous vector holding only one triangle:

gspr( order, uplo, N, α, x, sx, AP )
gspmv( order, uplo, N, α, AP, x, sx, β, y, sy )

What separates the generic packages from the typed ones is underneath rather than in the signature. A double-precision package can index a Float64Array directly; a generic one cannot assume the input is a typed array at all, so it resolves accessors for the input and output and reads and writes through those. That costs a call per element, but it is what allows these routines to run over ordinary arrays, complex arrays, and any other array-like stdlib supports and the top-level packages need exactly that in order to dispatch over dtype.

Open PRs reviewed

Level 1:

The above routines compute the dot product of two complex vectors, conjugated and unconjugated respectively xᴴ·y and xᵀ·y. The pair exists because only the conjugated form gives a real, non-negative result when x = y, which is what makes it usable as a norm; the unconjugated form is the one that matches the bilinear definition. They take the same parameters as the real-valued dot product and differ only in returning a complex value:

cdotc( N, x, strideX, y, strideY )

Ndarray Implementation:

cdotc.ndarray( N, x, strideX, offsetX, y, strideY, offsetY )

The above two routines return the index of the element with the largest |Re| + |Im|, single and double precision respectively. They take a vector and nothing else, since there is no scalar and no second operand:

icamax( N, x, strideX )

Ndarray Implementation:

icamax.ndarray( N, x, strideX, offset )

The i prefix marks the return value as an index rather than a value, and where the reference routines return a one-based index, stdlib returns a zero-based one, and -1 for an empty vector.

Level 2:

  • #10237: feat: add blas/base/zgemv The above routine performs the matrix-vector operation y = α*op(A)*x + β*y, where op(A) is A, Aᵀ or Aᴴ, for complex double precision:
zgemv( order, trans, M, N, α, A, LDA, x, strideX, β, y, strideY )

Ndarray Implementation:

zgemv.ndarray( trans, M, N, α, A, sa1, sa2, offsetA, x, strideX, offsetX, β, y, strideY, offsetY )

The above three are refactors rather than new packages, moving the symmetric and packed symmetric matrix-vector products onto the stride-based ndarray form by replacing order and LDA with a stride per dimension.

C and Fortran implementations

Above routine scales a complex single-precision vector x by a real single-precision constant α, in place. It is the mixed-type member of the scaling family: cscal takes a complex α, csscal takes a real one, which avoids two wasted multiplications and a sign flip per element. The reference LAPACK implementation contains parameters:

csscal( N, α, x, strideX )

Native Implementation:

csscal( N, α, x, strideX )

Ndarray Implementation:

csscal.ndarray( N, α, x, strideX, offsetX )

Open PRs reviewed:

2. Top-level blas/*

Unlike everything under base, they take an ndarray and a set of options, dispatch over dtype, and reduce along one or more dimensions.

Merged PRs:

The above two routines compute the L2-norm and the sum of absolute values along one or more ndarray dimensions:

nrm2( x[, options] )
asum( x[, options] )
x = [ 1.0, -2.0, 2.0 ]
 
nrm2( x )
// returns <ndarray>[ 3.0 ]

The options are dims (which dimensions to reduce over), dtype (the output dtype), and keepdims (whether reduced dimensions survive as singletons).

Open PRs:

The above routine scales an ndarray by a constant, following the same shape as the two merged packages above:

scal( x, alpha[, options] )

3. blas/base/ndarray/*

Every package here takes a single array-like argument holding ndarrays, with shape, strides, offset and dtype read from the objects themselves.

Merged PRs

These perform y = α*op(A)*x + β*y, where op(A) is A, Aᵀ or Aᴴ, α and β are scalars, x and y are vectors, and A is an M by N matrix. The corresponding base implementation contains parameters:

dgemv.ndarray( trans, M, N, α, A, sa1, sa2, oa, x, sx, ox, β, y, sy, oy )

Ndarray-object Implementation:

dgemv( [ A, x, y, trans, α, β ] )

A is a two-dimensional ndarray, x and y are one-dimensional, and trans, α and β are zero-dimensional. y is both input and output, and is returned by reference rather than copied.

syr performs the symmetric rank 1 operation A = α*x*xᵀ + A and syr2 the symmetric rank 2 operation A = α*x*yᵀ + α*y*xᵀ + A, where A is an N by N symmetric matrix of which only the triangle named by uplo is referenced and updated. The corresponding base implementations contain parameters:

dsyr.ndarray( uplo, N, α, x, sx, ox, A, sa1, sa2, oa )
dsyr2.ndarray( uplo, N, α, x, sx, ox, y, sy, oy, A, sa1, sa2, oa )

Ndarray-object Implementation:

dsyr( [ A, x, uplo, α ] )
dsyr2( [ A, x, y, uplo, α ] )

Because A is read and written in the same pass, the loop must consume x[i] and y[i] before the i-th row is modified.

dger.ndarray( M, N, α, x, sx, ox, y, sy, oy, A, sa1, sa2, oa )

symv computes y = α*A*x + β*y for a symmetric A. The sp* packages are the packed equivalents, where A is supplied as AP, a single contiguous vector holding only one triangle:

dspmv.ndarray( uplo, N, α, AP, oap, x, sx, ox, β, y, sy, oy )
dspr.ndarray( uplo, N, α, x, sx, ox, AP, sap, oap )

AP can be thought of as a longer global vector containing the packed triangle of some N by N symmetric matrix, with the offset selecting where it begins. Because the stored rows have unequal lengths, the index mapping here is arithmetic rather than a pair of stride multiplications.

trmv computes x = op(A)*x and trsv solves op(A)*x = b, where A is an N by N unit or non-unit, upper or lower triangular matrix. A unit triangular matrix is one whose diagonal elements are all 1, in which case the diagonal is never read:

dtrmv.ndarray( uplo, trans, diag, N, A, sa1, sa2, oa, x, sx, ox )
dtrsv.ndarray( uplo, trans, diag, N, A, sa1, sa2, oa, x, sx, ox )

The Level 3 routines, performing C = α*op(A)*op(B) + β*C, with op(A) an M by K matrix, op(B) a K by N matrix and C an M by N matrix:

dgemm.ndarray( transA, transB, M, N, K, α, A, sa1, sa2, oa, B, sb1, sb2, ob, β, C, sc1, sc2, oc )

Ndarray-object Implementation:

dgemm( [ A, B, C, transA, transB, α, β ] )

Since both A and B carry their own stride pair, a transpose is a stride swap performed before the kernel runs rather than a separate loop nest.

Level 1 packages, all following the f( [ x, ... ] ) form:

dsdot and sdsdot both take single-precision inputs but accumulate the dot product in double precision, accumulating a long single-precision sum in single precision loses badly to rounding.

Open PRs

rot applies a plane rotation to a pair of vectors; rotm applies a modified Givens rotation.

4. blas/ext/base/*

Merged PRs

triu copies the upper triangular part of a matrix A into another matrix B, and tril copies the lower triangular part. k is the diagonal the copy is measured from: k = 0 is the main diagonal, k > 0 a super-diagonal, and k < 0 a sub-diagonal, so a negative k for triu copies the upper triangle plus one or more sub-diagonals.

Native Implementation:

dtriu( order, M, N, k, A, LDA, B, LDB )
dtril( order, M, N, k, A, LDA, B, LDB )

Ndarray Implementation:

dtriu.ndarray( M, N, k, A, sa1, sa2, offsetA, B, sb1, sb2, offsetB )
dtril.ndarray( M, N, k, A, sa1, sa2, offsetA, B, sb1, sb2, offsetB )

triu2tril reflects the upper triangular part of A into the lower triangular part of B, and tril2triu does the reverse. Given a single stored triangle these reconstruct the full symmetric matrix in one pass:

dtriu2tril( order, M, N, k, A, LDA, B, LDB )
dtriu2tril.ndarray( M, N, k, A, sa1, sa2, offsetA, B, sb1, sb2, offsetB )

Since A and B carry independent strides, the reflection is not a second pass over the output, writing with sb1 and sb2 swapped is what performs it.

Open PRs

These five complete the family across all five dtypes.

5. blas/ext/base/ndarray/*

Merged PRs

The ndarray-object form of the triangular family

Ndarray-object Implementation:

dtriu( [ A, B, k ] )
dtril( [ A, B, k ] )

A and B are two-dimensional ndarrays and k is a zero-dimensional one. M, N, LDA, LDB and order are all recovered from the ndarray metadata, so a mismatched leading dimension is not expressible here. k is a generic-dtype scalar rather than int8, since unlike a transpose flag its useful range depends on the matrix dimensions.

dtriu2tril( [ A, B, k ] )
dtril2triu( [ A, B, k ] )

6. Enum support for Level 2 and Level 3 blas/base/* routines

if ( !isMatrixTriangle( uplo ) ) {
    throw new TypeError( ... );
}

The new form resolves instead, accepting either a string or the matching enumerated integer and returning a canonical value or null:

u = resolveStr( uplo );
if ( u === null ) {
    throw new TypeError( ... );
}

Merged PRs

Current state

Level 1 is effectively finished in JavaScript across real single, double and generic precisions, and the complex routines have caught up with them; Every Level 2 real routine in the general, symmetric, packed and triangular families is now implemented in JavaScript at all three forms and for double, single and generic dtypes, The complex side has started rather than finished: the general matrix-vector product is merged for single precision and in review for double, while the Hermitian family and the banded and packed complex variants are untouched. The banded and packed-triangular routines are also still open on the real side.

At Level 3 the three real matrix-matrix multiplication packages are complete,The rest of the level, the symmetric and triangular matrix-matrix routines and everything complex has not been started.

Note: open and draft PRs that are on the finish line are counted as implementation here, alongside merged ones.

What remains

The tracking issue has moved forward this summer, but plenty of packages are still open. Taking it level by level, as #2039 currently stands:

Level 1 real and complex packages:

  • Nothing outstanding, the last gaps, crotg and zrotg, both have PRs in review

Level 2 real single and double precision packages:

  • Nothing outstanding, the banded, packed triangular and spr2 routines all have PRs in review

Level 2 complex single and double precision packages:

  • 10/17 packages remaining in each precision

Level 3 real single and double precision packages:

  • Nothing remaining, all packages have PRs in review

Level 3 complex single and double precision packages:

  • 8/9 packages remaining in each precision

Note: I have left out the WebAssembly implementations and packages whose blockers are not yet cleared, since that has to happen first.

Challenges and lessons learned

Working on different implementations exposed me to a wide range of challenges, and each one improved the way I approach development. Comparing my early work with my current process, the biggest change is that I now pay much more attention to preparation, consistency, and edge cases.

I learned that good implementation begins before writing code. Studying documentation and similar packages helped me understand project conventions, while careful testing showed me how easily bugs can remain hidden under common inputs. I also became more disciplined about reviewing my work and writing clear, meaningful commit messages.

Some tasks that initially appeared minor turned out to be important foundations for later work. This taught me to understand the wider purpose and dependencies of a task instead of judging it only by its size.

Whenever I faced difficulties, discussions with my mentors helped me find the right direction. Gradually, these lessons became part of my regular workflow and made me a more patient, systematic, and thoughtful developer.

Conclusion

I came into this able to use linear algebra libraries and I am leaving able to contribute in. That is a smaller sentence than it deserves: knowing that y = α·A·x + β·y is a matrix-vector product is a different thing from knowing why β = 0 needs its own branch, why the loop order in a rank-2 update is constrained, why a packed matrix cannot be addressed with two strides, or why passing an option as a string quietly forecloses your entire bindings roadmap. Every abstraction I had been standing on turned out to have a reason underneath it, and this summer was mostly spent finding out what those reasons were.

I want to thank @kgryte, whose reviews were exacting in the most useful possible way and who was reliably available at hours that cannot have suited his timezone. A great deal of what I now consider ordinary practice came from his comments. Thanks also to Gunj Joshi (@gunjjoshi) for mentoring me through this project, and for pointing me at the right prior art more times than I can count. I would also like to thank @nakulkrishnakumar, @prajjwalbajpai, @pratikbhagwat, and @sachinpangal for making this journey both insightful and enjoyable. I’ll really miss our regular standups, the conversations, the shared struggles, and the feeling that we were all learning and moving forward together.

Clone this wiki locally