dlaqp2.f

Go to the documentation of this file.

*> \brief \b DLAQP2 computes a QR factorization with column pivoting of the matrix block.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAQP2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqp2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqp2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqp2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
*                          WORK )
*
*       .. Scalar Arguments ..
*       INTEGER            LDA, M, N, OFFSET
*       ..
*       .. Array Arguments ..
*       INTEGER            JPVT( * )
*       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
*      $                   WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DLAQP2 computes a QR factorization with column pivoting of
*> the block A(OFFSET+1:M,1:N).
*> The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrix A. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] OFFSET
*> \verbatim
*>          OFFSET is INTEGER
*>          The number of rows of the matrix A that must be pivoted
*>          but no factorized. OFFSET >= 0.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*>          On entry, the M-by-N matrix A.
*>          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
*>          the triangular factor obtained; the elements in block
*>          A(OFFSET+1:M,1:N) below the diagonal, together with the
*>          array TAU, represent the orthogonal matrix Q as a product of
*>          elementary reflectors. Block A(1:OFFSET,1:N) has been
*>          accordingly pivoted, but no factorized.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[in,out] JPVT
*> \verbatim
*>          JPVT is INTEGER array, dimension (N)
*>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
*>          to the front of A*P (a leading column); if JPVT(i) = 0,
*>          the i-th column of A is a free column.
*>          On exit, if JPVT(i) = k, then the i-th column of A*P
*>          was the k-th column of A.
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
*>          The scalar factors of the elementary reflectors.
*> \endverbatim
*>
*> \param[in,out] VN1
*> \verbatim
*>          VN1 is DOUBLE PRECISION array, dimension (N)
*>          The vector with the partial column norms.
*> \endverbatim
*>
*> \param[in,out] VN2
*> \verbatim
*>          VN2 is DOUBLE PRECISION array, dimension (N)
*>          The vector with the exact column norms.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (N)
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup laqp2
*
*> \par Contributors:
*  ==================
*>
*>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
*>    X. Sun, Computer Science Dept., Duke University, USA
*> \n
*>  Partial column norm updating strategy modified on April 2011
*>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
*>    University of Zagreb, Croatia.
*
*> \par References:
*  ================
*>
*> LAPACK Working Note 176
*
*> \htmlonly
*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
*> \endhtmlonly
*
*  =====================================================================
      SUBROUTINE dlaqp2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
     $                   WORK )
*
*  -- 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            LDA, M, N, OFFSET
*     ..
*     .. Array Arguments ..
      INTEGER            JPVT( * )
      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, MN, OFFPI, PVT
      DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlarf, dlarfg, dswap
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, min, sqrt
*     ..
*     .. External Functions ..
      INTEGER            IDAMAX
      DOUBLE PRECISION   DLAMCH, DNRM2
      EXTERNAL           idamax, dlamch, dnrm2
*     ..
*     .. Executable Statements ..
*
      mn = min( m-offset, n )
      tol3z = sqrt(dlamch('Epsilon'))
*
*     Compute factorization.
*
      DO 20 i = 1, mn
*
         offpi = offset + i
*
*        Determine ith pivot column and swap if necessary.
*
         pvt = ( i-1 ) + idamax( n-i+1, vn1( i ), 1 )
*
         IF( pvt.NE.i ) THEN
            CALL dswap( m, a( 1, pvt ), 1, a( 1, i ), 1 )
            itemp = jpvt( pvt )
            jpvt( pvt ) = jpvt( i )
            jpvt( i ) = itemp
            vn1( pvt ) = vn1( i )
            vn2( pvt ) = vn2( i )
         END IF
*
*        Generate elementary reflector H(i).
*
         IF( offpi.LT.m ) THEN
            CALL dlarfg( m-offpi+1, a( offpi, i ), a( offpi+1, i ), 1,
     $                   tau( i ) )
         ELSE
            CALL dlarfg( 1, a( m, i ), a( m, i ), 1, tau( i ) )
         END IF
*
         IF( i.LT.n ) THEN
*
*           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
*
            aii = a( offpi, i )
            a( offpi, i ) = one
            CALL dlarf( 'Left', m-offpi+1, n-i, a( offpi, i ), 1,
     $                  tau( i ), a( offpi, i+1 ), lda, work( 1 ) )
            a( offpi, i ) = aii
         END IF
*
*        Update partial column norms.
*
         DO 10 j = i + 1, n
            IF( vn1( j ).NE.zero ) THEN
*
*              NOTE: The following 4 lines follow from the analysis in
*              Lapack Working Note 176.
*
               temp = one - ( abs( a( offpi, j ) ) / vn1( j ) )**2
               temp = max( temp, zero )
               temp2 = temp*( vn1( j ) / vn2( j ) )**2
               IF( temp2 .LE. tol3z ) THEN
                  IF( offpi.LT.m ) THEN
                     vn1( j ) = dnrm2( m-offpi, a( offpi+1, j ), 1 )
                     vn2( j ) = vn1( j )
                  ELSE
                     vn1( j ) = zero
                     vn2( j ) = zero
                  END IF
               ELSE
                  vn1( j ) = vn1( j )*sqrt( temp )
               END IF
            END IF
   10    CONTINUE
*
   20 CONTINUE
*
      RETURN
*
*     End of DLAQP2
*
      END