SRC/stzrzf.f

   1 *> \brief \b STZRZF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download STZRZF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stzrzf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stzrzf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stzrzf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, LDA, LWORK, M, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       REAL               A( LDA, * ), TAU( * ), WORK( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
  37 *> to upper triangular form by means of orthogonal transformations.
  38 *>
  39 *> The upper trapezoidal matrix A is factored as
  40 *>
  41 *>    A = ( R  0 ) * Z,
  42 *>
  43 *> where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
  44 *> triangular matrix.
  45 *> \endverbatim
  46 *
  47 *  Arguments:
  48 *  ==========
  49 *
  50 *> \param[in] M
  51 *> \verbatim
  52 *>          M is INTEGER
  53 *>          The number of rows of the matrix A.  M >= 0.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] N
  57 *> \verbatim
  58 *>          N is INTEGER
  59 *>          The number of columns of the matrix A.  N >= M.
  60 *> \endverbatim
  61 *>
  62 *> \param[in,out] A
  63 *> \verbatim
  64 *>          A is REAL array, dimension (LDA,N)
  65 *>          On entry, the leading M-by-N upper trapezoidal part of the
  66 *>          array A must contain the matrix to be factorized.
  67 *>          On exit, the leading M-by-M upper triangular part of A
  68 *>          contains the upper triangular matrix R, and elements M+1 to
  69 *>          N of the first M rows of A, with the array TAU, represent the
  70 *>          orthogonal matrix Z as a product of M elementary reflectors.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] LDA
  74 *> \verbatim
  75 *>          LDA is INTEGER
  76 *>          The leading dimension of the array A.  LDA >= max(1,M).
  77 *> \endverbatim
  78 *>
  79 *> \param[out] TAU
  80 *> \verbatim
  81 *>          TAU is REAL array, dimension (M)
  82 *>          The scalar factors of the elementary reflectors.
  83 *> \endverbatim
  84 *>
  85 *> \param[out] WORK
  86 *> \verbatim
  87 *>          WORK is REAL array, dimension (MAX(1,LWORK))
  88 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  89 *> \endverbatim
  90 *>
  91 *> \param[in] LWORK
  92 *> \verbatim
  93 *>          LWORK is INTEGER
  94 *>          The dimension of the array WORK.  LWORK >= max(1,M).
  95 *>          For optimum performance LWORK >= M*NB, where NB is
  96 *>          the optimal blocksize.
  97 *>
  98 *>          If LWORK = -1, then a workspace query is assumed; the routine
  99 *>          only calculates the optimal size of the WORK array, returns
 100 *>          this value as the first entry of the WORK array, and no error
 101 *>          message related to LWORK is issued by XERBLA.
 102 *> \endverbatim
 103 *>
 104 *> \param[out] INFO
 105 *> \verbatim
 106 *>          INFO is INTEGER
 107 *>          = 0:  successful exit
 108 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 109 *> \endverbatim
 110 *
 111 *  Authors:
 112 *  ========
 113 *
 114 *> \author Univ. of Tennessee
 115 *> \author Univ. of California Berkeley
 116 *> \author Univ. of Colorado Denver
 117 *> \author NAG Ltd.
 118 *
 119 *> \date April 2012
 120 *
 121 *> \ingroup realOTHERcomputational
 122 *
 123 *> \par Contributors:
 124 *  ==================
 125 *>
 126 *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 127 *
 128 *> \par Further Details:
 129 *  =====================
 130 *>
 131 *> \verbatim
 132 *>
 133 *>  The N-by-N matrix Z can be computed by
 134 *>
 135 *>     Z =  Z(1)*Z(2)* ... *Z(M)
 136 *>
 137 *>  where each N-by-N Z(k) is given by
 138 *>
 139 *>     Z(k) = I - tau(k)*v(k)*v(k)**T
 140 *>
 141 *>  with v(k) is the kth row vector of the M-by-N matrix
 142 *>
 143 *>     V = ( I   A(:,M+1:N) )
 144 *>
 145 *>  I is the M-by-M identity matrix, A(:,M+1:N)
 146 *>  is the output stored in A on exit from DTZRZF,
 147 *>  and tau(k) is the kth element of the array TAU.
 148 *>
 149 *> \endverbatim
 150 *>
 151 *  =====================================================================
 152       SUBROUTINE STZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 153 *
 154 *  -- LAPACK computational routine (version 3.4.1) --
 155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 157 *     April 2012
 158 *
 159 *     .. Scalar Arguments ..
 160       INTEGER            INFO, LDA, LWORK, M, N
 161 *     ..
 162 *     .. Array Arguments ..
 163       REAL               A( LDA, * ), TAU( * ), WORK( * )
 164 *     ..
 165 *
 166 *  =====================================================================
 167 *
 168 *     .. Parameters ..
 169       REAL               ZERO
 170       PARAMETER          ( ZERO = 0.0E+0 )
 171 *     ..
 172 *     .. Local Scalars ..
 173       LOGICAL            LQUERY
 174       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKMIN, LWKOPT,
 175      $                   M1, MU, NB, NBMIN, NX
 176 *     ..
 177 *     .. External Subroutines ..
 178       EXTERNAL           XERBLA, SLARZB, SLARZT, SLATRZ
 179 *     ..
 180 *     .. Intrinsic Functions ..
 181       INTRINSIC          MAX, MIN
 182 *     ..
 183 *     .. External Functions ..
 184       INTEGER            ILAENV
 185       EXTERNAL           ILAENV
 186 *     ..
 187 *     .. Executable Statements ..
 188 *
 189 *     Test the input arguments
 190 *
 191       INFO = 0
 192       LQUERY = ( LWORK.EQ.-1 )
 193       IF( M.LT.0 ) THEN
 194          INFO = -1
 195       ELSE IF( N.LT.M ) THEN
 196          INFO = -2
 197       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 198          INFO = -4
 199       END IF
 200 *
 201       IF( INFO.EQ.0 ) THEN
 202          IF( M.EQ.0 .OR. M.EQ.N ) THEN
 203             LWKOPT = 1
 204             LWKMIN = 1
 205          ELSE
 206 *
 207 *           Determine the block size.
 208 *
 209             NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
 210             LWKOPT = M*NB
 211             LWKMIN = MAX( 1, M )
 212          END IF
 213          WORK( 1 ) = LWKOPT
 214 *
 215          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
 216             INFO = -7
 217          END IF
 218       END IF
 219 *
 220       IF( INFO.NE.0 ) THEN
 221          CALL XERBLA( 'STZRZF', -INFO )
 222          RETURN
 223       ELSE IF( LQUERY ) THEN
 224          RETURN
 225       END IF
 226 *
 227 *     Quick return if possible
 228 *
 229       IF( M.EQ.0 ) THEN
 230          RETURN
 231       ELSE IF( M.EQ.N ) THEN
 232          DO 10 I = 1, N
 233             TAU( I ) = ZERO
 234    10    CONTINUE
 235          RETURN
 236       END IF
 237 *
 238       NBMIN = 2
 239       NX = 1
 240       IWS = M
 241       IF( NB.GT.1 .AND. NB.LT.M ) THEN
 242 *
 243 *        Determine when to cross over from blocked to unblocked code.
 244 *
 245          NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) )
 246          IF( NX.LT.M ) THEN
 247 *
 248 *           Determine if workspace is large enough for blocked code.
 249 *
 250             LDWORK = M
 251             IWS = LDWORK*NB
 252             IF( LWORK.LT.IWS ) THEN
 253 *
 254 *              Not enough workspace to use optimal NB:  reduce NB and
 255 *              determine the minimum value of NB.
 256 *
 257                NB = LWORK / LDWORK
 258                NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1,
 259      $                 -1 ) )
 260             END IF
 261          END IF
 262       END IF
 263 *
 264       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
 265 *
 266 *        Use blocked code initially.
 267 *        The last kk rows are handled by the block method.
 268 *
 269          M1 = MIN( M+1, N )
 270          KI = ( ( M-NX-1 ) / NB )*NB
 271          KK = MIN( M, KI+NB )
 272 *
 273          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
 274             IB = MIN( M-I+1, NB )
 275 *
 276 *           Compute the TZ factorization of the current block
 277 *           A(i:i+ib-1,i:n)
 278 *
 279             CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
 280      $                   WORK )
 281             IF( I.GT.1 ) THEN
 282 *
 283 *              Form the triangular factor of the block reflector
 284 *              H = H(i+ib-1) . . . H(i+1) H(i)
 285 *
 286                CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
 287      $                      LDA, TAU( I ), WORK, LDWORK )
 288 *
 289 *              Apply H to A(1:i-1,i:n) from the right
 290 *
 291                CALL SLARZB( 'Right', 'No transpose', 'Backward',
 292      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
 293      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
 294      $                      WORK( IB+1 ), LDWORK )
 295             END IF
 296    20    CONTINUE
 297          MU = I + NB - 1
 298       ELSE
 299          MU = M
 300       END IF
 301 *
 302 *     Use unblocked code to factor the last or only block
 303 *
 304       IF( MU.GT.0 )
 305      $   CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
 306 *
 307       WORK( 1 ) = LWKOPT
 308 *
 309       RETURN
 310 *
 311 *     End of STZRZF
 312 *
 313       END