Annotation of rpl/lapack/lapack/zlaqp2.f, revision 1.9

1.1       bertrand    1:       SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
                      2:      $                   WORK )
                      3: *
1.9     ! bertrand    4: *  -- LAPACK auxiliary routine (version 3.3.1) --
1.1       bertrand    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9     ! bertrand    7: *  -- April 2011                                                      --
1.1       bertrand    8: *
                      9: *     .. Scalar Arguments ..
                     10:       INTEGER            LDA, M, N, OFFSET
                     11: *     ..
                     12: *     .. Array Arguments ..
                     13:       INTEGER            JPVT( * )
                     14:       DOUBLE PRECISION   VN1( * ), VN2( * )
                     15:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
                     16: *     ..
                     17: *
                     18: *  Purpose
                     19: *  =======
                     20: *
                     21: *  ZLAQP2 computes a QR factorization with column pivoting of
                     22: *  the block A(OFFSET+1:M,1:N).
                     23: *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
                     24: *
                     25: *  Arguments
                     26: *  =========
                     27: *
                     28: *  M       (input) INTEGER
                     29: *          The number of rows of the matrix A. M >= 0.
                     30: *
                     31: *  N       (input) INTEGER
                     32: *          The number of columns of the matrix A. N >= 0.
                     33: *
                     34: *  OFFSET  (input) INTEGER
                     35: *          The number of rows of the matrix A that must be pivoted
                     36: *          but no factorized. OFFSET >= 0.
                     37: *
                     38: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
                     39: *          On entry, the M-by-N matrix A.
                     40: *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
                     41: *          the triangular factor obtained; the elements in block
                     42: *          A(OFFSET+1:M,1:N) below the diagonal, together with the
                     43: *          array TAU, represent the orthogonal matrix Q as a product of
                     44: *          elementary reflectors. Block A(1:OFFSET,1:N) has been
                     45: *          accordingly pivoted, but no factorized.
                     46: *
                     47: *  LDA     (input) INTEGER
                     48: *          The leading dimension of the array A. LDA >= max(1,M).
                     49: *
                     50: *  JPVT    (input/output) INTEGER array, dimension (N)
                     51: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
                     52: *          to the front of A*P (a leading column); if JPVT(i) = 0,
                     53: *          the i-th column of A is a free column.
                     54: *          On exit, if JPVT(i) = k, then the i-th column of A*P
                     55: *          was the k-th column of A.
                     56: *
                     57: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
                     58: *          The scalar factors of the elementary reflectors.
                     59: *
                     60: *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)
                     61: *          The vector with the partial column norms.
                     62: *
                     63: *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)
                     64: *          The vector with the exact column norms.
                     65: *
                     66: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
                     67: *
                     68: *  Further Details
                     69: *  ===============
                     70: *
                     71: *  Based on contributions by
                     72: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
                     73: *    X. Sun, Computer Science Dept., Duke University, USA
                     74: *
                     75: *  Partial column norm updating strategy modified by
                     76: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
                     77: *    University of Zagreb, Croatia.
1.9     ! bertrand   78: *  -- April 2011                                                      --
1.1       bertrand   79: *  For more details see LAPACK Working Note 176.
                     80: *  =====================================================================
                     81: *
                     82: *     .. Parameters ..
                     83:       DOUBLE PRECISION   ZERO, ONE
                     84:       COMPLEX*16         CONE
                     85:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
                     86:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
                     87: *     ..
                     88: *     .. Local Scalars ..
                     89:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
                     90:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
                     91:       COMPLEX*16         AII
                     92: *     ..
                     93: *     .. External Subroutines ..
1.5       bertrand   94:       EXTERNAL           ZLARF, ZLARFG, ZSWAP
1.1       bertrand   95: *     ..
                     96: *     .. Intrinsic Functions ..
                     97:       INTRINSIC          ABS, DCONJG, MAX, MIN, SQRT
                     98: *     ..
                     99: *     .. External Functions ..
                    100:       INTEGER            IDAMAX
                    101:       DOUBLE PRECISION   DLAMCH, DZNRM2
                    102:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
                    103: *     ..
                    104: *     .. Executable Statements ..
                    105: *
                    106:       MN = MIN( M-OFFSET, N )
                    107:       TOL3Z = SQRT(DLAMCH('Epsilon'))
                    108: *
                    109: *     Compute factorization.
                    110: *
                    111:       DO 20 I = 1, MN
                    112: *
                    113:          OFFPI = OFFSET + I
                    114: *
                    115: *        Determine ith pivot column and swap if necessary.
                    116: *
                    117:          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
                    118: *
                    119:          IF( PVT.NE.I ) THEN
                    120:             CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
                    121:             ITEMP = JPVT( PVT )
                    122:             JPVT( PVT ) = JPVT( I )
                    123:             JPVT( I ) = ITEMP
                    124:             VN1( PVT ) = VN1( I )
                    125:             VN2( PVT ) = VN2( I )
                    126:          END IF
                    127: *
                    128: *        Generate elementary reflector H(i).
                    129: *
                    130:          IF( OFFPI.LT.M ) THEN
1.5       bertrand  131:             CALL ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
1.1       bertrand  132:      $                   TAU( I ) )
                    133:          ELSE
1.5       bertrand  134:             CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
1.1       bertrand  135:          END IF
                    136: *
                    137:          IF( I.LT.N ) THEN
                    138: *
1.9     ! bertrand  139: *           Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
1.1       bertrand  140: *
                    141:             AII = A( OFFPI, I )
                    142:             A( OFFPI, I ) = CONE
                    143:             CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
                    144:      $                  DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
                    145:      $                  WORK( 1 ) )
                    146:             A( OFFPI, I ) = AII
                    147:          END IF
                    148: *
                    149: *        Update partial column norms.
                    150: *
                    151:          DO 10 J = I + 1, N
                    152:             IF( VN1( J ).NE.ZERO ) THEN
                    153: *
                    154: *              NOTE: The following 4 lines follow from the analysis in
                    155: *              Lapack Working Note 176.
                    156: *
                    157:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
                    158:                TEMP = MAX( TEMP, ZERO )
                    159:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
                    160:                IF( TEMP2 .LE. TOL3Z ) THEN
                    161:                   IF( OFFPI.LT.M ) THEN
                    162:                      VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
                    163:                      VN2( J ) = VN1( J )
                    164:                   ELSE
                    165:                      VN1( J ) = ZERO
                    166:                      VN2( J ) = ZERO
                    167:                   END IF
                    168:                ELSE
                    169:                   VN1( J ) = VN1( J )*SQRT( TEMP )
                    170:                END IF
                    171:             END IF
                    172:    10    CONTINUE
                    173: *
                    174:    20 CONTINUE
                    175: *
                    176:       RETURN
                    177: *
                    178: *     End of ZLAQP2
                    179: *
                    180:       END

CVSweb interface <joel.bertrand@systella.fr>