LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ sgeqp3rk()

subroutine sgeqp3rk ( integer m,
integer n,
integer nrhs,
integer kmax,
real abstol,
real reltol,
real, dimension( lda, * ) a,
integer lda,
integer k,
real maxc2nrmk,
real relmaxc2nrmk,
integer, dimension( * ) jpiv,
real, dimension( * ) tau,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer info )

SGEQP3RK computes a truncated Householder QR factorization with column pivoting of a real m-by-n matrix A by using Level 3 BLAS and overwrites a real m-by-nrhs matrix B with Q**T * B.

Download SGEQP3RK + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SGEQP3RK 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 orthogonal 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.
!>
!> =====================================================================
!>
!> The matrices A and B are stored on input in the array A as
!> the left and right blocks A(1:M,1:N) and A(1:M, N+1:N+NRHS)
!> respectively.
!>
!>                                  N     NRHS
!>             array_A   =   M  [ mat_A, mat_B ]
!>
!> 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.
!>
!>  To factorize the whole matrix A, use the values
!>  KMAX >= min(M,N), ABSTOL < 0.0 and RELTOL < 0.0.
!>
!>  The routine returns:
!>     a) Q(K), R(K)_approx = ( R11(K), R12(K) ),
!>        R(K)_residual = R22(K), P(K), i.e. the resulting matrices
!>        of the factorization; P(K) is represented by JPIV,
!>        ( if K = min(M,N), R(K)_approx is the full factor R,
!>        and there is no residual matrix R(K)_residual);
!>     b) K, the number of columns that were factorized,
!>        i.e. factorization rank;
!>     c) MAXC2NRMK, the maximum column 2-norm of the residual
!>        matrix R(K)_residual = R22(K),
!>        ( if K = min(M,N), MAXC2NRMK = 0.0 );
!>     d) RELMAXC2NRMK equals MAXC2NRMK divided by MAXC2NRM, the maximum
!>        column 2-norm of the original matrix A, which is equal
!>        to abs(R(1,1)), ( if K = min(M,N), RELMAXC2NRMK = 0.0 );
!>     e) Q(K)**T * B, the matrix B with the orthogonal
!>        transformation Q(K)**T applied on the left.
!>
!> The N-by-N permutation matrix P(K) is stored in a compact form in
!> the integer array JPIV. For 1 <= j <= N, column j
!> of the matrix A was interchanged with column JPIV(j).
!>
!> The M-by-M orthogonal matrix Q is represented as a product
!> of elementary Householder reflectors
!>
!>     Q(K) = H(1) *  H(2) * . . . * H(K),
!>
!> where K is the number of columns that were factorized.
!>
!> Each H(j) has the form
!>
!>     H(j) = I - tau * v * v**T,
!>
!> where 1 <= j <= K and
!>   I    is an M-by-M identity matrix,
!>   tau  is a real scalar,
!>   v    is a real vector with v(1:j-1) = 0 and v(j) = 1.
!>
!> v(j+1:M) is stored on exit in A(j+1:M,j) and tau in TAU(j).
!>
!> See the Further Details section for more information.
!>
Parameters
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A. M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A. N >= 0.
!>
[in]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e. the number of
!>          columns of the matrix B. NRHS >= 0.
!>
[in]KMAX
!>          KMAX is INTEGER
!>
!>          The first factorization stopping criterion. KMAX >= 0.
!>
!>          The maximum number of columns of the matrix A to factorize,
!>          i.e. the maximum factorization rank.
!>
!>          a) If KMAX >= min(M,N), then this stopping criterion
!>                is not used, the routine factorizes columns
!>                depending on ABSTOL and RELTOL.
!>
!>          b) If KMAX = 0, then this stopping criterion is
!>                satisfied on input and the routine exits immediately.
!>                This means that the factorization is not performed,
!>                the matrices A and B are not modified, and
!>                the matrix A is itself the residual.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>
!>          The second factorization stopping criterion, cannot be NaN.
!>
!>          The absolute tolerance (stopping threshold) for
!>          maximum column 2-norm of the residual matrix R22(K).
!>          The algorithm converges (stops the factorization) when
!>          the maximum column 2-norm of the residual matrix R22(K)
!>          is less than or equal to ABSTOL. Let SAFMIN = SLAMCH('S').
!>
!>          a) If ABSTOL is NaN, then no computation is performed
!>                and an error message ( INFO = -5 ) is issued
!>                by XERBLA.
!>
!>          b) If ABSTOL < 0.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on KMAX and RELTOL.
!>                This includes the case ABSTOL = -Inf.
!>
!>          c) If 0.0 <= ABSTOL < 2*SAFMIN, then ABSTOL = 2*SAFMIN
!>                is used. This includes the case ABSTOL = -0.0.
!>
!>          d) If 2*SAFMIN <= ABSTOL then the input value
!>                of ABSTOL is used.
!>
!>          Let MAXC2NRM be the maximum column 2-norm of the
!>          whole original matrix A.
!>          If ABSTOL chosen above is >= MAXC2NRM, then this
!>          stopping criterion is satisfied on input and routine exits
!>          immediately after MAXC2NRM is computed. The routine
!>          returns MAXC2NRM in MAXC2NORMK,
!>          and 1.0 in RELMAXC2NORMK.
!>          This includes the case ABSTOL = +Inf. This means that the
!>          factorization is not performed, the matrices A and B are not
!>          modified, and the matrix A is itself the residual.
!>
[in]RELTOL
!>          RELTOL is REAL
!>
!>          The third factorization stopping criterion, cannot be NaN.
!>
!>          The tolerance (stopping threshold) for the ratio
!>          abs(R(K+1,K+1))/abs(R(1,1)) of the maximum column 2-norm of
!>          the residual matrix R22(K) to the maximum column 2-norm of
!>          the original matrix A. The algorithm converges (stops the
!>          factorization), when abs(R(K+1,K+1))/abs(R(1,1)) A is less
!>          than or equal to RELTOL. Let EPS = SLAMCH('E').
!>
!>          a) If RELTOL is NaN, then no computation is performed
!>                and an error message ( INFO = -6 ) is issued
!>                by XERBLA.
!>
!>          b) If RELTOL < 0.0, then this stopping criterion is not
!>                used, the routine factorizes columns depending
!>                on KMAX and ABSTOL.
!>                This includes the case RELTOL = -Inf.
!>
!>          c) If 0.0 <= RELTOL < EPS, then RELTOL = EPS is used.
!>                This includes the case RELTOL = -0.0.
!>
!>          d) If EPS <= RELTOL then the input value of RELTOL
!>                is used.
!>
!>          Let MAXC2NRM be the maximum column 2-norm of the
!>          whole original matrix A.
!>          If RELTOL chosen above is >= 1.0, then this stopping
!>          criterion is satisfied on input and routine exits
!>          immediately after MAXC2NRM is computed.
!>          The routine returns MAXC2NRM in MAXC2NORMK,
!>          and 1.0 in RELMAXC2NORMK.
!>          This includes the case RELTOL = +Inf. This means that the
!>          factorization is not performed, the matrices A and B are not
!>          modified, and the matrix A is itself the residual.
!>
!>          NOTE: We recommend that RELTOL satisfy
!>                min( max(M,N)*EPS, sqrt(EPS) ) <= RELTOL
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N+NRHS)
!>
!>          On entry:
!>
!>          a) The subarray A(1:M,1:N) contains the M-by-N matrix A.
!>          b) The subarray A(1:M,N+1:N+NRHS) contains the M-by-NRHS
!>             matrix B.
!>
!>                                  N     NRHS
!>              array_A   =   M  [ mat_A, mat_B ]
!>
!>          On exit:
!>
!>          a) The subarray A(1:M,1:N) contains parts of the factors
!>             of the matrix A:
!>
!>            1) If K = 0, A(1:M,1:N) contains the original matrix A.
!>            2) If K > 0, A(1:M,1:N) contains parts of the
!>            factors:
!>
!>              1. The elements below the diagonal of the subarray
!>                 A(1:M,1:K) together with TAU(1:K) represent the
!>                 orthogonal matrix Q(K) as a product of K Householder
!>                 elementary reflectors.
!>
!>              2. The elements on and above the diagonal of
!>                 the subarray A(1:K,1:N) contain K-by-N
!>                 upper-trapezoidal matrix
!>                 R(K)_approx = ( R11(K), R12(K) ).
!>                 NOTE: If K=min(M,N), i.e. full rank factorization,
!>                       then R_approx(K) is the full factor R which
!>                       is upper-trapezoidal. If, in addition, M>=N,
!>                       then R is upper-triangular.
!>
!>              3. The subarray A(K+1:M,K+1:N) contains (M-K)-by-(N-K)
!>                 rectangular matrix R(K)_residual = R22(K).
!>
!>          b) If NRHS > 0, the subarray A(1:M,N+1:N+NRHS) contains
!>             the M-by-NRHS product Q(K)**T * B.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A. LDA >= max(1,M).
!>          This is the leading dimension for both matrices, A and B.
!>
[out]K
!>          K is INTEGER
!>          Factorization rank of the matrix A, i.e. the rank of
!>          the factor R, which is the same as the number of non-zero
!>          rows of the factor R. 0 <= K <= min(M,KMAX,N).
!>
!>          K also represents the number of non-zero Householder
!>          vectors.
!>
!>          NOTE: If K = 0, a) the arrays A and B are not modified;
!>                          b) the array TAU(1:min(M,N)) is set to ZERO,
!>                             if the matrix A does not contain NaN,
!>                             otherwise the elements TAU(1:min(M,N))
!>                             are undefined;
!>                          c) the elements of the array JPIV are set
!>                             as follows: for j = 1:N, JPIV(j) = j.
!>
[out]MAXC2NRMK
!>          MAXC2NRMK is REAL
!>          The maximum column 2-norm of the residual matrix R22(K),
!>          when the factorization stopped at rank K. MAXC2NRMK >= 0.
!>
!>          a) If K = 0, i.e. the factorization was not performed,
!>             the matrix A was not modified and is itself a residual
!>             matrix, then MAXC2NRMK equals the maximum column 2-norm
!>             of the original matrix A.
!>
!>          b) If 0 < K < min(M,N), then MAXC2NRMK is returned.
!>
!>          c) If K = min(M,N), i.e. the whole matrix A was
!>             factorized and there is no residual matrix,
!>             then MAXC2NRMK = 0.0.
!>
!>          NOTE: MAXC2NRMK in the factorization step K would equal
!>                R(K+1,K+1) in the next factorization step K+1.
!>
[out]RELMAXC2NRMK
!>          RELMAXC2NRMK is REAL
!>          The ratio MAXC2NRMK / MAXC2NRM of the maximum column
!>          2-norm of the residual matrix R22(K) (when the factorization
!>          stopped at rank K) to the maximum column 2-norm of the
!>          whole original matrix A. RELMAXC2NRMK >= 0.
!>
!>          a) If K = 0, i.e. the factorization was not performed,
!>             the matrix A was not modified and is itself a residual
!>             matrix, then RELMAXC2NRMK = 1.0.
!>
!>          b) If 0 < K < min(M,N), then
!>                RELMAXC2NRMK = MAXC2NRMK / MAXC2NRM is returned.
!>
!>          c) If K = min(M,N), i.e. the whole matrix A was
!>             factorized and there is no residual matrix,
!>             then RELMAXC2NRMK = 0.0.
!>
!>         NOTE: RELMAXC2NRMK in the factorization step K would equal
!>               abs(R(K+1,K+1))/abs(R(1,1)) in the next factorization
!>               step K+1.
!>
[out]JPIV
!>          JPIV is INTEGER array, dimension (N)
!>          Column pivot indices. For 1 <= j <= N, column j
!>          of the matrix A was interchanged with column JPIV(j).
!>
!>          The elements of the array JPIV(1:N) are always set
!>          by the routine, for example, even  when no columns
!>          were factorized, i.e. when K = 0, the elements are
!>          set as JPIV(j) = j for j = 1:N.
!>
[out]TAU
!>          TAU is REAL array, dimension (min(M,N))
!>          The scalar factors of the elementary reflectors.
!>
!>          If 0 < K <= min(M,N), only the elements TAU(1:K) of
!>          the array TAU are modified by the factorization.
!>          After the factorization computed, if no NaN was found
!>          during the factorization, the remaining elements
!>          TAU(K+1:min(M,N)) are set to zero, otherwise the
!>          elements TAU(K+1:min(M,N)) are not set and therefore
!>          undefined.
!>          ( If K = 0, all elements of TAU are set to zero, if
!>          the matrix A does not contain NaN. )
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          LWORK >= 1, if MIN(M,N) = 0, and
!>          LWORK >= (3*N+NRHS-1), otherwise.
!>          For optimal performance LWORK >= (2*N + NB*( N+NRHS+1 )),
!>          where NB is the optimal block size for SGEQP3RK returned
!>          by ILAENV. Minimal block size MINNB=2.
!>
!>          NOTE: The decision, whether to use unblocked BLAS 2
!>          or blocked BLAS 3 code is based not only on the dimension
!>          LWORK of the availbale workspace WORK, but also also on the
!>          matrix A dimension N via crossover point NX returned
!>          by ILAENV. (For N less than NX, unblocked code should be
!>          used.)
!>
!>          If LWORK = -1, then a workspace query is assumed;
!>          the routine only calculates the optimal size of the WORK
!>          array, returns this value as the first entry of the WORK
!>          array, and no error message related to LWORK is issued
!>          by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (N-1).
!>          Is a work array. ( IWORK is used to store indices
!>          of  columns for norm downdating in the residual
!>          matrix in the blocked step auxiliary subroutine SLAQP3RK ).
!>
[out]INFO
!>          INFO is INTEGER
!>          1) INFO = 0: successful exit.
!>          2) INFO < 0: if INFO = -i, the i-th argument had an
!>                       illegal value.
!>          3) If INFO = j_1, where 1 <= j_1 <= N, then NaN was
!>             detected and the routine stops the computation.
!>             The j_1-th column of the matrix A or the j_1-th
!>             element of array TAU contains the first occurrence
!>             of NaN in the factorization step K+1 ( when K columns
!>             have been factorized ).
!>
!>             On exit:
!>             K                  is set to the number of
!>                                   factorized columns without
!>                                   exception.
!>             MAXC2NRMK          is set to NaN.
!>             RELMAXC2NRMK       is set to NaN.
!>             TAU(K+1:min(M,N))  is not set and contains undefined
!>                                   elements. If j_1=K+1, TAU(K+1)
!>                                   may contain NaN.
!>          4) If INFO = j_2, where N+1 <= j_2 <= 2*N, then no NaN
!>             was detected, but +Inf (or -Inf) was detected and
!>             the routine continues the computation until completion.
!>             The (j_2-N)-th column of the matrix A contains the first
!>             occurrence of +Inf (or -Inf) in the factorization
!>             step K+1 ( when K columns have been factorized ).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!> SGEQP3RK is based on the same BLAS3 Householder QR factorization
!> algorithm with column pivoting as in SGEQP3 routine which uses
!> SLARFG routine to generate Householder reflectors
!> for QR factorization.
!>
!> We can also write:
!>
!>   A = A_approx(K) + A_residual(K)
!>
!> The low rank approximation matrix A(K)_approx from
!> the truncated QR factorization of rank K of the matrix A is:
!>
!>   A(K)_approx = Q(K) * ( R(K)_approx ) * P(K)**T
!>                        (     0     0 )
!>
!>               = Q(K) * ( R11(K) R12(K) ) * P(K)**T
!>                        (      0      0 )
!>
!> The residual A_residual(K) of the matrix A is:
!>
!>   A_residual(K) = Q(K) * ( 0              0 ) * P(K)**T =
!>                          ( 0  R(K)_residual )
!>
!>                 = Q(K) * ( 0        0 ) * P(K)**T
!>                          ( 0   R22(K) )
!>
!> The truncated (rank K) factorization guarantees that
!> the maximum column 2-norm of A_residual(K) is less than
!> or equal to MAXC2NRMK up to roundoff error.
!>
!> NOTE: An approximation of the null vectors
!>       of A can be easily computed from R11(K)
!>       and R12(K):
!>
!>       Null( A(K) )_approx = P * ( inv(R11(K)) * R12(K) )
!>                                 (         -I           )
!>
!>
References:
[1] A Level 3 BLAS QR factorization algorithm with column pivoting developed in 1996. G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain. X. Sun, Computer Science Dept., Duke University, USA. C. H. Bischof, Math. and Comp. Sci. Div., Argonne National Lab, USA. A BLAS-3 version of the QR factorization with column pivoting. LAPACK Working Note 114 https://www.netlib.org/lapack/lawnspdf/lawn114.pdf and in SIAM J. Sci. Comput., 19(5):1486-1494, Sept. 1998. https://doi.org/10.1137/S1064827595296732

[2] A partial column norm updating strategy developed in 2006. Z. Drmac and Z. Bujanovic, Dept. of Math., University of Zagreb, Croatia. On the failure of rank revealing QR factorization software – a case study. LAPACK Working Note 176. http://www.netlib.org/lapack/lawnspdf/lawn176.pdf and in ACM Trans. Math. Softw. 35, 2, Article 12 (July 2008), 28 pages. https://doi.org/10.1145/1377612.1377616

Contributors:
!>
!>  November  2023, Igor Kozachenko, James Demmel,
!>                  EECS Department,
!>                  University of California, Berkeley, USA.
!>
!>

Definition at line 573 of file sgeqp3rk.f.

576 IMPLICIT NONE
577*
578* -- LAPACK computational routine --
579* -- LAPACK is a software package provided by Univ. of Tennessee, --
580* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
581*
582* .. Scalar Arguments ..
583 INTEGER INFO, K, KF, KMAX, LDA, LWORK, M, N, NRHS
584 REAL ABSTOL, MAXC2NRMK, RELMAXC2NRMK, RELTOL
585* ..
586* .. Array Arguments ..
587 INTEGER IWORK( * ), JPIV( * )
588 REAL A( LDA, * ), TAU( * ), WORK( * )
589* ..
590*
591* =====================================================================
592*
593* .. Parameters ..
594 INTEGER INB, INBMIN, IXOVER
595 parameter( inb = 1, inbmin = 2, ixover = 3 )
596 REAL ZERO, ONE, TWO
597 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
598* ..
599* .. Local Scalars ..
600 LOGICAL LQUERY, DONE
601 INTEGER IINFO, IOFFSET, IWS, J, JB, JBF, JMAXB, JMAX,
602 $ JMAXC2NRM, KP1, LWKOPT, MINMN, N_SUB, NB,
603 $ NBMIN, NX
604 REAL EPS, HUGEVAL, MAXC2NRM, SAFMIN
605* ..
606* .. External Subroutines ..
607 EXTERNAL slaqp2rk, slaqp3rk, xerbla
608* ..
609* .. External Functions ..
610 LOGICAL SISNAN
611 INTEGER ISAMAX, ILAENV
612 REAL SLAMCH, SNRM2, SROUNDUP_LWORK
613 EXTERNAL sisnan, slamch, snrm2, isamax, ilaenv,
614 $ sroundup_lwork
615* ..
616* .. Intrinsic Functions ..
617 INTRINSIC real, max, min
618* ..
619* .. Executable Statements ..
620*
621* Test input arguments
622* ====================
623*
624 info = 0
625 lquery = ( lwork.EQ.-1 )
626 IF( m.LT.0 ) THEN
627 info = -1
628 ELSE IF( n.LT.0 ) THEN
629 info = -2
630 ELSE IF( nrhs.LT.0 ) THEN
631 info = -3
632 ELSE IF( kmax.LT.0 ) THEN
633 info = -4
634 ELSE IF( sisnan( abstol ) ) THEN
635 info = -5
636 ELSE IF( sisnan( reltol ) ) THEN
637 info = -6
638 ELSE IF( lda.LT.max( 1, m ) ) THEN
639 info = -8
640 END IF
641*
642* If the input parameters M, N, NRHS, KMAX, LDA are valid:
643* a) Test the input workspace size LWORK for the minimum
644* size requirement IWS.
645* b) Determine the optimal block size NB and optimal
646* workspace size LWKOPT to be returned in WORK(1)
647* in case of (1) LWORK < IWS, (2) LQUERY = .TRUE.,
648* (3) when routine exits.
649* Here, IWS is the miminum workspace required for unblocked
650* code.
651*
652 IF( info.EQ.0 ) THEN
653 minmn = min( m, n )
654 IF( minmn.EQ.0 ) THEN
655 iws = 1
656 lwkopt = 1
657 ELSE
658*
659* Minimal workspace size in case of using only unblocked
660* BLAS 2 code in SLAQP2RK.
661* 1) SGEQP3RK and SLAQP2RK: 2*N to store full and partial
662* column 2-norms.
663* 2) SLAQP2RK: N+NRHS-1 to use in WORK array that is used
664* in SLARF1F subroutine inside SLAQP2RK to apply an
665* elementary reflector from the left.
666* TOTAL_WORK_SIZE = 3*N + NRHS - 1
667*
668 iws = 3*n + nrhs - 1
669*
670* Assign to NB optimal block size.
671*
672 nb = ilaenv( inb, 'SGEQP3RK', ' ', m, n, -1, -1 )
673*
674* A formula for the optimal workspace size in case of using
675* both unblocked BLAS 2 in SLAQP2RK and blocked BLAS 3 code
676* in SLAQP3RK.
677* 1) SGEQP3RK, SLAQP2RK, SLAQP3RK: 2*N to store full and
678* partial column 2-norms.
679* 2) SLAQP2RK: N+NRHS-1 to use in WORK array that is used
680* in SLARF1F subroutine to apply an elementary reflector
681* from the left.
682* 3) SLAQP3RK: NB*(N+NRHS) to use in the work array F that
683* is used to apply a block reflector from
684* the left.
685* 4) SLAQP3RK: NB to use in the auxilixary array AUX.
686* Sizes (2) and ((3) + (4)) should intersect, therefore
687* TOTAL_WORK_SIZE = 2*N + NB*( N+NRHS+1 ), given NBMIN=2.
688*
689 lwkopt = 2*n + nb*( n+nrhs+1 )
690 END IF
691 work( 1 ) = sroundup_lwork( lwkopt )
692*
693 IF( ( lwork.LT.iws ) .AND. .NOT.lquery ) THEN
694 info = -15
695 END IF
696 END IF
697*
698* NOTE: The optimal workspace size is returned in WORK(1), if
699* the input parameters M, N, NRHS, KMAX, LDA are valid.
700*
701 IF( info.NE.0 ) THEN
702 CALL xerbla( 'SGEQP3RK', -info )
703 RETURN
704 ELSE IF( lquery ) THEN
705 RETURN
706 END IF
707*
708* Quick return if possible for M=0 or N=0.
709*
710 IF( minmn.EQ.0 ) THEN
711 k = 0
712 maxc2nrmk = zero
713 relmaxc2nrmk = zero
714 work( 1 ) = sroundup_lwork( lwkopt )
715 RETURN
716 END IF
717*
718* ==================================================================
719*
720* Initialize column pivot array JPIV.
721*
722 DO j = 1, n
723 jpiv( j ) = j
724 END DO
725*
726* ==================================================================
727*
728* Initialize storage for partial and exact column 2-norms.
729* a) The elements WORK(1:N) are used to store partial column
730* 2-norms of the matrix A, and may decrease in each computation
731* step; initialize to the values of complete columns 2-norms.
732* b) The elements WORK(N+1:2*N) are used to store complete column
733* 2-norms of the matrix A, they are not changed during the
734* computation; initialize the values of complete columns 2-norms.
735*
736 DO j = 1, n
737 work( j ) = snrm2( m, a( 1, j ), 1 )
738 work( n+j ) = work( j )
739 END DO
740*
741* ==================================================================
742*
743* Compute the pivot column index and the maximum column 2-norm
744* for the whole original matrix stored in A(1:M,1:N).
745*
746 kp1 = isamax( n, work( 1 ), 1 )
747 maxc2nrm = work( kp1 )
748*
749* ==================================================================.
750*
751 IF( sisnan( maxc2nrm ) ) THEN
752*
753* Check if the matrix A contains NaN, set INFO parameter
754* to the column number where the first NaN is found and return
755* from the routine.
756*
757 k = 0
758 info = kp1
759*
760* Set MAXC2NRMK and RELMAXC2NRMK to NaN.
761*
762 maxc2nrmk = maxc2nrm
763 relmaxc2nrmk = maxc2nrm
764*
765* Array TAU is not set and contains undefined elements.
766*
767 work( 1 ) = sroundup_lwork( lwkopt )
768 RETURN
769 END IF
770*
771* ===================================================================
772*
773 IF( maxc2nrm.EQ.zero ) THEN
774*
775* Check is the matrix A is a zero matrix, set array TAU and
776* return from the routine.
777*
778 k = 0
779 maxc2nrmk = zero
780 relmaxc2nrmk = zero
781*
782 DO j = 1, minmn
783 tau( j ) = zero
784 END DO
785*
786 work( 1 ) = sroundup_lwork( lwkopt )
787 RETURN
788*
789 END IF
790*
791* ===================================================================
792*
793 hugeval = slamch( 'Overflow' )
794*
795 IF( maxc2nrm.GT.hugeval ) THEN
796*
797* Check if the matrix A contains +Inf or -Inf, set INFO parameter
798* to the column number, where the first +/-Inf is found plus N,
799* and continue the computation.
800*
801 info = n + kp1
802*
803 END IF
804*
805* ==================================================================
806*
807* Quick return if possible for the case when the first
808* stopping criterion is satisfied, i.e. KMAX = 0.
809*
810 IF( kmax.EQ.0 ) THEN
811 k = 0
812 maxc2nrmk = maxc2nrm
813 relmaxc2nrmk = one
814 DO j = 1, minmn
815 tau( j ) = zero
816 END DO
817 work( 1 ) = sroundup_lwork( lwkopt )
818 RETURN
819 END IF
820*
821* ==================================================================
822*
823 eps = slamch('Epsilon')
824*
825* Adjust ABSTOL
826*
827 IF( abstol.GE.zero ) THEN
828 safmin = slamch('Safe minimum')
829 abstol = max( abstol, two*safmin )
830 END IF
831*
832* Adjust RELTOL
833*
834 IF( reltol.GE.zero ) THEN
835 reltol = max( reltol, eps )
836 END IF
837*
838* ===================================================================
839*
840* JMAX is the maximum index of the column to be factorized,
841* which is also limited by the first stopping criterion KMAX.
842*
843 jmax = min( kmax, minmn )
844*
845* ===================================================================
846*
847* Quick return if possible for the case when the second or third
848* stopping criterion for the whole original matrix is satified,
849* i.e. MAXC2NRM <= ABSTOL or RELMAXC2NRM <= RELTOL
850* (which is ONE <= RELTOL).
851*
852 IF( maxc2nrm.LE.abstol .OR. one.LE.reltol ) THEN
853*
854 k = 0
855 maxc2nrmk = maxc2nrm
856 relmaxc2nrmk = one
857*
858 DO j = 1, minmn
859 tau( j ) = zero
860 END DO
861*
862 work( 1 ) = sroundup_lwork( lwkopt )
863 RETURN
864 END IF
865*
866* ==================================================================
867* Factorize columns
868* ==================================================================
869*
870* Determine the block size.
871*
872 nbmin = 2
873 nx = 0
874*
875 IF( ( nb.GT.1 ) .AND. ( nb.LT.minmn ) ) THEN
876*
877* Determine when to cross over from blocked to unblocked code.
878* (for N less than NX, unblocked code should be used).
879*
880 nx = max( 0, ilaenv( ixover, 'SGEQP3RK', ' ', m, n, -1,
881 $ -1 ))
882*
883 IF( nx.LT.minmn ) THEN
884*
885* Determine if workspace is large enough for blocked code.
886*
887 IF( lwork.LT.lwkopt ) THEN
888*
889* Not enough workspace to use optimal block size that
890* is currently stored in NB.
891* Reduce NB and determine the minimum value of NB.
892*
893 nb = ( lwork-2*n ) / ( n+1 )
894 nbmin = max( 2, ilaenv( inbmin, 'SGEQP3RK', ' ', m, n,
895 $ -1, -1 ) )
896*
897 END IF
898 END IF
899 END IF
900*
901* ==================================================================
902*
903* DONE is the boolean flag to rerpresent the case when the
904* factorization completed in the block factorization routine,
905* before the end of the block.
906*
907 done = .false.
908*
909* J is the column index.
910*
911 j = 1
912*
913* (1) Use blocked code initially.
914*
915* JMAXB is the maximum column index of the block, when the
916* blocked code is used, is also limited by the first stopping
917* criterion KMAX.
918*
919 jmaxb = min( kmax, minmn - nx )
920*
921 IF( nb.GE.nbmin .AND. nb.LT.jmax .AND. jmaxb.GT.0 ) THEN
922*
923* Loop over the column blocks of the matrix A(1:M,1:JMAXB). Here:
924* J is the column index of a column block;
925* JB is the column block size to pass to block factorization
926* routine in a loop step;
927* JBF is the number of columns that were actually factorized
928* that was returned by the block factorization routine
929* in a loop step, JBF <= JB;
930* N_SUB is the number of columns in the submatrix;
931* IOFFSET is the number of rows that should not be factorized.
932*
933 DO WHILE( j.LE.jmaxb )
934*
935 jb = min( nb, jmaxb-j+1 )
936 n_sub = n-j+1
937 ioffset = j-1
938*
939* Factorize JB columns among the columns A(J:N).
940*
941 CALL slaqp3rk( m, n_sub, nrhs, ioffset, jb, abstol,
942 $ reltol, kp1, maxc2nrm, a( 1, j ), lda,
943 $ done, jbf, maxc2nrmk, relmaxc2nrmk,
944 $ jpiv( j ), tau( j ),
945 $ work( j ), work( n+j ),
946 $ work( 2*n+1 ), work( 2*n+jb+1 ),
947 $ n+nrhs-j+1, iwork, iinfo )
948*
949* Set INFO on the first occurence of Inf.
950*
951 IF( iinfo.GT.n_sub .AND. info.EQ.0 ) THEN
952 info = 2*ioffset + iinfo
953 END IF
954*
955 IF( done ) THEN
956*
957* Either the submatrix is zero before the end of the
958* column block, or ABSTOL or RELTOL criterion is
959* satisfied before the end of the column block, we can
960* return from the routine. Perform the following before
961* returning:
962* a) Set the number of factorized columns K,
963* K = IOFFSET + JBF from the last call of blocked
964* routine.
965* NOTE: 1) MAXC2NRMK and RELMAXC2NRMK are returned
966* by the block factorization routine;
967* 2) The remaining TAUs are set to ZERO by the
968* block factorization routine.
969*
970 k = ioffset + jbf
971*
972* Set INFO on the first occurrence of NaN, NaN takes
973* prcedence over Inf.
974*
975 IF( iinfo.LE.n_sub .AND. iinfo.GT.0 ) THEN
976 info = ioffset + iinfo
977 END IF
978*
979* Return from the routine.
980*
981 work( 1 ) = sroundup_lwork( lwkopt )
982*
983 RETURN
984*
985 END IF
986*
987 j = j + jbf
988*
989 END DO
990*
991 END IF
992*
993* Use unblocked code to factor the last or only block.
994* J = JMAX+1 means we factorized the maximum possible number of
995* columns, that is in ELSE clause we need to compute
996* the MAXC2NORM and RELMAXC2NORM to return after we processed
997* the blocks.
998*
999 IF( j.LE.jmax ) THEN
1000*
1001* N_SUB is the number of columns in the submatrix;
1002* IOFFSET is the number of rows that should not be factorized.
1003*
1004 n_sub = n-j+1
1005 ioffset = j-1
1006*
1007 CALL slaqp2rk( m, n_sub, nrhs, ioffset, jmax-j+1,
1008 $ abstol, reltol, kp1, maxc2nrm, a( 1, j ), lda,
1009 $ kf, maxc2nrmk, relmaxc2nrmk, jpiv( j ),
1010 $ tau( j ), work( j ), work( n+j ),
1011 $ work( 2*n+1 ), iinfo )
1012*
1013* ABSTOL or RELTOL criterion is satisfied when the number of
1014* the factorized columns KF is smaller then the number
1015* of columns JMAX-J+1 supplied to be factorized by the
1016* unblocked routine, we can return from
1017* the routine. Perform the following before returning:
1018* a) Set the number of factorized columns K,
1019* b) MAXC2NRMK and RELMAXC2NRMK are returned by the
1020* unblocked factorization routine above.
1021*
1022 k = j - 1 + kf
1023*
1024* Set INFO on the first exception occurence.
1025*
1026* Set INFO on the first exception occurence of Inf or NaN,
1027* (NaN takes precedence over Inf).
1028*
1029 IF( iinfo.GT.n_sub .AND. info.EQ.0 ) THEN
1030 info = 2*ioffset + iinfo
1031 ELSE IF( iinfo.LE.n_sub .AND. iinfo.GT.0 ) THEN
1032 info = ioffset + iinfo
1033 END IF
1034*
1035 ELSE
1036*
1037* Compute the return values for blocked code.
1038*
1039* Set the number of factorized columns if the unblocked routine
1040* was not called.
1041*
1042 k = jmax
1043*
1044* If there exits a residual matrix after the blocked code:
1045* 1) compute the values of MAXC2NRMK, RELMAXC2NRMK of the
1046* residual matrix, otherwise set them to ZERO;
1047* 2) Set TAU(K+1:MINMN) to ZERO.
1048*
1049 IF( k.LT.minmn ) THEN
1050 jmaxc2nrm = k + isamax( n-k, work( k+1 ), 1 )
1051 maxc2nrmk = work( jmaxc2nrm )
1052 IF( k.EQ.0 ) THEN
1053 relmaxc2nrmk = one
1054 ELSE
1055 relmaxc2nrmk = maxc2nrmk / maxc2nrm
1056 END IF
1057*
1058 DO j = k + 1, minmn
1059 tau( j ) = zero
1060 END DO
1061*
1062 END IF
1063*
1064* END IF( J.LE.JMAX ) THEN
1065*
1066 END IF
1067*
1068 work( 1 ) = sroundup_lwork( lwkopt )
1069*
1070 RETURN
1071*
1072* End of SGEQP3RK
1073*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:57
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
slaqp2rk
subroutine slaqp2rk(m, n, nrhs, ioffset, kmax, abstol, reltol, kp1, maxc2nrm, a, lda, k, maxc2nrmk, relmaxc2nrmk, jpiv, tau, vn1, vn2, work, info)
SLAQP2RK computes truncated QR factorization with column pivoting of a real matrix block using Level ...
Definition slaqp2rk.f:334
slaqp3rk
subroutine slaqp3rk(m, n, nrhs, ioffset, nb, abstol, reltol, kp1, maxc2nrm, a, lda, done, kb, maxc2nrmk, relmaxc2nrmk, jpiv, tau, vn1, vn2, auxv, f, ldf, iwork, info)
SLAQP3RK computes a step of truncated QR factorization with column pivoting of a real m-by-n matrix A...
Definition slaqp3rk.f:392
Here is the call graph for this function:
Here is the caller graph for this function: