LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgeqp3rk()

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

DGEQP3RK 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 DGEQP3RK + dependencies [TGZ] [ZIP] [TXT]

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