SRC/sgeqlf.f

   1 *> \brief \b SGEQLF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGEQLF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqlf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqlf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqlf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGEQLF( 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 *> SGEQLF computes a QL factorization of a real M-by-N matrix A:
  37 *> A = Q * L.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] M
  44 *> \verbatim
  45 *>          M is INTEGER
  46 *>          The number of rows of the matrix A.  M >= 0.
  47 *> \endverbatim
  48 *>
  49 *> \param[in] N
  50 *> \verbatim
  51 *>          N is INTEGER
  52 *>          The number of columns of the matrix A.  N >= 0.
  53 *> \endverbatim
  54 *>
  55 *> \param[in,out] A
  56 *> \verbatim
  57 *>          A is REAL array, dimension (LDA,N)
  58 *>          On entry, the M-by-N matrix A.
  59 *>          On exit,
  60 *>          if m >= n, the lower triangle of the subarray
  61 *>          A(m-n+1:m,1:n) contains the N-by-N lower triangular matrix L;
  62 *>          if m <= n, the elements on and below the (n-m)-th
  63 *>          superdiagonal contain the M-by-N lower trapezoidal matrix L;
  64 *>          the remaining elements, with the array TAU, represent the
  65 *>          orthogonal matrix Q as a product of elementary reflectors
  66 *>          (see Further Details).
  67 *> \endverbatim
  68 *>
  69 *> \param[in] LDA
  70 *> \verbatim
  71 *>          LDA is INTEGER
  72 *>          The leading dimension of the array A.  LDA >= max(1,M).
  73 *> \endverbatim
  74 *>
  75 *> \param[out] TAU
  76 *> \verbatim
  77 *>          TAU is REAL array, dimension (min(M,N))
  78 *>          The scalar factors of the elementary reflectors (see Further
  79 *>          Details).
  80 *> \endverbatim
  81 *>
  82 *> \param[out] WORK
  83 *> \verbatim
  84 *>          WORK is REAL array, dimension (MAX(1,LWORK))
  85 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] LWORK
  89 *> \verbatim
  90 *>          LWORK is INTEGER
  91 *>          The dimension of the array WORK.  LWORK >= max(1,N).
  92 *>          For optimum performance LWORK >= N*NB, where NB is the
  93 *>          optimal blocksize.
  94 *>
  95 *>          If LWORK = -1, then a workspace query is assumed; the routine
  96 *>          only calculates the optimal size of the WORK array, returns
  97 *>          this value as the first entry of the WORK array, and no error
  98 *>          message related to LWORK is issued by XERBLA.
  99 *> \endverbatim
 100 *>
 101 *> \param[out] INFO
 102 *> \verbatim
 103 *>          INFO is INTEGER
 104 *>          = 0:  successful exit
 105 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 106 *> \endverbatim
 107 *
 108 *  Authors:
 109 *  ========
 110 *
 111 *> \author Univ. of Tennessee
 112 *> \author Univ. of California Berkeley
 113 *> \author Univ. of Colorado Denver
 114 *> \author NAG Ltd.
 115 *
 116 *> \date November 2011
 117 *
 118 *> \ingroup realGEcomputational
 119 *
 120 *> \par Further Details:
 121 *  =====================
 122 *>
 123 *> \verbatim
 124 *>
 125 *>  The matrix Q is represented as a product of elementary reflectors
 126 *>
 127 *>     Q = H(k) . . . H(2) H(1), where k = min(m,n).
 128 *>
 129 *>  Each H(i) has the form
 130 *>
 131 *>     H(i) = I - tau * v * v**T
 132 *>
 133 *>  where tau is a real scalar, and v is a real vector with
 134 *>  v(m-k+i+1:m) = 0 and v(m-k+i) = 1; v(1:m-k+i-1) is stored on exit in
 135 *>  A(1:m-k+i-1,n-k+i), and tau in TAU(i).
 136 *> \endverbatim
 137 *>
 138 *  =====================================================================
 139       SUBROUTINE SGEQLF( M, N, A, LDA, TAU, WORK, LWORK, INFO )
 140 *
 141 *  -- LAPACK computational routine (version 3.4.0) --
 142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 144 *     November 2011
 145 *
 146 *     .. Scalar Arguments ..
 147       INTEGER            INFO, LDA, LWORK, M, N
 148 *     ..
 149 *     .. Array Arguments ..
 150       REAL               A( LDA, * ), TAU( * ), WORK( * )
 151 *     ..
 152 *
 153 *  =====================================================================
 154 *
 155 *     .. Local Scalars ..
 156       LOGICAL            LQUERY
 157       INTEGER            I, IB, IINFO, IWS, K, KI, KK, LDWORK, LWKOPT,
 158      $                   MU, NB, NBMIN, NU, NX
 159 *     ..
 160 *     .. External Subroutines ..
 161       EXTERNAL           SGEQL2, SLARFB, SLARFT, XERBLA
 162 *     ..
 163 *     .. Intrinsic Functions ..
 164       INTRINSIC          MAX, MIN
 165 *     ..
 166 *     .. External Functions ..
 167       INTEGER            ILAENV
 168       EXTERNAL           ILAENV
 169 *     ..
 170 *     .. Executable Statements ..
 171 *
 172 *     Test the input arguments
 173 *
 174       INFO = 0
 175       LQUERY = ( LWORK.EQ.-1 )
 176       IF( M.LT.0 ) THEN
 177          INFO = -1
 178       ELSE IF( N.LT.0 ) THEN
 179          INFO = -2
 180       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 181          INFO = -4
 182       END IF
 183 *
 184       IF( INFO.EQ.0 ) THEN
 185          K = MIN( M, N )
 186          IF( K.EQ.0 ) THEN
 187             LWKOPT = 1
 188          ELSE
 189             NB = ILAENV( 1, 'SGEQLF', ' ', M, N, -1, -1 )
 190             LWKOPT = N*NB
 191          END IF
 192          WORK( 1 ) = LWKOPT
 193 *
 194          IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
 195             INFO = -7
 196          END IF
 197       END IF
 198 *
 199       IF( INFO.NE.0 ) THEN
 200          CALL XERBLA( 'SGEQLF', -INFO )
 201          RETURN
 202       ELSE IF( LQUERY ) THEN
 203          RETURN
 204       END IF
 205 *
 206 *     Quick return if possible
 207 *
 208       IF( K.EQ.0 ) THEN
 209          RETURN
 210       END IF
 211 *
 212       NBMIN = 2
 213       NX = 1
 214       IWS = N
 215       IF( NB.GT.1 .AND. NB.LT.K ) THEN
 216 *
 217 *        Determine when to cross over from blocked to unblocked code.
 218 *
 219          NX = MAX( 0, ILAENV( 3, 'SGEQLF', ' ', M, N, -1, -1 ) )
 220          IF( NX.LT.K ) THEN
 221 *
 222 *           Determine if workspace is large enough for blocked code.
 223 *
 224             LDWORK = N
 225             IWS = LDWORK*NB
 226             IF( LWORK.LT.IWS ) THEN
 227 *
 228 *              Not enough workspace to use optimal NB:  reduce NB and
 229 *              determine the minimum value of NB.
 230 *
 231                NB = LWORK / LDWORK
 232                NBMIN = MAX( 2, ILAENV( 2, 'SGEQLF', ' ', M, N, -1,
 233      $                 -1 ) )
 234             END IF
 235          END IF
 236       END IF
 237 *
 238       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
 239 *
 240 *        Use blocked code initially.
 241 *        The last kk columns are handled by the block method.
 242 *
 243          KI = ( ( K-NX-1 ) / NB )*NB
 244          KK = MIN( K, KI+NB )
 245 *
 246          DO 10 I = K - KK + KI + 1, K - KK + 1, -NB
 247             IB = MIN( K-I+1, NB )
 248 *
 249 *           Compute the QL factorization of the current block
 250 *           A(1:m-k+i+ib-1,n-k+i:n-k+i+ib-1)
 251 *
 252             CALL SGEQL2( M-K+I+IB-1, IB, A( 1, N-K+I ), LDA, TAU( I ),
 253      $                   WORK, IINFO )
 254             IF( N-K+I.GT.1 ) THEN
 255 *
 256 *              Form the triangular factor of the block reflector
 257 *              H = H(i+ib-1) . . . H(i+1) H(i)
 258 *
 259                CALL SLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB,
 260      $                      A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK )
 261 *
 262 *              Apply H**T to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
 263 *
 264                CALL SLARFB( 'Left', 'Transpose', 'Backward',
 265      $                      'Columnwise', M-K+I+IB-1, N-K+I-1, IB,
 266      $                      A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA,
 267      $                      WORK( IB+1 ), LDWORK )
 268             END IF
 269    10    CONTINUE
 270          MU = M - K + I + NB - 1
 271          NU = N - K + I + NB - 1
 272       ELSE
 273          MU = M
 274          NU = N
 275       END IF
 276 *
 277 *     Use unblocked code to factor the last or only block
 278 *
 279       IF( MU.GT.0 .AND. NU.GT.0 )
 280      $   CALL SGEQL2( MU, NU, A, LDA, TAU, WORK, IINFO )
 281 *
 282       WORK( 1 ) = IWS
 283       RETURN
 284 *
 285 *     End of SGEQLF
 286 *
 287       END