For instance, there a fast iterative solvers. They start with an initial guess for the solution (which can be an arbitrary vector) and successively generate more accurate approximations.

These iterative solvers don't need access to the inverse of A. Furthermore, many of those solvers don't even need access to A: They only need access to a subroutine which, given a vector u, returns Au.

Many matrices arising in applications (for instance one often reduces the problem of solving partial differential equations to the problem of solving a (very huge) linear system of equations) are sparse (have only few non-zero entries). For these it can make sense to write specialized code which specifically calculates Au.

If you want to read more about this, then a good starting point is https://en.wikipedia.org/wiki/Kaczmarz_method.

    A x = b
    A => QR  # QR decomposition
    Q R x = b
    R x = Q' b  # Q is inverted by transposing
    x = backsubstitute(R, Q' b)

Now you can let Ux = y, solve the system Ly = b for y, then solve Ux = y for x.

Solving each triangular system is fast: O(n^2).

If you need to do this with a different b, you can do it again, and you already have A factorized, which was the expensive part, i.e. O(n^3).

(Depending on the type of problem, you might want to factorize the original matrix in a different way, e.g. QR or Cholesky or whatever.)

Lets say you are using QR factorization. Then R\b is easy to compute by back-substitution, because R is triangular. Q\b is also easy, because Q is orthogonal, and Q's inverse is its transpose, so you just calculate transpose(Q)*b.

Same thing for LDL^T factorization, LU factorization, etc. Its easy because back-substitution with a triangular matrix is easy. In the case of LDL^T factorization, D is diagonal, and again D\b is easy (its a bunch of rank-1 division operations).

No.

>How else would you solve Ax = b for x?

Call dgetrs instead of dgetri.

You use the factorisation directly (not the inverse) to get x from b. This is done by calling dgetrf to get the factorisation, then calling dgetrs to get the solution. You don't need to call dgetri.

I'd do an example, but there are other replies to your comment that go through the mathematics.

These solvers construct a solution by adding up terms of the form c_n * A^n * y and truncating the series when the residuals are small enough. This is especially efficient in PDEs, where problem matrices are sparse so matrix-vector products are efficient (i.e. O(n)), and in practice only a few terms of the series are necessary.

Incidentally, while solving `Ax = b` in general is definitely equivalent to inverting `A` (simply by taking `b` to be basis vectors, so that the solutions `x` are the columns of the inverse matrix in that basis), it is certainly possible to solve it in particular cases without inverting `A`; for example, the solution to `Ax = b`, where `b` is the first column of `A` in a particular basis, is the first vector `x` of that basis—no inversion required (nor even computation of `b`)!

Numerical linear algebra has a long and vast history. You really do need specialized information that your undergrad linear algebra classes would not have imparted.

For instance, the condition number is important: https://en.wikipedia.org/wiki/Condition_number

min_x* ( ||A x* - b|| (+\lamba regularization(x*) )