001:       SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
002:      $                   WORK )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            LDA, M, N, OFFSET
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            JPVT( * )
014:       DOUBLE PRECISION   VN1( * ), VN2( * )
015:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZLAQP2 computes a QR factorization with column pivoting of
022: *  the block A(OFFSET+1:M,1:N).
023: *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  M       (input) INTEGER
029: *          The number of rows of the matrix A. M >= 0.
030: *
031: *  N       (input) INTEGER
032: *          The number of columns of the matrix A. N >= 0.
033: *
034: *  OFFSET  (input) INTEGER
035: *          The number of rows of the matrix A that must be pivoted
036: *          but no factorized. OFFSET >= 0.
037: *
038: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
039: *          On entry, the M-by-N matrix A.
040: *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
041: *          the triangular factor obtained; the elements in block
042: *          A(OFFSET+1:M,1:N) below the diagonal, together with the
043: *          array TAU, represent the orthogonal matrix Q as a product of
044: *          elementary reflectors. Block A(1:OFFSET,1:N) has been
045: *          accordingly pivoted, but no factorized.
046: *
047: *  LDA     (input) INTEGER
048: *          The leading dimension of the array A. LDA >= max(1,M).
049: *
050: *  JPVT    (input/output) INTEGER array, dimension (N)
051: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
052: *          to the front of A*P (a leading column); if JPVT(i) = 0,
053: *          the i-th column of A is a free column.
054: *          On exit, if JPVT(i) = k, then the i-th column of A*P
055: *          was the k-th column of A.
056: *
057: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
058: *          The scalar factors of the elementary reflectors.
059: *
060: *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)
061: *          The vector with the partial column norms.
062: *
063: *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)
064: *          The vector with the exact column norms.
065: *
066: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
067: *
068: *  Further Details
069: *  ===============
070: *
071: *  Based on contributions by
072: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
073: *    X. Sun, Computer Science Dept., Duke University, USA
074: *
075: *  Partial column norm updating strategy modified by
076: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
077: *    University of Zagreb, Croatia.
078: *    June 2006.
079: *  For more details see LAPACK Working Note 176.
080: *  =====================================================================
081: *
082: *     .. Parameters ..
083:       DOUBLE PRECISION   ZERO, ONE
084:       COMPLEX*16         CONE
085:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
086:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
087: *     ..
088: *     .. Local Scalars ..
089:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
090:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
091:       COMPLEX*16         AII
092: *     ..
093: *     .. External Subroutines ..
094:       EXTERNAL           ZLARF, ZLARFP, ZSWAP
095: *     ..
096: *     .. Intrinsic Functions ..
097:       INTRINSIC          ABS, DCONJG, MAX, MIN, SQRT
098: *     ..
099: *     .. 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
131:             CALL ZLARFP( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
132:      $                   TAU( I ) )
133:          ELSE
134:             CALL ZLARFP( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
135:          END IF
136: *
137:          IF( I.LT.N ) THEN
138: *
139: *           Apply H(i)' to A(offset+i:m,i+1:n) from the left.
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
181: