Documentation 
Numerical singular values of a matrix
This functionality does not run in MATLAB.
numeric::singularvalues(A, <Hard  HardwareFloats  Soft  SoftwareFloats>)
numeric::singularvalues(A) returns numerical singular values of the matrix A.
The singular values of an m×n matrix A are the p = min(m, n) real nonnegative square roots of the eigenvalues of A^{H} A (for p = n) or of A A^{H} (for p = m). The Hermitian transpose A^{H} is the complex conjugate of the transpose of A.
numeric::singularvalues returns a list of real singular values [d_{1}, …, d_{p}] sorted by numeric::sort, i.e., d_{1} ≥ … ≥ d_{p} ≥ 0.0.
All entries of A must be numerical. Numerical expressions such as etc. are accepted and converted to floats. Nonnumerical symbolic entries lead to an error.
Cat::Matrix objects, i.e., matrices A of a matrix domain such as Dom::Matrix(…) or Dom::SquareMatrix(…) are internally converted to arrays over expressions via expr(A).
Note: Singular values are approximated with an absolute precision of where r is largest singular value of A. Consequently, large singular values should be computed correctly to DIGITS decimal places. The numerical approximations of small singular values are less accurate. 
Singular values may also be computed via map ( numeric::eigenvalues( A A^{H}), sqrt ) or map ( numeric::eigenvalues( A^{H} A ), sqrt ), respectively. The use of numeric::singularvalues avoids the costs of the matrix multiplication. Further, the eigenvalue routine requires about twice as many DIGITS to compute small singular values with the same precision as numeric::singularvalues. Cf. Example 2.
The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.
The singular values of A and A^{H} coincide:
A := array(1..3, 1..2, [[1, 2*I], [2, 3],[3, sqrt(2)]]):
numeric::singularvalues(A)
The Hermitian transpose B = A^{H}:
B := array(1..2, 1..3, [[1, 2, 3], [2*I, 3, sqrt(2)]]):
numeric::singularvalues(B)
delete A, B:
We use numeric::eigenvalues to compute singular values:
A := matrix([[1/15, 2/15*I], [PI, 314159265358980/50000000000000*I], [2, 4*I]]):
The Hermitian transpose B = A^{H} can be computed by the methods conjugate and transpose of the matrix domain:
B := A::dom::conjugate(A::dom::transpose(A)):
Note that A^{H} A is positive semidefinite and cannot have negative eigenvalues. However, computing small eigenvalues is numerically illconditioned, and a small negative value occurs due to roundoff:
numeric::eigenvalues(B*A)
Consequently, an (incorrect) imaginary singular value is computed:
map(%, sqrt)
We have to increase DIGITS in order to compute this value more accurately:
DIGITS := 22: map(numeric::eigenvalues(B*A), sqrt)
With numeric::singularvalues, the standard precision suffices:
DIGITS := 10: numeric::singularvalues(A, SoftwareFloats)
numeric::singularvalues(A, HardwareFloats)
delete A, B:
We demonstrate the use of hardware floats. Hilbert matrices are notoriously illconditioned: the computation of the small singular values is subject to severe roundoff effects. In the following results, both with HardwareFloats as well as with SoftwareFloats, the small singular values are dominated by numerical roundoff. Consequently, the results with HardwareFloats differ from those with with SoftwareFloats:
numeric::singularvalues(linalg::hilbert(13))
[1.813830119, 0.396833076, 0.04902941942, 0.004348755075, 0.0002951777135, 0.0000156237036, 0.0000006466418563, 0.00000002076321421, 0.0000000005076551851, 9.141268657e12, 1.143562234e13, 8.877867351e16, 7.878607674e19]
A := linalg::hilbert(15): numeric::singularvalues(A, HardwareFloats); numeric::singularvalues(A, SoftwareFloats)
[1.845927746, 0.426627957, 0.05721209253, 0.005639834756, 0.0004364765944, 0.00002710853923, 0.000001361582242, 0.00000005528988481, 0.000000001802959758, 4.657785895e11, 9.321516341e13, 1.386205079e14, 1.463931556e16, 1.249693852e17, 6.620874158e18] [1.845927746, 0.426627957, 0.05721209253, 0.005639834756, 0.0004364765944, 0.00002710853923, 0.000001361582242, 0.00000005528988482, 0.000000001802959751, 4.657786546e11, 9.321608034e13, 1.39406179e14, 1.46659771e16, 9.68226589e19, 3.017915362e21]
delete A:
Hard, HardwareFloats, Soft, SoftwareFloats 
With Hard (or HardwareFloats), computations are done using fast hardware float arithmetic from within a MuPAD^{®} session. Hard and HardwareFloats are equivalent. With this option, the input data are converted to hardware floats and processed by compiled C code. The result is reconverted to MuPAD floats and returned to the MuPAD session. With Soft (or SoftwareFloats) computations are dome using software float arithmetic provided by the MuPAD kernel. Soft and SoftwareFloats are equivalent. SoftwareFloats is used by default if the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY. Compared to the SoftwareFloats used by the MuPAD kernel, the computation with HardwareFloats may be many times faster. Note, however, that the precision of hardware arithmetic is limited to about 15 digits. Further, the size of floatingpoint numbers may not be larger than approximately 10^{308} and not smaller than approximately 10^{ 308}. If no HardwareFloats or SoftwareFloats are requested explicitly, the following strategy is used: If the current value of DIGITS is smaller than 16 or if the matrix A is a hardware float array of domain type DOM_HFARRAY, then hardware arithmetic is tried. If this is successful, the result is returned. If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried. If the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY, or if one of the options Soft, SoftwareFloats or Symbolic is specified, MuPAD computes the result with its software arithmetic without trying to use hardware floats first. There may be several reasons for hardware arithmetic to fail:
If neither HardwareFloats nor SoftwareFloats is specified, the user is not informed whether hardware floats or software floats are used. If HardwareFloats are specified but fail due to one of the reasons above, a warning is issued that the (much slower) software floatingpoint arithmetic of the MuPAD kernel is used. Note that HardwareFloats can only be used if all input data can be converted to floatingpoint numbers. The trailing digits in floatingpoint results computed with HardwareFloats and SoftwareFloats may differ.

The code implements standard numerical algorithms from the Handbook of Automatic Computation by Wilkinson and Reinsch.