1*> \brief \b SLAQP2 computes a QR factorization with column pivoting of the matrix block.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLAQP2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqp2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqp2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqp2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE SLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
22*                          WORK )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            LDA, M, N, OFFSET
26*       ..
27*       .. Array Arguments ..
28*       INTEGER            JPVT( * )
29*       REAL               A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
30*      $                   WORK( * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> SLAQP2 computes a QR factorization with column pivoting of
40*> the block A(OFFSET+1:M,1:N).
41*> The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
42*> \endverbatim
43*
44*  Arguments:
45*  ==========
46*
47*> \param[in] M
48*> \verbatim
49*>          M is INTEGER
50*>          The number of rows of the matrix A. M >= 0.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*>          N is INTEGER
56*>          The number of columns of the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in] OFFSET
60*> \verbatim
61*>          OFFSET is INTEGER
62*>          The number of rows of the matrix A that must be pivoted
63*>          but no factorized. OFFSET >= 0.
64*> \endverbatim
65*>
66*> \param[in,out] A
67*> \verbatim
68*>          A is REAL array, dimension (LDA,N)
69*>          On entry, the M-by-N matrix A.
70*>          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
71*>          the triangular factor obtained; the elements in block
72*>          A(OFFSET+1:M,1:N) below the diagonal, together with the
73*>          array TAU, represent the orthogonal matrix Q as a product of
74*>          elementary reflectors. Block A(1:OFFSET,1:N) has been
75*>          accordingly pivoted, but no factorized.
76*> \endverbatim
77*>
78*> \param[in] LDA
79*> \verbatim
80*>          LDA is INTEGER
81*>          The leading dimension of the array A. LDA >= max(1,M).
82*> \endverbatim
83*>
84*> \param[in,out] JPVT
85*> \verbatim
86*>          JPVT is INTEGER array, dimension (N)
87*>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
88*>          to the front of A*P (a leading column); if JPVT(i) = 0,
89*>          the i-th column of A is a free column.
90*>          On exit, if JPVT(i) = k, then the i-th column of A*P
91*>          was the k-th column of A.
92*> \endverbatim
93*>
94*> \param[out] TAU
95*> \verbatim
96*>          TAU is REAL array, dimension (min(M,N))
97*>          The scalar factors of the elementary reflectors.
98*> \endverbatim
99*>
100*> \param[in,out] VN1
101*> \verbatim
102*>          VN1 is REAL array, dimension (N)
103*>          The vector with the partial column norms.
104*> \endverbatim
105*>
106*> \param[in,out] VN2
107*> \verbatim
108*>          VN2 is REAL array, dimension (N)
109*>          The vector with the exact column norms.
110*> \endverbatim
111*>
112*> \param[out] WORK
113*> \verbatim
114*>          WORK is REAL array, dimension (N)
115*> \endverbatim
116*
117*  Authors:
118*  ========
119*
120*> \author Univ. of Tennessee
121*> \author Univ. of California Berkeley
122*> \author Univ. of Colorado Denver
123*> \author NAG Ltd.
124*
125*> \ingroup realOTHERauxiliary
126*
127*> \par Contributors:
128*  ==================
129*>
130*>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
131*>    X. Sun, Computer Science Dept., Duke University, USA
132*> \n
133*>  Partial column norm updating strategy modified on April 2011
134*>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
135*>    University of Zagreb, Croatia.
136*
137*> \par References:
138*  ================
139*>
140*> LAPACK Working Note 176
141*
142*> \htmlonly
143*> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
144*> \endhtmlonly
145*
146*  =====================================================================
147      SUBROUTINE SLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
148     $                   WORK )
149*
150*  -- LAPACK auxiliary routine --
151*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
152*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153*
154*     .. Scalar Arguments ..
155      INTEGER            LDA, M, N, OFFSET
156*     ..
157*     .. Array Arguments ..
158      INTEGER            JPVT( * )
159      REAL               A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
160     $                   WORK( * )
161*     ..
162*
163*  =====================================================================
164*
165*     .. Parameters ..
166      REAL               ZERO, ONE
167      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
168*     ..
169*     .. Local Scalars ..
170      INTEGER            I, ITEMP, J, MN, OFFPI, PVT
171      REAL               AII, TEMP, TEMP2, TOL3Z
172*     ..
173*     .. External Subroutines ..
174      EXTERNAL           SLARF, SLARFG, SSWAP
175*     ..
176*     .. Intrinsic Functions ..
177      INTRINSIC          ABS, MAX, MIN, SQRT
178*     ..
179*     .. External Functions ..
180      INTEGER            ISAMAX
181      REAL               SLAMCH, SNRM2
182      EXTERNAL           ISAMAX, SLAMCH, SNRM2
183*     ..
184*     .. Executable Statements ..
185*
186      MN = MIN( M-OFFSET, N )
187      TOL3Z = SQRT(SLAMCH('Epsilon'))
188*
189*     Compute factorization.
190*
191      DO 20 I = 1, MN
192*
193         OFFPI = OFFSET + I
194*
195*        Determine ith pivot column and swap if necessary.
196*
197         PVT = ( I-1 ) + ISAMAX( N-I+1, VN1( I ), 1 )
198*
199         IF( PVT.NE.I ) THEN
200            CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
201            ITEMP = JPVT( PVT )
202            JPVT( PVT ) = JPVT( I )
203            JPVT( I ) = ITEMP
204            VN1( PVT ) = VN1( I )
205            VN2( PVT ) = VN2( I )
206         END IF
207*
208*        Generate elementary reflector H(i).
209*
210         IF( OFFPI.LT.M ) THEN
211            CALL SLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
212     $                   TAU( I ) )
213         ELSE
214            CALL SLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
215         END IF
216*
217         IF( I.LT.N ) THEN
218*
219*           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
220*
221            AII = A( OFFPI, I )
222            A( OFFPI, I ) = ONE
223            CALL SLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
224     $                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
225            A( OFFPI, I ) = AII
226         END IF
227*
228*        Update partial column norms.
229*
230         DO 10 J = I + 1, N
231            IF( VN1( J ).NE.ZERO ) THEN
232*
233*              NOTE: The following 4 lines follow from the analysis in
234*              Lapack Working Note 176.
235*
236               TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
237               TEMP = MAX( TEMP, ZERO )
238               TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
239               IF( TEMP2 .LE. TOL3Z ) THEN
240                  IF( OFFPI.LT.M ) THEN
241                     VN1( J ) = SNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
242                     VN2( J ) = VN1( J )
243                  ELSE
244                     VN1( J ) = ZERO
245                     VN2( J ) = ZERO
246                  END IF
247               ELSE
248                  VN1( J ) = VN1( J )*SQRT( TEMP )
249               END IF
250            END IF
251   10    CONTINUE
252*
253   20 CONTINUE
254*
255      RETURN
256*
257*     End of SLAQP2
258*
259      END
260