SRC/claqp2.f

   1 *> \brief \b CLAQP2 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 CLAQP2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claqp2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claqp2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claqp2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLAQP2( 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               VN1( * ), VN2( * )
  30 *       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> CLAQP2 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 complexOTHERauxiliary
 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 CLAQP2( 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               VN1( * ), VN2( * )
 163       COMPLEX            A( LDA, * ), TAU( * ), WORK( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       REAL               ZERO, ONE
 170       COMPLEX            CONE
 171       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0,
 172      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
 173 *     ..
 174 *     .. Local Scalars ..
 175       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
 176       REAL               TEMP, TEMP2, TOL3Z
 177       COMPLEX            AII
 178 *     ..
 179 *     .. External Subroutines ..
 180       EXTERNAL           CLARF, CLARFG, CSWAP
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          ABS, CONJG, MAX, MIN, SQRT
 184 *     ..
 185 *     .. External Functions ..
 186       INTEGER            ISAMAX
 187       REAL               SCNRM2, SLAMCH
 188       EXTERNAL           ISAMAX, SCNRM2, SLAMCH
 189 *     ..
 190 *     .. Executable Statements ..
 191 *
 192       MN = MIN( M-OFFSET, N )
 193       TOL3Z = SQRT(SLAMCH('Epsilon'))
 194 *
 195 *     Compute factorization.
 196 *
 197       DO 20 I = 1, MN
 198 *
 199          OFFPI = OFFSET + I
 200 *
 201 *        Determine ith pivot column and swap if necessary.
 202 *
 203          PVT = ( I-1 ) + ISAMAX( N-I+1, VN1( I ), 1 )
 204 *
 205          IF( PVT.NE.I ) THEN
 206             CALL CSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
 207             ITEMP = JPVT( PVT )
 208             JPVT( PVT ) = JPVT( I )
 209             JPVT( I ) = ITEMP
 210             VN1( PVT ) = VN1( I )
 211             VN2( PVT ) = VN2( I )
 212          END IF
 213 *
 214 *        Generate elementary reflector H(i).
 215 *
 216          IF( OFFPI.LT.M ) THEN
 217             CALL CLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
 218      $                   TAU( I ) )
 219          ELSE
 220             CALL CLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
 221          END IF
 222 *
 223          IF( I.LT.N ) THEN
 224 *
 225 *           Apply H(i)**H to A(offset+i:m,i+1:n) from the left.
 226 *
 227             AII = A( OFFPI, I )
 228             A( OFFPI, I ) = CONE
 229             CALL CLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
 230      $                  CONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
 231      $                  WORK( 1 ) )
 232             A( OFFPI, I ) = AII
 233          END IF
 234 *
 235 *        Update partial column norms.
 236 *
 237          DO 10 J = I + 1, N
 238             IF( VN1( J ).NE.ZERO ) THEN
 239 *
 240 *              NOTE: The following 4 lines follow from the analysis in
 241 *              Lapack Working Note 176.
 242 *
 243                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
 244                TEMP = MAX( TEMP, ZERO )
 245                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
 246                IF( TEMP2 .LE. TOL3Z ) THEN
 247                   IF( OFFPI.LT.M ) THEN
 248                      VN1( J ) = SCNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
 249                      VN2( J ) = VN1( J )
 250                   ELSE
 251                      VN1( J ) = ZERO
 252                      VN2( J ) = ZERO
 253                   END IF
 254                ELSE
 255                   VN1( J ) = VN1( J )*SQRT( TEMP )
 256                END IF
 257             END IF
 258    10    CONTINUE
 259 *
 260    20 CONTINUE
 261 *
 262       RETURN
 263 *
 264 *     End of CLAQP2
 265 *
 266       END