LAPACK 3.12.0
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.
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 "bad" 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 588 of file cgeqp3rk.f.

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