Efficiently solving nonlinear equations underpins numerous scientific and engineering disciplines, yet scaling these solutions for challenging system models remains a challenge. This article presents NonlinearSolve.jl —a suite of high-performance open source nonlinear equation solvers implemented natively in the Julia programming language. NonlinearSolve.jl distinguishes itself by offering a unified API that accommodates a diverse range of solver specifications alongside features such as automatic algorithm selection based on runtime analysis, support for static array kernels for improved GPU computation on smaller problems, and the utilization of sparse automatic differentiation and Jacobian-free Krylov methods for large-scale problem-solving. Through rigorous comparison with established tools such as PETSc SNES , Sundials KINSOL , and MINPACK , NonlinearSolve.jl demonstrates robustness and efficiency, achieving significant advancements in solving nonlinear equations while being implemented in a high-level programming language. The capabilities of NonlinearSolve.jl unlock new potentials in modeling and simulation across various domains, making it a valuable addition to the computational toolkit of researchers and practitioners alike.

This paper presents a portable, GPU-accelerated implementation of a QR-based singular value computation algorithm in Julia. The singular value decomposition (SVD) is a fundamental numerical tool in scientific computing and machine learning, providing optimal low-rank matrix approximations. Its importance has increased even more in large-scale machine learning pipelines, including large language models (LLMs), where it enables low-rank adaptation (LoRA). The implemented algorithm is based on the classic two-stage QR reduction, consisting of successive matrix reduction to band form and bidiagonal form. Our implementation leverages Julia’s multiple dispatch and metaprogramming capabilities, integrating with the GPUArrays and KernelAbstractions frameworks to provide a unified type and hardware-agnostic function. It supports diverse GPU architectures and data types, and is, to our knowledge, the first GPU-accelerated singular value implementation to support Apple Metal GPUs and half precision. Performance results on multiple GPU backends and data types demonstrate that portability does not require sacrificing performance: the unified function outperforms most linear algebra libraries (MAGMA, SLATE, rocSOLVER, oneMKL) for matrix sizes larger than 1024 × 1024, and achieves 80%-90% of the performance of cuSOLVER for large matrices.

Deep generative models have recently been applied to physical systems governed by partial differential equations (PDEs), offering scalable simulation and uncertainty-aware inference. However, enforcing physical constraints, such as conservation laws (linear and nonlinear) and physical consistencies, remains challenging. Existing methods often rely on soft penalties or architectural biases that fail to guarantee hard constraints. In this work, we propose Physics-Constrained Flow Matching (PCFM), a zero-shot inference framework that enforces arbitrary nonlinear constraints in pretrained flow-based generative models. PCFM continuously guides the sampling process through physics-based corrections applied to intermediate solution states, while remaining aligned with the learned flow and satisfying physical constraints. Empirically, PCFM outperforms both unconstrained and constrained baselines on a range of PDEs, including those with shocks, discontinuities, and sharp features, while ensuring exact constraint satisfaction at the final solution. Our method provides a flexible framework for enforcing hard constraints in both scientific and general-purpose generative models, especially in applications where constraint satisfaction is essential.

Neural Differential Equations have become an important modeling framework due to their ability to adapt to new problems automatically. Training a neural differential equation is effectively a search over a space of plausible dynamical systems. Controlling the computational cost for these models is difficult since it relies on the number of steps the adaptive solver takes. Most prior works have used higher-order methods to reduce prediction timings while greatly increasing training time or reducing both training and prediction timings by relying on specific training algorithms, which are harder to use as a drop-in replacement. In this manuscript, *we use internal cost heuristics of adaptive differential equation solvers at stochastic time-points to guide the training towards learning a dynamical system that is easier to integrate*. We ``close the blackbox'' and allow the use of our method with any sensitivity method. We perform experimental studies to compare our method to global regularization to show that we attain similar performance numbers without compromising on the flexibility of implementation. *We develop two sampling strategies to trade-off between performance and training time*. Our method reduces the number of function evaluations to 0.556x - 0.733x and accelerates predictions by 1.3x - 2x.

We propose a more open, efficient, expressive, and ergonomic future for AI development, machine learning, and data science based on the Julia programming language. Our thesis is that the current tapestry of high level codes with library calls creates programmer indirections that can work well for the “one off”, but can slow general progress. We provide examples from Machine Learning, Automatic Differentiation, and Data Handling Technologies.

Summary form only given. The Julia programming language is gaining enormous popularity. Julia was designed to be easy and fast. Most importantly, Julia shatters deeply established notions widely held in the applied community. Julia shows the fascinating dance between specialization and abstraction. Specialization allows for custom treatment. We can pick just the right algorithm for the right circumstance and this can happen at runtime based on argument types (code selection via multiple dispatch). Abstraction recognizes what remains the same after differences are stripped away and ignored as irrelevant. The recognition of abstraction allows for code reuse (generic programming). A simple idea that yields incredible power. Julia is many things to many people. In this talk we describe how Julia was built on the heels of our parallel computing experience with Star-P which began as an MIT research project and was a software product of Interactive Supercomputing. Our experience taught us that bolting parallelism onto an existing language that was not designed for performance or parallelism is difficult at best, and impossible at worst. One of our (not so secret) motivations to build Julia was to have the language we wanted for parallel numerical computing.

Approximating ideal program outputs is a common technique for solving computationally difficult problems, for adhering to processing or timing constraints, and for performance optimization in situations where perfect precision is not necessary. To this end, programmers often use approximation algorithms, iterative methods, data resampling, and other heuristics. However, programming such variable accuracy algorithms presents difficult challenges since the optimal algorithms and parameters may change with different accuracy requirements and usage environments. This problem is further compounded when multiple variable accuracy algorithms are nested together due to the complex way that accuracy requirements can propagate across algorithms and because of the size of the set of allowable compositions. As a result, programmers often deal with this issue in an ad-hoc manner that can sometimes violate sound programming practices such as maintaining library abstractions. In this paper, we propose language extensions that expose trade-offs between time and accuracy to the compiler. The compiler performs fully automatic compile-time and installtime autotuning and analyses in order to construct optimized algorithms to achieve any given target accuracy. We present novel compiler techniques and a structured genetic tuning algorithm to search the space of candidate algorithms and accuracies in the presence of recursion and sub-calls to other variable accuracy code. These techniques benefit both the library writer, by providing an easy way to describe and search the parameter and algorithmic choice space, and the library user, by allowing high level specification of accuracy requirements which are then met automatically without the need for the user to understand any algorithm-specific parameters. Additionally, we present a new suite of benchmarks, written in our language, to examine the efficacy of our techniques. Our experimental results show that by relaxing accuracy requirements, we can easily obtain performance improvements ranging from 1.1× to orders of magnitude of speedup.

It is often impossible to obtain a one-size-fits-all solution for high performance algorithms when considering different choices for data distributions, parallelism, transformations, and blocking. The best solution to these choices is often tightly coupled to different architectures, problem sizes, data, and available system resources. In some cases, completely different algorithms may provide the best performance. Current compiler and programming language techniques are able to change some of these parameters, but today there is no simple way for the programmer to express or the compiler to choose different algorithms to handle different parts of the data. Existing solutions normally can handle only coarse-grained, library level selections or hand coded cutoffs between base cases and recursive cases.

Algorithmic choice is essential in any problem domain to realizing optimal computational performance. Multigrid is a prime example: not only is it possible to make choices at the highest grid resolution, but a program can switch techniques as the problem is recursively attacked on coarser grid levels to take advantage of algorithms with different scaling behaviors. Additionally, users with different convergence criteria must experiment with parameters to yield a tuned algorithm that meets their accuracy requirements. Even after a tuned algorithm has been found, users often have to start all over when migrating from one machine to another.

The high performance MATLAB® user now has more choices than ever. Interactive Supercomputing's Star-P embraces this new world where, as an example, a MATLAB user who never wants to leave MATLAB might sit next to a C++ programmer at the office and both surf the Internet for the latest high speed FFT written in yet another language. The MATLAB of the past now becomes one browser into a bigger computational world. HPC users need this bigger world. Other "browsers" can be imagined. The open Star-P platform gives users options never before available to programmers who have traditionally enjoyed living exclusively inside a MATLAB environment.

Free probability provides tools and techniques for analyzing the eigen-spectra of large Hermitian random matrices. These stochastic eigen-analysis techniques have been invaluable in providing insight into the structure of sample covariance matrices. We briefly outline how these techniques can be used to analytically predict the spectrum of large sample covariance matrices. An eigen-inference application is briefly discussed.

Star-P is a unique technology offered by Interactive Supercomputing after nurturing at MIT. Star-P through its clever abstractions is solving the ease of use problem that has plagued supercomputing. Given that there have been around 30 parallel MATLABs including three other major offerings, Star-P must be way ahead to compete in the marketplace.Some of the innovative features of Star-P are the ability to program in MATLAB, hook in task parallel codes written using a processor free abstraction, hook in existing parallel codes, and obtain the performance that represents the HPC promise. All this is through a client/server interface. Other clients such as Python or R could be possible. The MATLAB, Python, or R becomes the "browser." If we make it look easy, it is because decades of parallel computing experience has taught us that it is not.This talk demonstrates the abstractions and innovations that make this possible.

The MVDR beamformer is the most extensively used array processing algorithm and involves inverting the sample covariance matrix. In the snapshot deficient scenario, when the number of sensors is greater than or approximately equal to the number of snapshots, the eigenvalues of the resulting sample covariance matrix are poorly conditioned. Diagonal loading is then applied to the sample covariance matrix. Expressions for the bias of the resulting MVDR beamformer outputs in the sidelobe region are presented that are exact for asymptotically large arrays. Numerical simulations confirm the accuracy of these asymptotic expressions when predicting the bias of the outputs of moderately large arrays.

MATLAB is one of the most widely used mathematical computing environments in technical computing. It is an interactive environment that provides high-performance computational routines and an easy-to-use, C-like scripting language. It started out as an interactive interface to EISPACK and LINPACK and has remained a serial program. In 1995, C. Moler of Mathworks argued that there was no market at the time for a parallel MATLAB. But times have changed and we are seeing increasing interest in developing a parallel MATLAB, from both academic and commercial sectors. In a recent survey, 27 parallel MATLAB projects have been identified. We expand upon that survey and discuss the approaches the projects have taken to parallelize MATLAB. Also, we describe innovative features in some of the parallel MATLAB projects. Then we will conclude with an idea of a "right" parallel MATLAB. Finally we will give an example of what we think is a "right" parallel MATLAB: MATLAB*P.

Article Free Access Share on Modeling and rendering of weathered stone Authors: Julie Dorsey Laboratory for Computer Science, Massachusetts Institute of Technology Laboratory for Computer Science, Massachusetts Institute of TechnologyView Profile , Alan Edelman Laboratory for Computer Science, Massachusetts Institute of Technology Laboratory for Computer Science, Massachusetts Institute of TechnologyView Profile , Henrik Wann Jensen Laboratory for Computer Science, Massachusetts Institute of Technology Laboratory for Computer Science, Massachusetts Institute of TechnologyView Profile , Justin Legakis Laboratory for Computer Science, Massachusetts Institute of Technology Laboratory for Computer Science, Massachusetts Institute of TechnologyView Profile , Hans Køhling Pedersen Laboratory for Computer Science, Massachusetts Institute of Technology Laboratory for Computer Science, Massachusetts Institute of TechnologyView Profile Authors Info & Claims SIGGRAPH '99: Proceedings of the 26th annual conference on Computer graphics and interactive techniquesJuly 1999 Pages 225–234https://doi.org/10.1145/311535.311560Online:01 July 1999Publication History 176citation1,637DownloadsMetricsTotal Citations176Total Downloads1,637Last 12 Months6Last 6 weeks1 Get Citation AlertsNew Citation Alert added!This alert has been successfully added and will be sent to:You will be notified whenever a record that you have chosen has been cited.To manage your alert preferences, click on the button below.Manage my AlertsNew Citation Alert!Please log in to your account Save to BinderSave to BinderCreate a New BinderNameCancelCreateExport CitationPublisher SiteeReaderPDF

We present a linear algebraic formulation for a class of index transformations such as Gray code encoding and decoding, matrix transpose, bit reversal, vector reversal, shuffles, and other index or dimension permutations. This formulation unifies, simplifies, and can be used to derive algorithms for hypercube multiprocessors. We show how all the widely known properties of Gray codes, and some not so well-known properties as well, can be derived using this framework. Using this framework, we relate hypercube communications algorithms to Gauss-Jordan elimination on a matrix of 0's and 1's.>

The authors have designed and implemented a hypercube algorithm for direct N-body solvers on the Connection Machine CM-2. The algorithm is optimal in the sense that, as long as there is sufficient data, it uses the full communication bandwidth of a hypercube of any dimension. When the number of bodies per node is large enough, the communication time for the implementation is negligible, i.e., less than 2%. In particular, this means that one obtains close to optimal speedup in the regime. To obtain this performance, 'rotated and translated Gray codes' which result in time-wise edge disjoint Hamiltonian paths on the hypercube are used. Timings are presented for a collection of interacting point vortices in two dimensions. The computation of the velocities of 14,000 vortices in 32-bit precision takes 2 seconds on a 16K CM-2.>