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

◆ zgeqp3rk()

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

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

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

Purpose:
!>
!> ZGEQP3RK performs two tasks simultaneously:
!>
!> Task 1: The routine computes a truncated (rank K) or full rank
!> Householder QR factorization with column pivoting of a complex
!> 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 complex M-by-NRHS
!> matrix B with  Q(K)**H * 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)**H * B, the matrix B with the unitary
!>        transformation Q(K)**H 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 unitary 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**H,
!>
!> where 1 <= j <= K and
!>   I    is an M-by-M identity matrix,
!>   tau  is a complex scalar,
!>   v    is a complex 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( 10*max(M,N)*EPS, sqrt(EPS) ) <= RELTOL
!>
[in,out]A
!>          A is COMPLEX*16 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
!>                 unitary 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)**H * 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 COMPLEX*16 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 COMPLEX*16 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 >= N+NRHS-1, otherwise.
!>          For optimal performance LWORK >= NB*( N+NRHS+1 ),
!>          where NB is the optimal block size for ZGEQP3RK 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]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (2*N)
!>
[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 ZLAQP3RK ).
!>
[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:
!> ZGEQP3RK is based on the same BLAS3 Householder QR factorization
!> algorithm with column pivoting as in ZGEQP3 routine which uses
!> ZLARFG 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 579 of file zgeqp3rk.f.

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