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

◆ cgeqp3rk()

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

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

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