SRC/dorgql.f

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