(double precision real only)

- add file to CMakeLists.txt
- fix syntax errors etc.

- workspace queries are signalled by LWORK=-1, not INFO=-1
- fix undeclared, uninitialized variable errors

- fix workspace queries
- fix DORCSD2BY1 parameter containing the number of rows of X11
- avoid scaling the matrix B if B is zero
- handle the case where both matrices A, B are zero
- check INFO after calling DORGQR

Fix the computation of an n-by-n orthogonal matrix Q by DORGRQ by
ensuring that the R elementary reflectors from DGERQF can be found in
the last R rows of Q.

Ensure the number of columns indicated to DORCSD2BY1 is less than or
equal to the number of rows during workspace queries.

The variable FACTOR was copy-pasted from xGEEQUB.

- compute Q only when needed
- fix argument order to DORCSD2BY1

DLANGE does not use WORK when computing the Frobenius norm.
Nevertheless, this change may avoid future bugs.

Return the optimal workspace size even if both matrices A, B are zero.

The multiplication uses DGEMM so the lower triangular part of
R1( 1:RANK, : ) must be set to zero.

Ease debugging

This bug was introduced in commit 446ac9a2.

This bug was discovered by the random test from commit dcd04cd.
Minimal triggering example with m=1, n=2, p=1 with non-random entries:
  A = [1, 1],
  B = [1, 0].

Fix harmless out-of-bounds accesses to avoid spurious address sanitizer
(ASAN) failures.

* fix typos
* explain optimal workspace size computation

* fix an argument type
* fix name of called function

Remove integer variable `R` indicating the rank because
* integer `L` was already introduced for this purpose
* the purpose of `R` was never documented, and
* it can be confused with the triangular matrices occuring throghout the
  computation.

For CUNCSD2BY1, memory consumption is not at maximum with maximum matrix
dimensions.

In an attempt to avoid a large backward error, xGGQRCS scales on of the
input matrices so that both input matrices have comparable norm. This
changes (not perturb, changes!) the generalized singular values. The
previously committed test focusing on the singular values computed by
xGGQRCS highlighted the need to fix the singular values.

With this commit, it became obvious that
* the singular values can be computed to very high relative accuracy
  (small forward error),
* the matrices computed by xGGQRCS change significantly with scaling
  (large backward error).

Consequently, matrix scaling cannot be used and row sorting must be
applied to the matrix `(A, B)` before the initial QR decomposition.

xLASRT sorts an array of real values. xLASRTI sorts an array of indices
referencing real values in an array.

SLASRTI is based on SLASRT.

The indexing error could lead to an infinite loop.

The GSVD decomposes a pair of matrices A, B into
* `A = U1 D1 X`,
* `B = U2 D2 X`,

where `X` has full rank. Alternatively, one can compute the equivalent
decomposition
* `A = U1 D1 R Q^*`,
* `B = U2 D2 R Q^*`,

where `R` is upper triangular and `Q` orthogonal. The second form has
several advantages from the point of view of _numerical_ linear algebra.
None of these apply to the GSVD solvers based on QR and CS
decomposition.  Consequently, xGGQRCS returns from now on only the
product `X = RQ^*`.

This change makes xGGQRCS
* faster,
* more flexible (`X` is computed or not),
* easier to implement (no need to assemble `R Q^*` for tests), and
* numerically stable!

xGGQRCS becomes numerically stable because factoring `X` requires an
explicit LQ decomposition but the result may not be backward stable
unless `A` and `B` are similar in norm. Not computing the LQ
decomposition obviously solves this problem.

This commit is a proof-of-concept; {d,c,z}GGQRCS should be updated
accordingly.

To achieve a backward stable solver, input matrices A, B were previously
scaled to be similar norm but this changes the generalized singular
values and I seem unable to compensate for the scaling effects by
* recomputing the singular values and
* rescaling the right-hand side GSVD matrix X.

Row sorting cannot replace matrix scaling, see the recently added tests.

The singular value will not change but the corresponding row in the
right-hand side GSVD matrix X must be scaled.

* highlight conditional backward stability
* recommend use of SGGSVD3 when necessary
* emphasize accuracy of computed singular values

Row sorting has no discernible effect on forward or backward errors.

The singular values must be corrected for matrix scaling. Let w be a
power of two such that norm(wA) equals approximately norm(B). Then the
radians representation of the singular values must corrected with

x' = arctan(w * tan(x)),

where x is the computed radians value with matrix scaling and x' is the
computed radians value without matrix scaling. If w is less than one
and if x is large (close to pi/2), then roughly speaking

* tan(x) has a large derivative (i.e. it is badly conditioned),
* w * tan(x) will be smaller than tan(x),
* arctan(w * tan(x)) has derivative near one, because its argument
  is closer to zero.

This commit modifies the GSVD computation such that w >= 1 by choosing
to compute the GSVD of (A, B) or (B, A) depending on the relative norms.
Then

* tan(x) has a small derivative if x is sufficiently small,
* arctan(w * tan(x)) has derivative near zero, and
* if x is large (near pi/2), then x' will be near pi/2 as well.

This commit fixes the large relative forward errors detected by
xGGQRCS_test_singular_values<float>' when the singular values of the
matrix pair (A,B) where drawn from [0, pi/2000).

Previously SGGQRCS returned angles instead of sine, cosine values but
these computations (while numerically stable) may completely destroy the
backward stability of the solver. Specifically, in one test* a computed
angle x near π/2 was the best single-precision approximation to the true
value but cosine(x) was nowhere near the true cosine value, thereby
making it impossible to achieve a backward stable reassembly of the
input matrices from the GSVD.

For this reason, SGGQRCS returns sine and cosine values to the caller
instead of angles.

TODO: Documentation

*Test name  xGGQRCS_test_singular_accuracy_vs_radians_accuracy

Port the most recent SGGQRCS code to double precision.

Port the most recent SGGQRCS code to single-precision complex.

Port the most recent CGGQRCS code to double-precision complex.

For tall matrices, xGGSVD3 is faster than xGGQRCS. This commit tries to
improve performance by speeding up the matrix-matrix multiplication
`V^* R( 1:L, N )`, where `R` is upper triangular with, with `xTRMM()`.

Revert because
* the commit was incomplete and forgot to remove the superfluous SLASCL
  call zeroing the lower triangular matrix part; with this called
  removed...
* there is hardly any difference for matrices with m, n, p <= 256,
* if there is a difference, it may be in favor of xGEMM.

Fix out-of-bounds access when computing only the singular values.

The last element of the matrix X overlapped with the first element of
the matrix V^* because WORK(1) cannot be used by X.

This commit is a port of recent changes mades to CGGQRCS.

* fix typos
* remove spurious direction value for an internal variable
* fix formatting of internal variables section

Ensure the leading dimension of the assembled matrix `G = [A; B]` is
always at least one.