SRC/slaqp2.f

   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 *> \date September 2012
 126 *
 127 *> \ingroup realOTHERauxiliary
 128 *
 129 *> \par Contributors:
 130 *  ==================
 131 *>
 132 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
 133 *>    X. Sun, Computer Science Dept., Duke University, USA
 134 *> \n
 135 *>  Partial column norm updating strategy modified on April 2011
 136 *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
 137 *>    University of Zagreb, Croatia.
 138 *
 139 *> \par References:
 140 *  ================
 141 *>
 142 *> LAPACK Working Note 176
 143 *
 144 *> \htmlonly
 145 *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a>
 146 *> \endhtmlonly
 147 *
 148 *  =====================================================================
 149       SUBROUTINE SLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
 150      $                   WORK )
 151 *
 152 *  -- LAPACK auxiliary routine (version 3.4.2) --
 153 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 154 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 155 *     September 2012
 156 *
 157 *     .. Scalar Arguments ..
 158       INTEGER            LDA, M, N, OFFSET
 159 *     ..
 160 *     .. Array Arguments ..
 161       INTEGER            JPVT( * )
 162       REAL               A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
 163      $                   WORK( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       REAL               ZERO, ONE
 170       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 171 *     ..
 172 *     .. Local Scalars ..
 173       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
 174       REAL               AII, TEMP, TEMP2, TOL3Z
 175 *     ..
 176 *     .. External Subroutines ..
 177       EXTERNAL           SLARF, SLARFG, SSWAP
 178 *     ..
 179 *     .. Intrinsic Functions ..
 180       INTRINSIC          ABS, MAX, MIN, SQRT
 181 *     ..
 182 *     .. External Functions ..
 183       INTEGER            ISAMAX
 184       REAL               SLAMCH, SNRM2
 185       EXTERNAL           ISAMAX, SLAMCH, SNRM2
 186 *     ..
 187 *     .. Executable Statements ..
 188 *
 189       MN = MIN( M-OFFSET, N )
 190       TOL3Z = SQRT(SLAMCH('Epsilon'))
 191 *
 192 *     Compute factorization.
 193 *
 194       DO 20 I = 1, MN
 195 *
 196          OFFPI = OFFSET + I
 197 *
 198 *        Determine ith pivot column and swap if necessary.
 199 *
 200          PVT = ( I-1 ) + ISAMAX( N-I+1, VN1( I ), 1 )
 201 *
 202          IF( PVT.NE.I ) THEN
 203             CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
 204             ITEMP = JPVT( PVT )
 205             JPVT( PVT ) = JPVT( I )
 206             JPVT( I ) = ITEMP
 207             VN1( PVT ) = VN1( I )
 208             VN2( PVT ) = VN2( I )
 209          END IF
 210 *
 211 *        Generate elementary reflector H(i).
 212 *
 213          IF( OFFPI.LT.M ) THEN
 214             CALL SLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
 215      $                   TAU( I ) )
 216          ELSE
 217             CALL SLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
 218          END IF
 219 *
 220          IF( I.LT.N ) THEN
 221 *
 222 *           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
 223 *
 224             AII = A( OFFPI, I )
 225             A( OFFPI, I ) = ONE
 226             CALL SLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
 227      $                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
 228             A( OFFPI, I ) = AII
 229          END IF
 230 *
 231 *        Update partial column norms.
 232 *
 233          DO 10 J = I + 1, N
 234             IF( VN1( J ).NE.ZERO ) THEN
 235 *
 236 *              NOTE: The following 4 lines follow from the analysis in
 237 *              Lapack Working Note 176.
 238 *
 239                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
 240                TEMP = MAX( TEMP, ZERO )
 241                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
 242                IF( TEMP2 .LE. TOL3Z ) THEN
 243                   IF( OFFPI.LT.M ) THEN
 244                      VN1( J ) = SNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
 245                      VN2( J ) = VN1( J )
 246                   ELSE
 247                      VN1( J ) = ZERO
 248                      VN2( J ) = ZERO
 249                   END IF
 250                ELSE
 251                   VN1( J ) = VN1( J )*SQRT( TEMP )
 252                END IF
 253             END IF
 254    10    CONTINUE
 255 *
 256    20 CONTINUE
 257 *
 258       RETURN
 259 *
 260 *     End of SLAQP2
 261 *
 262       END