-
-
Notifications
You must be signed in to change notification settings - Fork 1.3k
GSoC 2026 ‐ Prajjwal Bajpai
Hi there! I am Prajjwal Bajpai from Greater Noida, Uttar Pradesh, India. I am a final year student at the Indian Institute of Technology(BHU), Varanasi, Uttar Pradesh, India. I have been interested in mathematics since school and have been coding for more than three years, and last year, I began my open-source journey with stdlib, where I have learned a great deal. Apart from software, I am also interested in cars and motorsports.
My project aimed to do the JavaScript implementation of the dgesvd (DGESVD computes the singular value decomposition (SVD)
for general matrices) routine under Singular Value Decomposition (SVD): Standard SVD driver, A = UΣV^H. The implementation was carried out by starting with the low-level routines and progressing upward through the dependency tree to the high-level, user-facing routines.
The implementation of a routine generally consisted of the following steps:
- Base implementation
-
ndarrayAPI(stdlib) - Support for both row-major(C-style) and column-major(Fortran-style) memory layouts
- Testing the implementation against Fortran output or through mathematical verification
- Benchmarking the routine
- Writing the documentations for the routine.
For the base implementation of a routine, I primarily referred to the Netlib's LAPACK documentation and LAPACKE C Interface to LAPACK. LAPACKE provides a C wrapper around the original Fortran routines while supporting both column-major (Fortran-style) and row-major (C-style) memory layouts.
In addition, I made updates to some existing LAPACK/BLAS routine PRs to ensure that they met the required standards.
Prior to GSoC, I built an LAPACK dependency graph. The dependency graph on the Netlib site is not always accurate and it includes BLAS routines. To tackle this, I parsed each file and checked for any valid LAPACK routine name followed by ‘(‘ which indicates that routine is being called and is therefore a dependency. This approached fixed the issue of encountering cycles during the topological sort of routines and produced a valid working order. All the graphs and orders are stored in this GitHub repository.
TO elaborate on the aforementioned process of implementing a LAPACK routine,the base API required strides and offsets for 1-D and 2-D arrays in addition to the routine's standard input parameters. Scalar output values were preferably stored in an out array to keep the implementation clean and convenient. This base implementation was kept private, with <routine_name>.js and ndarray.js providing the public-facing APIs on top of it.
The sources I used for reference primarily support column-major storage for 2-D matrices. Although interfaces such as LAPACKE provide an option for row-major storage, they internally transform the input to column-major format before calling the underlying LAPACK routine and then transform the output back to row-major format. Stdlib provides native support for both row-major and column-major storage layouts for LAPACK routines, providing greater flexibility for users and avoiding these additional transformations.
To support both row-major and column-major layouts, matrices are represented in linear memory using strides and offsets. For example, a matrix A:
is represented in row-major format as,
A = [ 1, 2, 3, 4, 5, 6, 7, 8, 9 ] // (`strideA1 = 3`, `strideA2 = 1`, `offsetA = 0`)and in column-major format as,
A = [ 1, 4, 7, 2, 5, 8, 3, 6, 9 ] // (`strideA1 = 1`, `strideA2 = 3`, `offsetA = 0`)To expose these strides and offsets in the ndarray API, the function signature had to be adapted from the original Fortran implementation. For example, the Fortran dlarft routine has the following signature:
SUBROUTINE DLARFT ( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )The corresponding base and ndarray API in JavaScript uses the following signature:
function dlarft( direct, storev, N, K, V, strideV1, strideV2, offsetV, TAU, strideTAU, offsetTAU, T, strideT1, strideT2, offsetT )whereas the dlarft.js API more closely matches the original Fortran signature:
function dlarft( order, direct, storev, N, K, V, LDV, TAU, T, LDT )To further optimize these routines, techniques such as loop reordering and loop tiling were employed, with the goal of maintaining comparable performance for both column-major and row-major arrays.
For testing, I compared the outputs of our implementation against those of the corresponding LAPACK routine across various test cases. The outputs were stored in JSON format and used as fixtures for tests written with tape. Achieving 100% test coverage for the higher-level routines required creating 128 × 128 2-D arrays to trigger the blocked algorithms. A blocked algorithm is a specialized execution path that processes data in submatrices (blocks) rather than element by element, improving cache reuse and overall computational efficiency.
However, creating multiple JSON fixtures of this size would have resulted in unnecessarily large test files. Therefore, we adopted mathematical verification techniques, such as checking the orthogonality of the resulting matrices, to test these code paths without relying on large output fixtures.
LAPACK Routines
-
Performs a single dqd (differential quotient-difference) transformation with a zero shift in the dqds algorithm, using a ping-pong data representation for improved computational efficiency. It incorporates safeguards against underflow and overflow while updating intermediate quantities required for accurate singular value computations of bidiagonal matrices.4
-
Performs a single dqds (differential quotient-difference with shifts) transformation in ping-pong form as part of LAPACK's singular value computation algorithm for bidiagonal matrices. It applies a specified shift while updating the qd array and tracking key intermediate quantities, with separate execution paths for IEEE and non-IEEE arithmetic to ensure numerical stability and robustness.
-
Computes an approximation to the smallest eigenvalue of a bidiagonal matrix by estimating an appropriate shift (
TAU) from the results of the previous dqds transformation. It employs multiple shift-selection strategies based on the current state of the computation, improving both the convergence rate and numerical stability of the dqds algorithm used in singular value decomposition. -
Sorts the elements of a real-valued array in either increasing or decreasing order, as specified by the input parameter. It uses a hybrid sorting strategy that combines Quick Sort for large partitions with Insertion Sort for small subarrays, providing efficient performance while maintaining low overhead for small datasets.
-
Constructs an
MXNreal orthogonal matrixQwith orthonormal rows from the elementary Householder reflectors produced by the LQ factorization routine DGELQF. The matrix is formed asQ=H(k)H(k−1)⋯H(2)H(1), where eachH(i)is an elementary reflector. DORGL2 implements the unblocked (Level-2 BLAS) algorithm, making it suitable for smaller matrices or as the base case of the blocked routine DORGLQ. -
Applies the orthogonal matrix
Q, generated during the LQ factorization by DGELQF, to a general matrix C. Depending on the input parameters, it computes one ofQ*C,Q^T*C,C*Q, orC*Q^T, whereQ=H(k)H(k−1)⋯H(2)H(1), and eachH(i)is an elementary Householder reflector. DORML2 implements the unblocked (Level-2 BLAS) algorithm, making it suitable for smaller matrices and serving as the base case for the blocked routine DORMLQ. -
Applies a sequence of Givens (plane) rotations to a real matrix
A, either from the left or the right. Depending on theSIDEparameter, it performsA←PAorA←AP^T, wherePis an orthogonal matrix implicitly represented as a sequence of plane rotations. The routine supports variable, top, and bottom pivoting strategies, as well as forward and backward application orders, enabling efficient orthogonal transformations without explicitly formingP. -
Computes the LQ factorization of a real m×n matrix A using an unblocked (Level-2 BLAS) algorithm, decomposing it as
A=LQ, whereLis a lower triangular (or lower trapezoidal) matrix andQis an orthogonal matrix represented implicitly as a product of Householder reflectors,Q=H(k)H(k−1)⋯H(2)H(1),k=min(M,N), with each reflector defined byH(i)=I−τi*vi*vi^T. The routine stores the reflectors compactly in the input matrix and the scalar factors in the TAU array, allowing Q to be formed or applied later by routines such as DORGL2 and DORML2. -
Computes the QR factorization of a real m×n matrix A using an unblocked (Level-2 BLAS) algorithm, decomposing it as
A=QR, whereQis an orthogonal matrix andRis an upper triangular (or upper trapezoidal) matrix. The orthogonal factor is represented implicitly as a product of Householder reflectors,Q=H(1)H(2)⋯H(k),k=min(M,N), with each reflector defined byH(i)=I−τi*vi*vi^T. The reflectors are stored compactly in the lower triangular part of the input matrix, while their scalar factors are stored in theTAUarray, allowingQto be efficiently formed or applied by subsequent LAPACK routines. -
Constructs the triangular factor
Tassociated with a block Householder reflector formed from k elementary reflectors. Together with the reflector matrixV, it provides the compact WY representationH=I−V*T*V^T(orH=I−V^T*T*Vwhen the reflectors are stored row-wise), enabling multiple Householder reflectors to be applied as a single block transformation. This compact representation is fundamental to LAPACK's blocked QR, LQ, RQ, and related factorizations, significantly improving performance by leveraging Level-3 BLAS operations. -
Reduces the first nb rows and columns of a real m×n matrix A to bidiagonal form by applying orthogonal transformations,
B=Q^T*A*P, whereBis upper bidiagonal whenM≥Nand lower bidiagonal whenM<N. In addition to storing the Householder reflectors definingQandP, the routine computes the auxiliary matricesXandY, which enable efficient blocked updates of the remaining unreduced submatrix in the higher-level DGEBRD algorithm using Level-3 BLAS operations. -
Reduces a real m×n matrix
Ato bidiagonal form using an unblocked (Level-2 BLAS) algorithm by computing the orthogonal transformationB=Q^T*A*P, whereBis upper bidiagonal ifM≥Nand lower bidiagonal otherwise. The orthogonal matricesQandPare represented implicitly as products of Householder reflectors, whose scalar factors are stored in theTAUQandTAUParrays. As the unblocked counterpart of DGEBRD, DGEBD2 is primarily used for smaller matrices or as the base case within the blocked bidiagonal reduction algorithm. -
Performs a single iteration of the dqds (differential quotient-difference with shifts) algorithm by checking for deflation, selecting an appropriate shift
τ, and invoking the dqds transformation. If the chosen shift produces non-positive intermediate values, the routine adaptively modifies the shift and retries until a numerically stable update is obtained, ensuring robust and accurate convergence during the computation of singular values of bidiagonal matrices. -
Reduces a real
M×NmatrixAto bidiagonal form using a blocked (Level-3 BLAS) algorithm, computing the orthogonal transformationB=Q^T*A*P, whereBis upper bidiagonal ifM≥Nand lower bidiagonal otherwise. The orthogonal matricesQandPare represented implicitly as products of Householder reflectors, while the routine employs DLABRD to reduce matrix panels and DGEBD2 for the remaining unreduced portion, enabling high performance through cache-efficient Level-3 BLAS updates. -
Computes all the eigenvalues of a symmetric positive definite tridiagonal matrix associated with a qd array
Zto high relative accuracy using the dqds (differential quotient-difference with shifts) algorithm. IfLandUare the unit lower and upper bidiagonal matrices encoded byZ, the routine computes the eigenvalues ofT=L*U, without explicitly formingT. It repeatedly invokes DLASQ3 for deflation and shift selection, and returns the computed eigenvalues in sorted order together with diagnostic information such as the iteration count and shift statistics. -
Generates an m×n real orthogonal matrix
Qwith orthonormal columns from the Householder reflectors produced by the QR factorization routine DGEQRF. The matrix is constructed asQ=H(1)H(2)⋯H(k), where eachH(i)is an elementary Householder reflector. DORGQR employs a blocked (Level-3 BLAS) algorithm that uses DLARFT to construct the triangular factor of block reflectors and DLARFB to apply them efficiently, while falling back to DORG2R for the unblocked portion of the computation. -
Generates an m×n real orthogonal matrix
Qwith orthonormal rows from the Householder reflectors produced by the LQ factorization routine DGELQF. The matrix is constructed asQ=H(k)H(k−1)⋯H(2)H(1), where eachH(i)is an elementary Householder reflector. DORGLQ employs a blocked (Level-3 BLAS) algorithm that uses DLARFT to construct the triangular factor of block reflectors and DLARFB to apply them efficiently, while falling back to DORGL2 for the unblocked portion of the computation. -
Applies the orthogonal matrix
Q, generated during the QR factorization by DGEQRF, to a general matrixC. Depending on the input parameters, it computesQ*C,Q^T*C,C*Q, orC*Q^T, whereQ=H(1)H(2)⋯H(k), and eachH(i)is an elementary Householder reflector. DORMQR employs a blocked (Level-3 BLAS) algorithm that uses DLARFT and DLARFB to efficiently apply blocks of reflectors, while reverting to DORM2R for the unblocked portion of the computation. -
Applies the orthogonal matrix
Q, generated during the LQ factorization by DGELQF, to a general matrixC. Depending on the input parameters, it computesQ*C,Q^T*C,C*Q, orC*Q^T, whereQ=H(k)H(k−1)⋯H(2)H(1), and eachH(i)is an elementary Householder reflector. DORMLQ employs a blocked (Level-3 BLAS) algorithm that uses DLARFT and DLARFB to efficiently apply blocks of reflectors, while reverting to DORML2 for the unblocked portion of the computation. -
Computes the singular values of a real
N×Nbidiagonal matrixBwith high relative accuracy by applying the dqds (differential quotient-difference with shifts) algorithm. Given the diagonal and off-diagonal elements ofB, it computes the singular valuesσi = λi*(B^T*B), without explicitly formingB^T*B, thereby improving numerical stability and avoiding unnecessary computation. Internally, the routine scales the input to prevent overflow and underflow, invokes DLASQ2 to perform the dqds iterations, and returns the singular values in decreasing order. -
Computes the QR factorization of a real
M×NmatrixAusing a blocked (Level-3 BLAS) algorithm, decomposing it asA=Q*R, whereQis an orthogonal matrix andRis an upper triangular (or upper trapezoidal) matrix. The orthogonal factor is represented implicitly as a product of Householder reflectors,Q=H(1)H(2)⋯H(k),k=min(M,N), with each reflector defined byH(i)=I−τi*vi*vi^T. The routine factors the matrix panel-by-panel using DGEQR2, constructs block reflectors with DLARFT, and applies them using DLARFB, enabling cache-efficient execution through Level-3 BLAS operations. -
Computes the LQ factorization of a real
M×NmatrixAusing a blocked (Level-3 BLAS) algorithm, decomposing it asA=L*Q, whereLis a lower triangular (or lower trapezoidal) matrix andQis an orthogonal matrix. The orthogonal factor is represented implicitly as a product of Householder reflectors,Q=H(k)H(k−1)⋯H(2)H(1),k=min(M,N), with each reflector defined byH(i)=I−τi*vi*vi^T. The routine factors the matrix panel-by-panel using DGELQ2, constructs block reflectors with DLARFT, and applies them using DLARFB, enabling cache-efficient execution through Level-3 BLAS operations. -
Generates one of the orthogonal matrices
QorP^Tproduced during the bidiagonal reduction performed by DGEBRD, whereA=Q*B*P^T, andBis a bidiagonal matrix. Depending on theVECTparameter, it constructs eitherQ=H(1)H(2)⋯H(k)orP^T=G(k)⋯G(2)G(1), from the stored Householder reflectors. The routine employs the blocked routines DORGQR or DORGLQ to efficiently generate the requested orthogonal matrix using Level-3 BLAS operations. -
Applies one of the orthogonal matrices
QorP, generated during the bidiagonal reduction performed by DGEBRD, to a general matrixC. Depending on the input parameters, it computesQ*C,Q^T*C,C*Q,C*Q^T,P*C,P^T*C,C*P, orC*P^T, where the bidiagonal reduction satisfiesA=Q*B*P^T, withQandPrepresented implicitly as products of Householder reflectors. The routine internally invokes DORMQR or DORMLQ to efficiently apply the requested orthogonal transformation using blocked (Level-3 BLAS) operations.
Cleanup and maintenance
I am still working on the dbdsqr routine and cleaning up some of the existing PRs. Apart from these, the JavaScript implementations of all routines that are listed above are complete.
Although I was not able to reach the stretch goal of implementing the dgesvdq routine mentioned in my proposal, I have nearly reached the primary target of this project, i.e., the dgesvd routine. Two routines remain: dgesvd itself and dbdsqr.
I will be working on completing them post GSoC as well. However, my current priority is to tidy up the existing GSoC PRs and get them reviewed and merged.
There were a number of challenges that I faced during this project, many of them unexpected, but each one taught me something valuable.
While implementing the dgelq2 routine, the column-major outputs matched those of the Fortran implementation, but the row-major outputs did not. After extensive investigation, along with guidance from my mentor, I realized that the issue was caused by differences in the order of multiplication resulting from optimizations for row-major and column-major layouts in the dger BLAS routine. These very small floating-point errors accumulated over multiple passes and eventually led to noticeable differences in the output. This experience taught me the importance of patience when debugging complex numerical issues and of seeking help when needed.
This investigation also led me to discover mathematical verification as an alternative testing method, which could potentially be applied to other LAPACK routines in the future.
Another challenge I encountered was with ULP-based testing in higher-level routines. ULP-based comparisons can be unreliable for values close to zero, as accumulated floating-point errors may change the sign of a value, resulting in a very large ULP difference despite the values being numerically close. With guidance from my mentors, I developed a method to account for this issue.
Apart from these challenges, I frequently introduced minor issues in my PRs that were identified during the review process. While these mistakes were sometimes frustrating, the review process helped me improve significantly.
Overall, I learned more during my time consistently contributing to this project than I had through any previous development experience. The mistakes I made and the issues I encountered helped me grow as a developer and gain a deeper understanding of numerical computing, low-level implementation details, and the mathematical foundations behind these routines.
Contributing to stdlib through Google Summer of Code 2026 has been one of the most valuable experiences of my journey as a software developer. This project gave me the opportunity to work extensively with numerical algorithms and understand how mathematical routines are translated into efficient, production-quality software. Along the way, I gained a much deeper understanding of floating-point arithmetic, memory layouts, performance optimization, and numerical testing. The challenges I encountered throughout the project also helped me improve my problem-solving skills, become more patient with complex debugging, and develop a stronger understanding of the low-level details involved in numerical computing.
I would like to sincerely thank my mentor, Aayush Khanna, for his constant guidance and support throughout the project. I would also like to thank Athan Reines, Philipp Burckhardt, Mara Averick, Karan Anand, and Gunj Joshi for their guidance and support. I am also grateful to Pratik Bhagwat for the collaboration throughout the project as we worked on LAPACK-related GSoC projects together, and to other fellow GSoC 2026 contributors Kaustubh Patange, Nakul Krishnakumar, and Sachin Pangal. I look forward to continuing my contributions to stdlib beyond GSoC and exploring numerical computing and high-performance scientific software further.