I'm only interested in finding the largest singular value of an $m\times n$ real matrix $A$. I don't need the corresponding (left and right) singular vectors.
For matrices of general form and not especially small rank, we should not expect a direct method to locate eigenvalues or singular values "analytically". Some discussion of this is here, but the short explanation is that it would be equivalent to a formula for polynomial roots of arbitrary degree. We instead describe an "easy" iterative approach.
The now-deleted Answer by KCPP said simply, "use power method which is the simplest". We show that this is the germ of an elementary algorithm to approximate the largest singular value of a (not necessarily square) matrix $A$. [This numerical topic has been previously discussed at the Computational Science sister site, along with the more sophisticated Lanczos algorithm.]
The nonzero singular values of $A$ are square roots of the nonzero eigenvalues of $A^*A$. So the power method can be applied to Hermitian matrix $A^*A$ to iteratively approximate its maximum (real) eigenvalue, thus producing as the square root of that the maximum singular value of $A$.
For simplicity of notation, I'll assume $A$ itself is a real matrix of size $m\times n$ with $m\ge n$ as in Vini's Answer. However the basic idea would apply (forming Rayleigh quotients) using complex arithmetic and $A^*$ if instead the matrix was not "real and tall".
Generate a "random" initial vector $x_0$, say by picking components from a normal distribution. The goal here is simply to have a good chance that the initial vector is not "independent" of the maximum value eigenspace of $A^T A$. One then iterates:
$$ x_{i+1} := \frac{A^T(A x_i)}{||A^T(A x_i)||} $$
Convergence of this iteration is geometric in the ratio of next-largest eigenvalue of $A^T A$ to the largest such eigenvalue. If the approximation of the right singular vector were of interest, the sequence of $x_i$ would serve this purpose (as $Ax_i$ would approximate the left singular vector). Setting that aside, the Rayleigh quotients $||A^T(A x_i)||/||x_i||$ are lower bounds on the largest eigenvalue $\lambda_{\max}$ of $A^T A$, and their convergence is ordinarily easy to detect (although misconvergence with unlucky initial vectors is possible).
Finally we take a square root of the converged approximation to $\lambda_{\max}$ as our estimate of the largest singular value of $A$.
The standard way to solve the SVD problem is to reduce the $m \times n$ matrix $A$ (we are assuming that $m \ge n$, otherwise work with the transpose of $A$) to a bidiagonal matrix using orthogonal transformations such that $$ U^t A V = B, \;\;\; U^tU = I, \;\;\; V^tV = I $$ where $B$ is an $n \times n$ bidiagonal matrix. If all the singular values are required then the QR-algorithm or the divide and conquer algorithm is used on the matrix $B$; see the standard textbook "Matrix Computations" by Gene Golub and Charles Van Loan. However, if only a few (or one) singular values is required, bisection is used.
