001:       SUBROUTINE SORGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, K, LDA, LWORK, M, N
009: *     ..
010: *     .. Array Arguments ..
011:       REAL               A( LDA, * ), TAU( * ), WORK( * )
012: *     ..
013: *
014: *  Purpose
015: *  =======
016: *
017: *  SORGLQ generates an M-by-N real matrix Q with orthonormal rows,
018: *  which is defined as the first M rows of a product of K elementary
019: *  reflectors of order N
020: *
021: *        Q  =  H(k) . . . H(2) H(1)
022: *
023: *  as returned by SGELQF.
024: *
025: *  Arguments
026: *  =========
027: *
028: *  M       (input) INTEGER
029: *          The number of rows of the matrix Q. M >= 0.
030: *
031: *  N       (input) INTEGER
032: *          The number of columns of the matrix Q. N >= M.
033: *
034: *  K       (input) INTEGER
035: *          The number of elementary reflectors whose product defines the
036: *          matrix Q. M >= K >= 0.
037: *
038: *  A       (input/output) REAL array, dimension (LDA,N)
039: *          On entry, the i-th row must contain the vector which defines
040: *          the elementary reflector H(i), for i = 1,2,...,k, as returned
041: *          by SGELQF in the first k rows of its array argument A.
042: *          On exit, the M-by-N matrix Q.
043: *
044: *  LDA     (input) INTEGER
045: *          The first dimension of the array A. LDA >= max(1,M).
046: *
047: *  TAU     (input) REAL array, dimension (K)
048: *          TAU(i) must contain the scalar factor of the elementary
049: *          reflector H(i), as returned by SGELQF.
050: *
051: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
052: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
053: *
054: *  LWORK   (input) INTEGER
055: *          The dimension of the array WORK. LWORK >= max(1,M).
056: *          For optimum performance LWORK >= M*NB, where NB is
057: *          the optimal blocksize.
058: *
059: *          If LWORK = -1, then a workspace query is assumed; the routine
060: *          only calculates the optimal size of the WORK array, returns
061: *          this value as the first entry of the WORK array, and no error
062: *          message related to LWORK is issued by XERBLA.
063: *
064: *  INFO    (output) INTEGER
065: *          = 0:  successful exit
066: *          < 0:  if INFO = -i, the i-th argument has an illegal value
067: *
068: *  =====================================================================
069: *
070: *     .. Parameters ..
071:       REAL               ZERO
072:       PARAMETER          ( ZERO = 0.0E+0 )
073: *     ..
074: *     .. Local Scalars ..
075:       LOGICAL            LQUERY
076:       INTEGER            I, IB, IINFO, IWS, J, KI, KK, L, LDWORK,
077:      $                   LWKOPT, NB, NBMIN, NX
078: *     ..
079: *     .. External Subroutines ..
080:       EXTERNAL           SLARFB, SLARFT, SORGL2, XERBLA
081: *     ..
082: *     .. Intrinsic Functions ..
083:       INTRINSIC          MAX, MIN
084: *     ..
085: *     .. External Functions ..
086:       INTEGER            ILAENV
087:       EXTERNAL           ILAENV
088: *     ..
089: *     .. Executable Statements ..
090: *
091: *     Test the input arguments
092: *
093:       INFO = 0
094:       NB = ILAENV( 1, 'SORGLQ', ' ', M, N, K, -1 )
095:       LWKOPT = MAX( 1, M )*NB
096:       WORK( 1 ) = LWKOPT
097:       LQUERY = ( LWORK.EQ.-1 )
098:       IF( M.LT.0 ) THEN
099:          INFO = -1
100:       ELSE IF( N.LT.M ) THEN
101:          INFO = -2
102:       ELSE IF( K.LT.0 .OR. K.GT.M ) THEN
103:          INFO = -3
104:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
105:          INFO = -5
106:       ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
107:          INFO = -8
108:       END IF
109:       IF( INFO.NE.0 ) THEN
110:          CALL XERBLA( 'SORGLQ', -INFO )
111:          RETURN
112:       ELSE IF( LQUERY ) THEN
113:          RETURN
114:       END IF
115: *
116: *     Quick return if possible
117: *
118:       IF( M.LE.0 ) THEN
119:          WORK( 1 ) = 1
120:          RETURN
121:       END IF
122: *
123:       NBMIN = 2
124:       NX = 0
125:       IWS = M
126:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
127: *
128: *        Determine when to cross over from blocked to unblocked code.
129: *
130:          NX = MAX( 0, ILAENV( 3, 'SORGLQ', ' ', M, N, K, -1 ) )
131:          IF( NX.LT.K ) THEN
132: *
133: *           Determine if workspace is large enough for blocked code.
134: *
135:             LDWORK = M
136:             IWS = LDWORK*NB
137:             IF( LWORK.LT.IWS ) THEN
138: *
139: *              Not enough workspace to use optimal NB:  reduce NB and
140: *              determine the minimum value of NB.
141: *
142:                NB = LWORK / LDWORK
143:                NBMIN = MAX( 2, ILAENV( 2, 'SORGLQ', ' ', M, N, K, -1 ) )
144:             END IF
145:          END IF
146:       END IF
147: *
148:       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
149: *
150: *        Use blocked code after the last block.
151: *        The first kk rows are handled by the block method.
152: *
153:          KI = ( ( K-NX-1 ) / NB )*NB
154:          KK = MIN( K, KI+NB )
155: *
156: *        Set A(kk+1:m,1:kk) to zero.
157: *
158:          DO 20 J = 1, KK
159:             DO 10 I = KK + 1, M
160:                A( I, J ) = ZERO
161:    10       CONTINUE
162:    20    CONTINUE
163:       ELSE
164:          KK = 0
165:       END IF
166: *
167: *     Use unblocked code for the last or only block.
168: *
169:       IF( KK.LT.M )
170:      $   CALL SORGL2( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA,
171:      $                TAU( KK+1 ), WORK, IINFO )
172: *
173:       IF( KK.GT.0 ) THEN
174: *
175: *        Use blocked code
176: *
177:          DO 50 I = KI + 1, 1, -NB
178:             IB = MIN( NB, K-I+1 )
179:             IF( I+IB.LE.M ) THEN
180: *
181: *              Form the triangular factor of the block reflector
182: *              H = H(i) H(i+1) . . . H(i+ib-1)
183: *
184:                CALL SLARFT( 'Forward', 'Rowwise', N-I+1, IB, A( I, I ),
185:      $                      LDA, TAU( I ), WORK, LDWORK )
186: *
187: *              Apply H' to A(i+ib:m,i:n) from the right
188: *
189:                CALL SLARFB( 'Right', 'Transpose', 'Forward', 'Rowwise',
190:      $                      M-I-IB+1, N-I+1, IB, A( I, I ), LDA, WORK,
191:      $                      LDWORK, A( I+IB, I ), LDA, WORK( IB+1 ),
192:      $                      LDWORK )
193:             END IF
194: *
195: *           Apply H' to columns i:n of current block
196: *
197:             CALL SORGL2( IB, N-I+1, IB, A( I, I ), LDA, TAU( I ), WORK,
198:      $                   IINFO )
199: *
200: *           Set columns 1:i-1 of current block to zero
201: *
202:             DO 40 J = 1, I - 1
203:                DO 30 L = I, I + IB - 1
204:                   A( L, J ) = ZERO
205:    30          CONTINUE
206:    40       CONTINUE
207:    50    CONTINUE
208:       END IF
209: *
210:       WORK( 1 ) = IWS
211:       RETURN
212: *
213: *     End of SORGLQ
214: *
215:       END
216: