001:       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
002:      $                   WORK )
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            LDA, M, N, OFFSET
010: *     ..
011: *     .. Array Arguments ..
012:       INTEGER            JPVT( * )
013:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
014:      $                   WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DLAQP2 computes a QR factorization with column pivoting of
021: *  the block A(OFFSET+1:M,1:N).
022: *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  M       (input) INTEGER
028: *          The number of rows of the matrix A. M >= 0.
029: *
030: *  N       (input) INTEGER
031: *          The number of columns of the matrix A. N >= 0.
032: *
033: *  OFFSET  (input) INTEGER
034: *          The number of rows of the matrix A that must be pivoted
035: *          but no factorized. OFFSET >= 0.
036: *
037: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
038: *          On entry, the M-by-N matrix A.
039: *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is 
040: *          the triangular factor obtained; the elements in block
041: *          A(OFFSET+1:M,1:N) below the diagonal, together with the
042: *          array TAU, represent the orthogonal matrix Q as a product of
043: *          elementary reflectors. Block A(1:OFFSET,1:N) has been
044: *          accordingly pivoted, but no factorized.
045: *
046: *  LDA     (input) INTEGER
047: *          The leading dimension of the array A. LDA >= max(1,M).
048: *
049: *  JPVT    (input/output) INTEGER array, dimension (N)
050: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
051: *          to the front of A*P (a leading column); if JPVT(i) = 0,
052: *          the i-th column of A is a free column.
053: *          On exit, if JPVT(i) = k, then the i-th column of A*P
054: *          was the k-th column of A.
055: *
056: *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
057: *          The scalar factors of the elementary reflectors.
058: *
059: *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)
060: *          The vector with the partial column norms.
061: *
062: *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)
063: *          The vector with the exact column norms.
064: *
065: *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
066: *
067: *  Further Details
068: *  ===============
069: *
070: *  Based on contributions by
071: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
072: *    X. Sun, Computer Science Dept., Duke University, USA
073: *
074: *  Partial column norm updating strategy modified by
075: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
076: *    University of Zagreb, Croatia.
077: *    June 2006.
078: *  For more details see LAPACK Working Note 176.
079: *  =====================================================================
080: *
081: *     .. Parameters ..
082:       DOUBLE PRECISION   ZERO, ONE
083:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
084: *     ..
085: *     .. Local Scalars ..
086:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
087:       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
088: *     ..
089: *     .. External Subroutines ..
090:       EXTERNAL           DLARF, DLARFP, DSWAP
091: *     ..
092: *     .. Intrinsic Functions ..
093:       INTRINSIC          ABS, MAX, MIN, SQRT
094: *     ..
095: *     .. External Functions ..
096:       INTEGER            IDAMAX
097:       DOUBLE PRECISION   DLAMCH, DNRM2
098:       EXTERNAL           IDAMAX, DLAMCH, DNRM2
099: *     ..
100: *     .. Executable Statements ..
101: *
102:       MN = MIN( M-OFFSET, N )
103:       TOL3Z = SQRT(DLAMCH('Epsilon'))
104: *
105: *     Compute factorization.
106: *
107:       DO 20 I = 1, MN
108: *
109:          OFFPI = OFFSET + I
110: *
111: *        Determine ith pivot column and swap if necessary.
112: *
113:          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
114: *
115:          IF( PVT.NE.I ) THEN
116:             CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
117:             ITEMP = JPVT( PVT )
118:             JPVT( PVT ) = JPVT( I )
119:             JPVT( I ) = ITEMP
120:             VN1( PVT ) = VN1( I )
121:             VN2( PVT ) = VN2( I )
122:          END IF
123: *
124: *        Generate elementary reflector H(i).
125: *
126:          IF( OFFPI.LT.M ) THEN
127:             CALL DLARFP( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
128:      $                   TAU( I ) )
129:          ELSE
130:             CALL DLARFP( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
131:          END IF
132: *
133:          IF( I.LE.N ) THEN
134: *
135: *           Apply H(i)' to A(offset+i:m,i+1:n) from the left.
136: *
137:             AII = A( OFFPI, I )
138:             A( OFFPI, I ) = ONE
139:             CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
140:      $                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
141:             A( OFFPI, I ) = AII
142:          END IF
143: *
144: *        Update partial column norms.
145: *
146:          DO 10 J = I + 1, N
147:             IF( VN1( J ).NE.ZERO ) THEN
148: *
149: *              NOTE: The following 4 lines follow from the analysis in
150: *              Lapack Working Note 176.
151: *
152:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
153:                TEMP = MAX( TEMP, ZERO )
154:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
155:                IF( TEMP2 .LE. TOL3Z ) THEN
156:                   IF( OFFPI.LT.M ) THEN
157:                      VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
158:                      VN2( J ) = VN1( J )
159:                   ELSE
160:                      VN1( J ) = ZERO
161:                      VN2( J ) = ZERO
162:                   END IF
163:                ELSE
164:                   VN1( J ) = VN1( J )*SQRT( TEMP )
165:                END IF
166:             END IF
167:    10    CONTINUE
168: *
169:    20 CONTINUE
170: *
171:       RETURN
172: *
173: *     End of DLAQP2
174: *
175:       END
176: