GSoC 2026 ‐ Pratik Bhagwat

About me

Hey there! I'm Pratik Bhagwat, a final-year Electronics and Telecommunication Engineering student at N. K. Orchid College of Engineering & Technology, Maharashtra, India. I started out in web development, but over time I got drawn to the mathematical and numerical side of computing, which led me to stdlib in late 2025. I began with its statistics routines and later moved into the LAPACK part of the library, which became my project on linear equation solvers.

Project overview

Almost every piece of numerical software eventually runs into the same question: given a matrix A and a vector b, find the x that satisfies Ax = b. It hides inside curve fitting, physics simulations, computer graphics, and machine learning, and heavyweight libraries like NumPy, MATLAB, and R all lean on LAPACK to answer it quickly and accurately. My project set out to give stdlib that same ability, but in pure JavaScript.

stdlib already had a growing @stdlib/lapack/base/* namespace, yet the routines needed to actually assemble a working solver from end to end were missing. Rather than pick routines at random, I organized the work around three solver families: general dense systems via LU factorization, symmetric-indefinite systems via a packed Bunch-Kaufman factorization, and tridiagonal/banded systems. Each of these needs its own stack of factorization, solve, norm, and utility routines, so the real work was building those stacks piece by piece, and, as tends to happen, I got further with some of them than with others.

Project recap

A solver in LAPACK is never one function. Solving Ax = b means chaining a factorization, a triangular or structured solve, and a supporting cast of matrix norms, NaN checks, random generators, and permutation helpers. Over the summer this settled into a routine of its own, which broke down into roughly three stages:

Planning the order of work. Before writing anything, I figured out how the routines depended on each other so I wasn't building on top of missing pieces.

• For each of my three solver families I traced the dependency chain all the way down to the leaf routines.
• I then worked bottom-up: leaves first, then the routines that used them, and the drivers last. That way every PR I opened could be tested against real dependencies instead of stubs.

Porting a routine. Each routine went through the same shape of implementation.

• I started from the Netlib reference Fortran and made sure I understood what it actually did before translating it.
• I wrote it as two layers: a main API matching the usual LAPACK/BLAS signature, and an ndarray API taking explicit strides and offsets so it fits the rest of stdlib.
• For row-major support I avoided writing the algorithm twice by treating row-major as the transpose case, for example swapping DL/DU in dlagtm or flipping uplo in the symmetric-packed routines.
• Where LAPACK mutates a scalar in place (like the tolerance in dlagts), I passed it as a single-element typed array, the same trick the rest of @stdlib/lapack uses.

Getting the numbers exactly right. This is where most of my time actually went.

• I generated test fixtures by compiling Netlib LAPACK with gfortran and running it on the same inputs, then checked my output against it.
• I made a point of covering the annoying cases: NaN, Inf, signed zero, empty or degenerate matrices, and every parameter branch (each norm/uplo/diag combo, each job value in dlagts, the 1×1 vs 2×2 pivot paths in dsptrf).
• When something came out almost right, it was always arithmetic ordering, not a logic bug, for instance a reciprocal multiply where LAPACK does a divide, or a different summation order in a Frobenius norm. Matching LAPACK's exact order, and reusing helpers like dlassq and dlacpy, is what closed the gap.

To make this concrete, here's what a finished routine looks like in use. dlantr computes a norm of a triangular or trapezoidal matrix. A flat buffer, though, does not describe a matrix on its own: the very same values represent different matrices depending on whether they are read in row-major or column-major order, so order is part of the matrix's meaning, not a free reinterpretation of a fixed buffer. To take the norm of one matrix under both layouts, the caller therefore supplies two buffers, each laying that same matrix out for its own layout, and dlantr returns the same value for both:

var Float64Array = require( '@stdlib/array/float64' );
var dlantr = require( '@stdlib/lapack/base/dlantr' );

// The upper-triangular matrix we want the norm of:
//   [ 1 2 3 ]
//   [ 0 5 6 ]
//   [ 0 0 9 ]

// Stored in row-major order (consecutive values run along each row):
var rowMajor = new Float64Array([
    1.0, 2.0, 3.0,
    0.0, 5.0, 6.0,
    0.0, 0.0, 9.0
]);

// The same matrix stored in column-major order (consecutive values run down each column):
var columnMajor = new Float64Array([
    1.0, 0.0, 0.0,
    2.0, 5.0, 0.0,
    3.0, 6.0, 9.0
]);

// Workspace (only referenced for the infinity-norm):
var work = new Float64Array( 3 );

// order, norm, uplo, diag, M, N, A, LDA, work
var out = dlantr( 'row-major', 'fro', 'upper', 'non-unit', 3, 3, rowMajor, 3, work );
// out => ~12.49

out = dlantr( 'column-major', 'fro', 'upper', 'non-unit', 3, 3, columnMajor, 3, work );
// out => ~12.49

Because each buffer already lays the matrix out for its declared layout, both calls describe the identical upper-triangular matrix and return the same Frobenius norm. One implementation serves both layouts, rather than a single buffer silently changing meaning when order flips. Swapping 'fro' for 'max', 'one', or 'inf' picks a different norm through the same code path. That's the shape I was aiming for with every routine: one implementation, and you configure it with arguments.

Completed work

The core of the project is a set of routines under @stdlib/lapack/base/, grouped below by review status. A few contributions outside the solver scope follow at the end.

Merged

• lapack/base/dlaisnan (#12183)

  Tests whether two double-precision arguments are unequal x != y, which returns true whenever an input is NaN. It exists so a compiler cannot optimize a NaN check away, and it is the primitive that disnan is built on.

• lapack/base/disnan (#12289)

  Reports whether a double-precision value is NaN by calling dlaisnan(x, x). A small but foundational IEEE check relied on across the library.

• lapack/base/dlarf (#12331)

  Applies a single elementary (Householder) reflector H = I − tau·v·vᵀ to a general matrix from either the left or the right. Householder reflectors are the basic tool behind QR factorizations and many other orthogonal transformations.

• lapack/base/dlaruv (#12443)

  Generates a vector of up to 128 uniformly distributed random numbers on the interval (0, 1) using a multiplicative congruential method. It is the low-level generator that dlarnv drives.

Under review

Matrix norms. Each computes the one-norm, infinity-norm, Frobenius norm, or maximum absolute value of a matrix with a particular structure:

• lapack/base/dlangt (#12522): a general tridiagonal matrix supplied as its three diagonals.
• lapack/base/dlanst (#12957): a real symmetric tridiagonal matrix, stored as its main diagonal and off-diagonal.
• lapack/base/dlanhs (#13178): an upper Hessenberg matrix (upper triangular plus a single sub-diagonal).
• lapack/base/dlansb (#13182): a symmetric band matrix held in banded storage with a given number of super-diagonals.
• lapack/base/dlansy (#12837): a real symmetric matrix held in full storage.
• lapack/base/dlansp (#12640): a real symmetric matrix held in packed storage.

Random generation and permutations.

• lapack/base/dlarnv (#13107)

  Fills a vector with random numbers from a chosen distribution (uniform on (0, 1), uniform on (−1, 1), or standard normal), building on dlaruv.

• lapack/base/dlapmr (#12367)

  Rearranges the rows of a matrix according to a permutation vector, applied either forwards or backwards.

General dense (LU) and equilibration.

• lapack/base/dgetf2 (#12897)

  Computes the LU factorization A = P·L·U of a general m-by-n matrix using partial pivoting with row interchanges. This is the unblocked, level-2 BLAS variant of the factorization; the higher-level dgetrf computes the same result for larger matrices using a blocked algorithm.

• lapack/base/dgeequ (#13155)

  Computes row and column scaling factors that equilibrate a general matrix to reduce its condition number.

• lapack/base/dgeequb (#13205)

  Computes the same kind of scaling factors as dgeequ, but restricted to powers of two, so the scaling itself introduces no rounding error.

Symmetric packed (Bunch-Kaufman).

• lapack/base/dsptrf (#13288)

  Factorizes a real symmetric indefinite matrix held in packed storage using the Bunch-Kaufman diagonal pivoting method, producing A = U·D·Uᵀ (or L·D·Lᵀ) with a mix of 1×1 and 2×2 pivot blocks.

• lapack/base/dlaqsp (#12950)

  Applies precomputed scaling factors to equilibrate a symmetric matrix in packed storage, reporting whether equilibration was actually performed.

Tridiagonal.

• lapack/base/dgtsv (#13222)

  Solves the system A·X = B for one or more right-hand sides, where A is a general tridiagonal matrix, using Gaussian elimination with partial pivoting. This is the user-facing driver of the tridiagonal family.

• lapack/base/dlagtf (#13140)

  Factorizes T − λI (a tridiagonal matrix shifted by a scalar) as P·L·U with partial pivoting, producing the factors that dlagts consumes.

• lapack/base/dlagts (#13800)

  Solves (T − λI)x = y, or its transpose, from the factorization produced by dlagtf, across four job modes that also control overflow-avoiding perturbation.

• lapack/base/dlagtm (#13012)

  Computes the matrix-matrix product B := α·op(A)·X + β·B, where A is tridiagonal and op(A) is either A or Aᵀ.

• lapack/base/dptts2 (#13703)

  Solves a symmetric positive-definite tridiagonal system A·X = B using the L·D·Lᵀ factorization produced by dpttrf. It is the computational kernel behind dpttrs, and one of the routines I matched bit-for-bit against Netlib.

Ready to submit

• lapack/base/dlantr

  Computes the one-norm, infinity-norm, Frobenius norm, or maximum absolute value of a triangular or trapezoidal matrix, handling both memory layouts and all four norm types through a single implementation. Finished and verified locally against Netlib; the pull request is ready to open.

Work outside the solver scope

Several contributions outside the core solver work, in stdlib's statistics namespace:

• stats/base/ndarray/dvariancech (#14090)
• stats/base/ndarray/nanvariance (#14359)
• stats/base/ndarray/dnanvariance (#14386)
• stats/base/ndarray/nanvariancech (#14540)
• stats/base/ndarray/nanvariancepn (#14541)
• stats/base/ndarray/nanvariancetk (#14559)

Current state

The work above sits in one of four states. Four routines are already merged and available in stdlib: dlaisnan, disnan, dlarf, and dlaruv. Eighteen more are under review (the matrix norms, the random and permutation helpers, the LU kernel and both equilibration routines, the symmetric-packed factorization and its scaling helper, and the full tridiagonal stack, all listed above), meaning their implementations and tests are complete and the pull requests are open but not yet merged. One routine, dlantr (the norm of a triangular or trapezoidal matrix), is finished and verified locally but has no pull request opened yet. Everything past that is still remaining, meaning not yet implemented, and is spelled out in the next section.

Of the three solver families, the tridiagonal family is the furthest along: its simple driver dgtsv is written and in review, together with its supporting norms, factorization, and solve routines. The LU and symmetric-packed families have their factorization kernels in place but not yet their top-level drivers. I'm still iterating on review feedback for the open pull requests.

What remains

LAPACK is an enormous library, with roughly 1,700 routines, so completing every aspect of linear-equation solving was never the aim. Instead, I focused on building complete solver pipelines across three families, and what is left in each case is the short chain of routines that turns the factorization kernels I have already written into a full A·X = B solve.

• Which routines remain for each solver family.

  • General dense (LU). The simple driver dgesv factorizes with dgetrf and then solves with dgetrs. I already have the unblocked kernel dgetf2 written and in review, and stdlib provides dlaswp and dgemm, so what remains is the level-3 BLAS triangular solve dtrsm (which stdlib does not yet have), the blocked factorization dgetrf, the solve dgetrs, and the driver dgesv itself.
  • Symmetric-indefinite packed (Bunch-Kaufman). The simple driver dspsv pairs the factorization dsptrf with the solve dsptrs. The factorization is written and in review, as is the scaling helper dlaqsp, and the solve depends only on BLAS routines stdlib already provides, so what remains is dsptrs and the driver dspsv.
  • Tridiagonal. The general branch is essentially done: the driver dgtsv is written and in review, along with dlangt, dlagtm, dlagtf, and dlagts, so that pipeline is complete once merged. For positive-definite tridiagonal systems, the plan reuses the existing dpttrf factorization and adds a solve, dpttrs, whose computational kernel dptts2 is written and in review; what is left there is the thin dpttrs wrapper on top of it.

• Which completed routines are awaiting review. The eighteen routines enumerated under Completed work (the matrix norms, the random and permutation helpers, the LU kernel and equilibration routines, the symmetric-packed factorization and its scaling helper, and the tridiagonal stack) are implemented and tested but not yet merged. Carrying them through review is itself a meaningful part of finishing the pipelines, since every merged routine becomes a verified dependency the drivers can build on rather than a moving target.

• Which top-level drivers are still missing. The proposal targets three drivers: dgesv (general dense), dspsv (symmetric packed), and dgtsv (general tridiagonal). Two are still unwritten, dgesv and dspsv, while dgtsv is written and in review. The positive-definite tridiagonal family has no separate driver in the plan; it becomes usable as soon as the dpttrs solve lands alongside the dpttrf factorization already in stdlib. Each remaining driver is thin on its own, so the real effort in each case is finishing the factorization and solve routines beneath it first.

• Whether any work was descoped from the original proposal. The proposal targeted 42 double-precision routines forming three full solver pipelines, complete with equilibration and iterative refinement, plus a banded-LU branch. Only one routine was deliberately removed from scope: dlabad, whose pull request I closed because stdlib assumes IEEE 754 compliance, which makes the routine unnecessary. Nothing else was cut, but a few groups were deferred as later follow-ups rather than finished within the twelve weeks: the iterative-refinement routines (dgerfs, dsprfs, dgtrfs, dptrfs); the dense, packed, and banded triangular solvers, including their scaled variants (dtrtrs, dtptrs, dtbtrs, dlatbs, dlatps, dlaln2); the banded-LU pipeline (dgbtf2, dgbtrs); and the standalone tridiagonal solve dgttrs, whose one-shot equivalent dgtsv is already done. The priority throughout was getting the simple-driver solve path of each family correct and idiomatic first.

• A tracking issue. The remaining routines, along with the ones already merged or under review, are tracked in issue, which serves as the running checklist for finishing these solver pipelines.

Beyond wiring up the drivers, I want to explore C implementations to improve performance and eventually expand into areas such as eigenvalue computations and QR algorithms.

Measuring against the plan

Here's how the finished work lines up with what I actually proposed.

What the proposal expected to complete. The proposal committed to 42 double-precision LAPACK routines over a twelve-week schedule, organized as three complete solver pipelines: a general dense (LU) pipeline ending in the dgesv driver, a symmetric-packed (Bunch-Kaufman) pipeline ending in dspsv, and a tridiagonal pipeline ending in dgtsv. Each pipeline was to be fleshed out with equilibration, iterative refinement, and the relevant triangular solvers, and the plan also included a banded-LU branch (dgbtf2, dgbtrs). The approach was bottom-up at roughly three to four routines a week, so that leaf routines would clear review while the drivers above them were being written.

What was actually completed. Measured against that plan, twenty-two of the forty-two proposed routines are finished, together with dptts2 (the solve kernel I added beneath dpttrs), twenty-three in all, comprising the four merged, eighteen under review, and one local summarized under Current state above. Together they cover the leaf and mid-level layers of all three families: the NaN and random-number primitives, the full set of matrix norms, the LU kernel dgetf2 and both equilibration routines, the symmetric-packed factorization dsptrf with its scaling helper, and the entire tridiagonal stack up to and including the dgtsv driver. What did not land within the timeline were the upper layers: the remaining drivers dgesv and dspsv, the solves dgetrs, dsptrs, and dpttrs, the iterative-refinement routines, and the banded branch, all itemized in What remains above.

Why the timeline and scope changed. The main factor was simply that a single LAPACK routine in stdlib is a full package rather than one function. The dual-layout APIs, tests, documentation, and reference-generated fixtures around the algorithm (detailed under Challenges and lessons learned) take far more time than the porting. Carrying one routine through all of that and then through several rounds of review routinely stretched well beyond writing the code, so real throughput sat below the planned three-to-four routines a week. Rather than trade that quality away to hit a routine count, I chose to finish each family from the leaves upward and take the tridiagonal pipeline all the way to its driver, which pushed the remaining drivers and the expert-level extras past the end of the program.

What I learned about estimating LAPACK work. If I were planning this again, I'd count in reviewed-and-merged routines, not written ones. The algorithm is usually the smallest piece, and how fast things actually move depends far more on the surrounding package work and on review turnaround than on the porting. I'd also weight the schedule by dependency depth. Leaf routines are cheap and easy to run in parallel, but each driver quietly carries the cost of everything under it, so a plan that looks linear on paper is really back-loaded toward the final weeks. An honest version of this proposal would have committed to fewer routines with deeper, driver-complete pipelines, and treated forty-two as a stretch goal rather than a baseline.

Challenges and lessons learned

The thing that surprised me most was how much goes into a single LAPACK routine before it's actually done. The algorithm is honestly the easy part. On top of it there are two APIs to support, both row- and column-major layouts to handle, benchmarks to write, TypeScript declarations to keep in sync, the exact norm equations to document in the README, and fixtures that have to hit all the edge cases. I underestimated all of that early on. Working bottom-up along the dependency chains helped here, and it's a big part of why the tridiagonal stack came together as well as it did.

Numerical work also turned out to be much less forgiving than I expected. I'd write a routine, watch it match the reference to eleven decimal places, and assume I was finished, only to find it was off in the very last bit. It was never a real bug in the logic. It'd be something small, like the order of a subtraction, or a spot where I'd written a reciprocal multiply and LAPACK does a plain divide. Chasing those down is what drove home that "looks right" and "is right" aren't the same thing. Widening a tolerance would have hidden it; the actual fix is to match the reference's arithmetic exactly. Column-major versus row-major storage, and LAPACK's 1-based, sign-encoded index conventions, took a while to get used to too, but after enough routines it stopped being something I had to think about.

In a codebase as convention-heavy as stdlib, most of the back-and-forth wasn't about the math, since the fixtures already covered that. It was about fitting in: the shared validation patterns, indexing the way the rest of the library does it. Once I started keeping review comments as a checklist and applying them to the next routine before opening the PR, things moved a lot faster. It also changed how I think about progress. A finished, tested routine sitting in review isn't "not done yet", it's real work making its way through the pipeline.

Conclusion

What I set out to do was give stdlib a real path from a matrix and a right-hand side to a solution, all in plain JavaScript, and a good chunk of that foundation is now in place. Over the summer I got a piece of all three solver families working: the NaN and random-generation primitives, the matrix-norm routines, the LU factorization kernel and equilibration, the symmetric-packed Bunch-Kaufman factorization, and, furthest along, the whole tridiagonal stack up to its driver. Some of these match the Netlib reference bit-for-bit; the rest I checked against reference-generated fixtures with appropriate numerical comparisons. The part I'll probably get the most use out of, though, is what came out of the process itself: a fixture-generation workflow I can reuse, and a much clearer sense of what stdlib expects from a numerical routine.

I owe a lot of this to the people around me this summer. Thank you to the org admins, ⁠Athan Reines and ⁠Philipp Burckhardt, for running such a welcoming and well-organized program. A special thank you to my mentor, ⁠Karan Anand, whose reviews were detailed, patient, and consistently educational. He put real effort into helping me understand not just what to change but why, and a huge part of how much I grew this summer traces back to that guidance. I’m also grateful to my fellow contributors, ⁠Kaustubh Patange, ⁠Nakul Krishnakumar, and ⁠Sachin Pangal for the shared momentum and encouragement throughout. I also want to thank ⁠Prajjwal Bajpai, whose project also focused on LAPACK, for his help throughout the work, especially in understanding and generating fixtures. His guidance was valuable in helping me better understand the validation process for numerical routines. This was my first time contributing to something this large and this careful, and it has raised the bar for how I think about correctness and collaboration. I plan to continue working on the LAPACK routines and contributing to stdlib after the program ends as well.

Thanks everyone for reading! Cheers!