dlaqp2rk()

subroutine dlaqp2rk	(	integer	m,
		integer	n,
		integer	nrhs,
		integer	ioffset,
		integer	kmax,
		double precision	abstol,
		double precision	reltol,
		integer	kp1,
		double precision	maxc2nrm,
		double precision, dimension( lda, * )	a,
		integer	lda,
		integer	k,
		double precision	maxc2nrmk,
		double precision	relmaxc2nrmk,
		integer, dimension( * )	jpiv,
		double precision, dimension( * )	tau,
		double precision, dimension( * )	vn1,
		double precision, dimension( * )	vn2,
		double precision, dimension( * )	work,
		integer	info )

DLAQP2RK computes truncated QR factorization with column pivoting of a real matrix block using Level 2 BLAS and overwrites a real m-by-nrhs matrix B with Q**T * B.

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

Purpose:

!>
!> DLAQP2RK computes a truncated (rank K) or full rank Householder QR
!> factorization with column pivoting of a real matrix
!> block A(IOFFSET+1:M,1:N) as
!>
!>   A * P(K) = Q(K) * R(K).
!>
!> The routine uses Level 2 BLAS. The block A(1:IOFFSET,1:N)
!> is accordingly pivoted, but not factorized.
!>
!> The routine also overwrites the right-hand-sides matrix block B
!> stored in A(IOFFSET+1:M,N+1:N+NRHS) with Q(K)**T * B.
!> 

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]	IOFFSET	!> IOFFSET is INTEGER !> The number of rows of the matrix A that must be pivoted !> but not factorized. IOFFSET >= 0. !> !> IOFFSET also represents the number of columns of the whole !> original matrix A_orig that have been factorized !> in the previous steps. !>
[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-IOFFSET,N), then this stopping !> criterion is not used, factorize 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 and the arrays TAU, IPIV !> are not modified. !>
[in]	ABSTOL	!> ABSTOL is DOUBLE PRECISION, cannot be NaN. !> !> The second factorization stopping criterion. !> !> The absolute tolerance (stopping threshold) for !> maximum column 2-norm of the residual matrix. !> The algorithm converges (stops the factorization) when !> the maximum column 2-norm of the residual matrix !> is less than or equal to ABSTOL. !> !> a) 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. !> !> b) If 0.0 <= ABSTOL then the input value !> of ABSTOL is used. !>
[in]	RELTOL	!> RELTOL is DOUBLE PRECISION, cannot be NaN. !> !> The third factorization stopping criterion. !> !> The tolerance (stopping threshold) for the ratio of the !> maximum column 2-norm of the residual matrix to the maximum !> column 2-norm of the original matrix A_orig. The algorithm !> converges (stops the factorization), when this ratio is !> less than or equal to RELTOL. !> !> a) 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. !> !> d) If 0.0 <= RELTOL then the input value of RELTOL !> is used. !>
[in]	KP1	!> KP1 is INTEGER !> The index of the column with the maximum 2-norm in !> the whole original matrix A_orig determined in the !> main routine DGEQP3RK. 1 <= KP1 <= N_orig_mat. !>
[in]	MAXC2NRM	!> MAXC2NRM is DOUBLE PRECISION !> The maximum column 2-norm of the whole original !> matrix A_orig computed in the main routine DGEQP3RK. !> MAXC2NRM >= 0. !>
[in,out]	A	!> A is DOUBLE PRECISION array, dimension (LDA,N+NRHS) !> On entry: !> the M-by-N matrix A and M-by-NRHS matrix B, as in !> !> N NRHS !> array_A = M [ mat_A, mat_B ] !> !> On exit: !> 1. The elements in block A(IOFFSET+1:M,1:K) below !> the diagonal together with the array TAU represent !> the orthogonal matrix Q(K) as a product of elementary !> reflectors. !> 2. The upper triangular block of the matrix A stored !> in A(IOFFSET+1:M,1:K) is the triangular factor obtained. !> 3. The block of the matrix A stored in A(1:IOFFSET,1:N) !> has been accordingly pivoted, but not factorized. !> 4. The rest of the array A, block A(IOFFSET+1:M,K+1:N+NRHS). !> The left part A(IOFFSET+1:M,K+1:N) of this block !> contains the residual of the matrix A, and, !> if NRHS > 0, the right part of the block !> A(IOFFSET+1:M,N+1:N+NRHS) contains the block of !> the right-hand-side matrix B. Both these blocks have been !> updated by multiplication from the left by Q(K)**T. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !>
[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-IOFFSET,KMAX,N). !> !> K also represents the number of non-zero Householder !> vectors. !>
[out]	MAXC2NRMK	!> MAXC2NRMK is DOUBLE PRECISION !> The maximum column 2-norm of the residual matrix, !> when the factorization stopped at rank K. MAXC2NRMK >= 0. !>
[out]	RELMAXC2NRMK	!> RELMAXC2NRMK is DOUBLE PRECISION !> The ratio MAXC2NRMK / MAXC2NRM of the maximum column !> 2-norm of the residual matrix (when the factorization !> stopped at rank K) to the maximum column 2-norm of the !> whole original matrix A. RELMAXC2NRMK >= 0. !>
[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). !>
[out]	TAU	!> TAU is DOUBLE PRECISION array, dimension (min(M-IOFFSET,N)) !> The scalar factors of the elementary reflectors. !>
[in,out]	VN1	!> VN1 is DOUBLE PRECISION array, dimension (N) !> The vector with the partial column norms. !>
[in,out]	VN2	!> VN2 is DOUBLE PRECISION array, dimension (N) !> The vector with the exact column norms. !>
[out]	WORK	!> WORK is DOUBLE PRECISION array, dimension (N-1) !> Used in DLARF1F subroutine to apply an elementary !> reflector from the left. !>
[out]	INFO	!> INFO is INTEGER !> 1) INFO = 0: successful exit. !> 2) 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. !> 3) 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.

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 330 of file dlaqp2rk.f.

      IMPLICIT NONE
*
*  -- LAPACK auxiliary routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, IOFFSET, KP1, K, KMAX, LDA, M, N, NRHS
      DOUBLE PRECISION   ABSTOL, MAXC2NRM, MAXC2NRMK, RELMAXC2NRMK,
     $                   RELTOL
*     ..
*     .. Array Arguments ..
      INTEGER            JPIV( * )
      DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
     $                   WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ITEMP, J, JMAXC2NRM, KK, KP, MINMNFACT,
     $                   MINMNUPDT
      DOUBLE PRECISION   HUGEVAL, TEMP, TEMP2, TOL3Z
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlarf1f, dlarfg, dswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt
*     ..
*     .. External Functions ..
      LOGICAL            DISNAN
      INTEGER            IDAMAX
      DOUBLE PRECISION   DLAMCH, DNRM2
      EXTERNAL           disnan, dlamch, idamax, dnrm2
*     ..
*     .. Executable Statements ..
*
*     Initialize INFO
*
      info = 0
*
*     MINMNFACT in the smallest dimension of the submatrix
*     A(IOFFSET+1:M,1:N) to be factorized.
*
*     MINMNUPDT is the smallest dimension
*     of the subarray A(IOFFSET+1:M,1:N+NRHS) to be udated, which
*     contains the submatrices A(IOFFSET+1:M,1:N) and
*     B(IOFFSET+1:M,1:NRHS) as column blocks.
*
      minmnfact = min( m-ioffset, n )
      minmnupdt = min( m-ioffset, n+nrhs )
      kmax = min( kmax, minmnfact )
      tol3z = sqrt( dlamch( 'Epsilon' ) )
      hugeval = dlamch( 'Overflow' )
*
*     Compute the factorization, KK is the lomn loop index.
*
      DO kk = 1, kmax
*
         i = ioffset + kk
*
         IF( i.EQ.1 ) THEN
*
*           ============================================================
*
*           We are at the first column of the original whole matrix A,
*           therefore we use the computed KP1 and MAXC2NRM from the
*           main routine.
*
 
            kp = kp1
*
*           ============================================================
*
         ELSE
*
*           ============================================================
*
*           Determine the pivot column in KK-th step, i.e. the index
*           of the column with the maximum 2-norm in the
*           submatrix A(I:M,K:N).
*
            kp = ( kk-1 ) + idamax( n-kk+1, vn1( kk ), 1 )
*
*           Determine the maximum column 2-norm and the relative maximum
*           column 2-norm of the submatrix A(I:M,KK:N) in step KK.
*           RELMAXC2NRMK  will be computed later, after somecondition
*           checks on MAXC2NRMK.
*
            maxc2nrmk = vn1( kp )
*
*           ============================================================
*
*           Check if the submatrix A(I:M,KK:N) contains NaN, and set
*           INFO parameter to the column number, where the first NaN
*           is found and return from the routine.
*           We need to check the condition only if the
*           column index (same as row index) of the original whole
*           matrix is larger than 1, since the condition for whole
*           original matrix is checked in the main routine.
*
            IF( disnan( maxc2nrmk ) ) THEN
*
*              Set K, the number of factorized columns.
*              that are not zero.
*
                k = kk - 1
                info = k + kp
*
*               Set RELMAXC2NRMK to NaN.
*
                relmaxc2nrmk = maxc2nrmk
*
*               Array TAU(K+1:MINMNFACT) is not set and contains
*               undefined elements.
*
               RETURN
            END IF
*
*           ============================================================
*
*           Quick return, if the submatrix A(I:M,KK:N) is
*           a zero matrix.
*           We need to check the condition only if the
*           column index (same as row index) of the original whole
*           matrix is larger than 1, since the condition for whole
*           original matrix is checked in the main routine.
*
            IF( maxc2nrmk.EQ.zero ) THEN
*
*              Set K, the number of factorized columns.
*              that are not zero.
*
               k = kk - 1
               relmaxc2nrmk = zero
*
*              Set TAUs corresponding to the columns that were not
*              factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to ZERO.
*
               DO j = kk, minmnfact
                  tau( j ) = zero
               END DO
*
*              Return from the routine.
*
               RETURN
*
            END IF
*
*           ============================================================
*
*           Check if the submatrix A(I:M,KK:N) contains Inf,
*           set INFO parameter to the column number, where
*           the first Inf is found plus N, and continue
*           the computation.
*           We need to check the condition only if the
*           column index (same as row index) of the original whole
*           matrix is larger than 1, since the condition for whole
*           original matrix is checked in the main routine.
*
            IF( info.EQ.0 .AND. maxc2nrmk.GT.hugeval ) THEN
               info = n + kk - 1 + kp
            END IF
*
*           ============================================================
*
*           Test for the second and third stopping criteria.
*           NOTE: There is no need to test for ABSTOL >= ZERO, since
*           MAXC2NRMK is non-negative. Similarly, there is no need
*           to test for RELTOL >= ZERO, since RELMAXC2NRMK is
*           non-negative.
*           We need to check the condition only if the
*           column index (same as row index) of the original whole
*           matrix is larger than 1, since the condition for whole
*           original matrix is checked in the main routine.
 
            relmaxc2nrmk =  maxc2nrmk / maxc2nrm
*
            IF( maxc2nrmk.LE.abstol .OR. relmaxc2nrmk.LE.reltol ) THEN
*
*              Set K, the number of factorized columns.
*
               k = kk - 1
*
*              Set TAUs corresponding to the columns that were not
*              factorized to ZERO, i.e. set TAU(KK:MINMNFACT) to ZERO.
*
               DO j = kk, minmnfact
                  tau( j ) = zero
               END DO
*
*              Return from the routine.
*
               RETURN
*
            END IF
*
*           ============================================================
*
*           End ELSE of IF(I.EQ.1)
*
         END IF
*
*        ===============================================================
*
*        If the pivot column is not the first column of the
*        subblock A(1:M,KK:N):
*        1) swap the KK-th column and the KP-th pivot column
*           in A(1:M,1:N);
*        2) copy the KK-th element into the KP-th element of the partial
*           and exact 2-norm vectors VN1 and VN2. ( Swap is not needed
*           for VN1 and VN2 since we use the element with the index
*           larger than KK in the next loop step.)
*        3) Save the pivot interchange with the indices relative to the
*           the original matrix A, not the block A(1:M,1:N).
*
         IF( kp.NE.kk ) THEN
            CALL dswap( m, a( 1, kp ), 1, a( 1, kk ), 1 )
            vn1( kp ) = vn1( kk )
            vn2( kp ) = vn2( kk )
            itemp = jpiv( kp )
            jpiv( kp ) = jpiv( kk )
            jpiv( kk ) = itemp
         END IF
*
*        Generate elementary reflector H(KK) using the column A(I:M,KK),
*        if the column has more than one element, otherwise
*        the elementary reflector would be an identity matrix,
*        and TAU(KK) = ZERO.
*
         IF( i.LT.m ) THEN
            CALL dlarfg( m-i+1, a( i, kk ), a( i+1, kk ), 1,
     $                   tau( kk ) )
         ELSE
            tau( kk ) = zero
         END IF
*
*        Check if TAU(KK) contains NaN, set INFO parameter
*        to the column number where NaN is found and return from
*        the routine.
*        NOTE: There is no need to check TAU(KK) for Inf,
*        since DLARFG cannot produce TAU(KK) or Householder vector
*        below the diagonal containing Inf. Only BETA on the diagonal,
*        returned by DLARFG can contain Inf, which requires
*        TAU(KK) to contain NaN. Therefore, this case of generating Inf
*        by DLARFG is covered by checking TAU(KK) for NaN.
*
         IF( disnan( tau(kk) ) ) THEN
            k = kk - 1
            info = kk
*
*           Set MAXC2NRMK and  RELMAXC2NRMK to NaN.
*
            maxc2nrmk = tau( kk )
            relmaxc2nrmk = tau( kk )
*
*           Array TAU(KK:MINMNFACT) is not set and contains
*           undefined elements, except the first element TAU(KK) = NaN.
*
            RETURN
         END IF
*
*        Apply H(KK)**T to A(I:M,KK+1:N+NRHS) from the left.
*        ( If M >= N, then at KK = N there is no residual matrix,
*         i.e. no columns of A to update, only columns of B.
*         If M < N, then at KK = M-IOFFSET, I = M and we have a
*         one-row residual matrix in A and the elementary
*         reflector is a unit matrix, TAU(KK) = ZERO, i.e. no update
*         is needed for the residual matrix in A and the
*         right-hand-side-matrix in B.
*         Therefore, we update only if
*         KK < MINMNUPDT = min(M-IOFFSET, N+NRHS)
*         condition is satisfied, not only KK < N+NRHS )
*
         IF( kk.LT.minmnupdt ) THEN
            CALL dlarf1f( 'Left', m-i+1, n+nrhs-kk, a( i, kk ), 1,
     $                  tau( kk ), a( i, kk+1 ), lda, work( 1 ) )
         END IF
*
         IF( kk.LT.minmnfact ) THEN
*
*           Update the partial column 2-norms for the residual matrix,
*           only if the residual matrix A(I+1:M,KK+1:N) exists, i.e.
*           when KK < min(M-IOFFSET, N).
*
            DO j = kk + 1, n
               IF( vn1( j ).NE.zero ) THEN
*
*                 NOTE: The following lines follow from the analysis in
*                 Lapack Working Note 176.
*
                  temp = one - ( abs( a( i, j ) ) / vn1( j ) )**2
                  temp = max( temp, zero )
                  temp2 = temp*( vn1( j ) / vn2( j ) )**2
                  IF( temp2 .LE. tol3z ) THEN
*
*                    Compute the column 2-norm for the partial
*                    column A(I+1:M,J) by explicitly computing it,
*                    and store it in both partial 2-norm vector VN1
*                    and exact column 2-norm vector VN2.
*
                     vn1( j ) = dnrm2( m-i, a( i+1, j ), 1 )
                     vn2( j ) = vn1( j )
*
                  ELSE
*
*                    Update the column 2-norm for the partial
*                    column A(I+1:M,J) by removing one
*                    element A(I,J) and store it in partial
*                    2-norm vector VN1.
*
                     vn1( j ) = vn1( j )*sqrt( temp )
*
                  END IF
               END IF
            END DO
*
         END IF
*
*     End factorization loop
*
      END DO
*
*     If we reached this point, all colunms have been factorized,
*     i.e. no condition was triggered to exit the routine.
*     Set the number of factorized columns.
*
      k = kmax
*
*     We reached the end of the loop, i.e. all KMAX columns were
*     factorized, we need to set MAXC2NRMK and RELMAXC2NRMK before
*     we return.
*
      IF( k.LT.minmnfact ) THEN
*
         jmaxc2nrm = k + idamax( n-k, vn1( k+1 ), 1 )
         maxc2nrmk = vn1( jmaxc2nrm )
*
         IF( k.EQ.0 ) THEN
            relmaxc2nrmk = one
         ELSE
            relmaxc2nrmk = maxc2nrmk / maxc2nrm
         END IF
*
      ELSE
         maxc2nrmk = zero
         relmaxc2nrmk = zero
      END IF
*
*     We reached the end of the loop, i.e. all KMAX columns were
*     factorized, set TAUs corresponding to the columns that were
*     not factorized to ZERO, i.e. TAU(K+1:MINMNFACT) set to ZERO.
*
      DO j = k + 1, minmnfact
         tau( j ) = zero
      END DO
*
      RETURN
*
*     End of DLAQP2RK
*

Here is the call graph for this function:

Here is the caller graph for this function: