A strange-sounding rule from people who could write it themselves
You have spent whole rungs learning the machinery: LU factorization with pivoting, the QR algorithm for eigenvalues, the FFT in O(N log N), conjugate gradients for huge sparse systems. So it can feel like a betrayal to be told the first rule of professional scientific computing: don't roll your own. The point is not that you can't — you now genuinely understand these methods — but that re-implementing them in production code is almost always the wrong use of your effort, and quietly dangerous. The people who built the standard libraries say this loudest, because they know exactly how many ways there are to get the details subtly wrong.
Think of it like surgery. Understanding anatomy is essential — but understanding anatomy is not the same as deciding to operate on yourself. The textbook version of Gaussian elimination you wrote to learn it is a fine teaching tool and a fine prototype. The version inside LAPACK has had four decades of numerical analysts, compiler experts and hardware engineers stress-testing every corner. That gap is the whole subject of this guide.
The stack, layer by layer
What you actually stand on is a numerical software stack: a tower of layers, each one trusting the layer below and hiding it from the layer above. At the very bottom sits BLAS, the Basic Linear Algebra Subprograms — the tiny vocabulary of vector and matrix operations (dot products, matrix-vector and matrix-matrix products) that everything else is phrased in. One layer up is LAPACK, which expresses the big algorithms — solving A x = b, least squares, eigenvalues, the SVD — entirely as sequences of BLAS calls. Above that sit the friendly environments you type into: NumPy and SciPy in Python, the built-in array language of Julia, MATLAB, R. When you write a one-line solve in any of them, you are very often calling straight down into the same LAPACK and BLAS underneath.
your code: numpy.linalg.solve(A, b) <- one line you write ------------------------------------------------------------- SciPy / NumPy / Julia / MATLAB <- arrays, slicing, plots LAPACK dgesv, dgeqrf, dsyev, dgesdd <- the algorithms BLAS ddot, dgemv, dgemm <- the vocabulary hardware SIMD registers, cache, cores <- where flops happen
There is a beautiful engineering reason for splitting it this way. BLAS comes in three levels: Level 1 is vector-vector work (O(n) data, O(n) flops), Level 2 is matrix-vector (O(n^2) of each), and Level 3 is matrix-matrix (O(n^2) data but O(n^3) flops). Only Level 3 does enough arithmetic per number fetched from memory to keep a modern processor busy — recall from the HPC rung that a memory-bound kernel is limited by cache misses, not flops. So LAPACK is deliberately rewritten to spend as much time as possible inside Level-3 matrix-matrix calls. Vendors then ship a hand-tuned BLAS for their exact chip, and the same LAPACK code suddenly runs ten times faster. That clean seam — algorithm above, hardware-specific speed below — is why the split has survived for forty years.
What forty years of bug-hunting actually buys you
The deepest reason not to roll your own is not speed — it is correctness in the floating-point sense you learned earlier. Recall the master equation of this whole field: accuracy = conditioning × stability. Conditioning belongs to the problem; you can't change it. Stability belongs to the algorithm, and it is brutally easy to lose. A library routine like LAPACK's solver is backward stable: the answer it returns is the exact solution of a problem only a rounding-error's distance from the one you asked. That is the strongest honesty guarantee a finite-precision algorithm can offer, and proving it took genuine research. Your hand-rolled loop that skips pivoting, or sums in a careless order, gives you no such guarantee — and it fails silently, returning a plausible-looking wrong answer rather than crashing.
Be honest about the limit, though, because this is where beginners over-trust the library. Backward stability does NOT mean the answer is accurate. If your problem is ill-conditioned — say the matrix has a condition number near 10^8 — then even a perfect, backward-stable solver loses about 8 of your ~16 double-precision digits, purely because the problem amplifies the unavoidable input rounding. The library did its job flawlessly; the problem was hard. A good library hands you the condition number so you can see this coming; a hand-rolled solver usually leaves you blind to it. The library buys you a trustworthy algorithm, not a trustworthy problem.
Where the libraries end and your judgement begins
"Don't roll your own" is a rule about implementation, not about understanding — and it has real edges. The library will hand you the right answer to the question you literally asked, but it cannot know whether that was the right question, nor whether you fed it nonsense. Everything you learned about the pitfalls still applies above the library line. If you ask for a high-degree polynomial fit on equally spaced points, a perfect library still suffers Runge's phenomenon and the fit diverges at the edges; the cure is to choose Chebyshev nodes or splines, a modelling decision the library cannot make for you. If you shrink a finite-difference step h forever expecting more accuracy, the round-off floor eventually makes the answer worse, library or not.
There are also genuine times to write your own. If your problem has special structure no library exposes — an unusual sparsity pattern, a custom preconditioner, a matrix-free operator you only know how to apply, not store — you may legitimately assemble your own outer loop. But notice the craft: you build the novel part and still call BLAS, LAPACK or a sparse package for every standard step inside it. You roll the wrapper, never the kernel. And when you do, you treat your code with the same suspicion you'd give a stranger's — which is exactly what verification, the subject of the next guide, is for.
How to choose a tool without losing the thread
The three environments you'll meet most carry the same engine but feel different in the hand. NumPy and SciPy give you Python's huge ecosystem with thin, fast wrappers over the C and Fortran libraries — most of your time is spent in compiled code, so the slow Python interpreter rarely bites, provided you keep operations vectorized over arrays rather than looping element by element. Julia was designed from scratch for this work: it just-in-time compiles your own high-level code down to near-C speed, so the painful old split between a slow prototyping language and a fast "real" language softens. MATLAB is the classic numerical workbench, polished and widely taught, sitting on the very same LAPACK. There is no universally best choice — only the right one for your team, your libraries, and your problem.
Pull it together and the rule stops sounding like a betrayal. You learned the algorithms so you would know which library routine to reach for, what its result means, and where it can fail — not so you would retype them. The rest of this rung builds on exactly this footing: knowing the stack underneath you is trustworthy, you can spend your scarce attention on the questions a library will never answer for you — Did I solve the equations right? Did I solve the right equations? Could someone reproduce this? And what is my honest error bar?