*> The truncation criteria (i.e. when to stop the factorization)
*> can be any of the following:
*>
*>   1) The input parameter KMAX, the maximum number of columns
*>      KMAX to factorize, i.e. the factorization rank is limited
*>      to KMAX. If KMAX >= min(M,N), the criterion is not used.
*>
*>   2) The input parameter ABSTOL, the absolute tolerance for
*>      the maximum column 2-norm of the residual matrix R22(K). This
*>      means that the factorization stops if this norm is less or
*>      equal to ABSTOL. If ABSTOL < 0.0, the criterion is not used.
*>
*>   3) The input parameter RELTOL, the tolerance for the maximum
*>      column 2-norm matrix of the residual matrix R22(K) divided
*>      by the maximum column 2-norm of the original matrix A, which
*>      is equal to abs(R(1,1)). This means that the factorization stops
*>      when the ratio of the maximum column 2-norm of R22(K) to
*>      the maximum column 2-norm of A is less than or equal to RELTOL.
*>      If RELTOL < 0.0, the criterion is not used.
*>
*>   4) In case both stopping criteria ABSTOL or RELTOL are not used,
*>      and when the residual matrix R22(K) is a zero matrix in some
*>      factorization step K. ( This stopping criterion is implicit. )
*>
*>  The algorithm stops when any of these conditions is first
*>  satisfied, otherwise the whole matrix A is factorized.


*> =====================================================================
*> {S,D,C,Z}GEQP3RK performs two tasks simultaneously:
*>
*> Task 1: The routine computes a truncated (rank K) or full rank
*> Householder QR factorization with column pivoting of a real
*> M-by-N matrix A using Level 3 BLAS. K is the number of columns
*> that were factorized, i.e. factorization rank of the
*> factor R, K <= min(M,N).
*>
*>  A * P(K) = Q(K) * R(K)  =
*>
*>        = Q(K) * ( R11(K) R12(K) ) = Q(K) * (   R(K)_approx    )
*>                 ( 0      R22(K) )          ( 0  R(K)_residual ),
*>
*> where:
*>
*>  P(K)            is an N-by-N permutation matrix;
*>  Q(K)            is an M-by-M unitary matrix;
*>  R(K)_approx   = ( R11(K), R12(K) ) is a rank K approximation of the
*>                    full rank factor R with K-by-K upper-triangular
*>                    R11(K) and K-by-N rectangular R12(K). The diagonal
*>                    entries of R11(K) appear in non-increasing order
*>                    of absolute value, and absolute values of all of
*>                    them exceed the maximum column 2-norm of R22(K)
*>                    up to roundoff error.
*>  R(K)_residual = R22(K) is the residual of a rank K approximation
*>                    of the full rank factor R. It is a
*>                    an (M-K)-by-(N-K) rectangular matrix;
*>  0               is a an (M-K)-by-K zero matrix.
*>
*> Task 2: At the same time, the routine overwrites a real M-by-NRHS
*> matrix B with  Q(K)**T * B  using Level 3 BLAS.
*>
*> =====================================================================