001:       SUBROUTINE STZRZF( M, N, 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, LDA, LWORK, M, N
009: *     ..
010: *     .. Array Arguments ..
011:       REAL               A( LDA, * ), TAU( * ), WORK( * )
012: *     ..
013: *
014: *  Purpose
015: *  =======
016: *
017: *  STZRZF reduces the M-by-N ( M<=N ) real upper trapezoidal matrix A
018: *  to upper triangular form by means of orthogonal transformations.
019: *
020: *  The upper trapezoidal matrix A is factored as
021: *
022: *     A = ( R  0 ) * Z,
023: *
024: *  where Z is an N-by-N orthogonal matrix and R is an M-by-M upper
025: *  triangular matrix.
026: *
027: *  Arguments
028: *  =========
029: *
030: *  M       (input) INTEGER
031: *          The number of rows of the matrix A.  M >= 0.
032: *
033: *  N       (input) INTEGER
034: *          The number of columns of the matrix A.  N >= M.
035: *
036: *  A       (input/output) REAL array, dimension (LDA,N)
037: *          On entry, the leading M-by-N upper trapezoidal part of the
038: *          array A must contain the matrix to be factorized.
039: *          On exit, the leading M-by-M upper triangular part of A
040: *          contains the upper triangular matrix R, and elements M+1 to
041: *          N of the first M rows of A, with the array TAU, represent the
042: *          orthogonal matrix Z as a product of M elementary reflectors.
043: *
044: *  LDA     (input) INTEGER
045: *          The leading dimension of the array A.  LDA >= max(1,M).
046: *
047: *  TAU     (output) REAL array, dimension (M)
048: *          The scalar factors of the elementary reflectors.
049: *
050: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
051: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
052: *
053: *  LWORK   (input) INTEGER
054: *          The dimension of the array WORK.  LWORK >= max(1,M).
055: *          For optimum performance LWORK >= M*NB, where NB is
056: *          the optimal blocksize.
057: *
058: *          If LWORK = -1, then a workspace query is assumed; the routine
059: *          only calculates the optimal size of the WORK array, returns
060: *          this value as the first entry of the WORK array, and no error
061: *          message related to LWORK is issued by XERBLA.
062: *
063: *  INFO    (output) INTEGER
064: *          = 0:  successful exit
065: *          < 0:  if INFO = -i, the i-th argument had an illegal value
066: *
067: *  Further Details
068: *  ===============
069: *
070: *  Based on contributions by
071: *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
072: *
073: *  The factorization is obtained by Householder's method.  The kth
074: *  transformation matrix, Z( k ), which is used to introduce zeros into
075: *  the ( m - k + 1 )th row of A, is given in the form
076: *
077: *     Z( k ) = ( I     0   ),
078: *              ( 0  T( k ) )
079: *
080: *  where
081: *
082: *     T( k ) = I - tau*u( k )*u( k )',   u( k ) = (   1    ),
083: *                                                 (   0    )
084: *                                                 ( z( k ) )
085: *
086: *  tau is a scalar and z( k ) is an ( n - m ) element vector.
087: *  tau and z( k ) are chosen to annihilate the elements of the kth row
088: *  of X.
089: *
090: *  The scalar tau is returned in the kth element of TAU and the vector
091: *  u( k ) in the kth row of A, such that the elements of z( k ) are
092: *  in  a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in
093: *  the upper triangular part of A.
094: *
095: *  Z is given by
096: *
097: *     Z =  Z( 1 ) * Z( 2 ) * ... * Z( m ).
098: *
099: *  =====================================================================
100: *
101: *     .. Parameters ..
102:       REAL               ZERO
103:       PARAMETER          ( ZERO = 0.0E+0 )
104: *     ..
105: *     .. Local Scalars ..
106:       LOGICAL            LQUERY
107:       INTEGER            I, IB, IWS, KI, KK, LDWORK, LWKOPT, M1, MU, NB,
108:      $                   NBMIN, NX
109: *     ..
110: *     .. External Subroutines ..
111:       EXTERNAL           SLARZB, SLARZT, SLATRZ, XERBLA
112: *     ..
113: *     .. Intrinsic Functions ..
114:       INTRINSIC          MAX, MIN
115: *     ..
116: *     .. External Functions ..
117:       INTEGER            ILAENV
118:       EXTERNAL           ILAENV
119: *     ..
120: *     .. Executable Statements ..
121: *
122: *     Test the input arguments
123: *
124:       INFO = 0
125:       LQUERY = ( LWORK.EQ.-1 )
126:       IF( M.LT.0 ) THEN
127:          INFO = -1
128:       ELSE IF( N.LT.M ) THEN
129:          INFO = -2
130:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
131:          INFO = -4
132:       END IF
133: *
134:       IF( INFO.EQ.0 ) THEN
135:          IF( M.EQ.0 .OR. M.EQ.N ) THEN
136:             LWKOPT = 1
137:          ELSE
138: *
139: *           Determine the block size.
140: *
141:             NB = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
142:             LWKOPT = M*NB
143:          END IF
144:          WORK( 1 ) = LWKOPT
145: *
146:          IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN
147:             INFO = -7
148:          END IF
149:       END IF
150: *
151:       IF( INFO.NE.0 ) THEN
152:          CALL XERBLA( 'STZRZF', -INFO )
153:          RETURN
154:       ELSE IF( LQUERY ) THEN
155:          RETURN
156:       END IF
157: *
158: *     Quick return if possible
159: *
160:       IF( M.EQ.0 ) THEN
161:          RETURN
162:       ELSE IF( M.EQ.N ) THEN
163:          DO 10 I = 1, N
164:             TAU( I ) = ZERO
165:    10    CONTINUE
166:          RETURN
167:       END IF
168: *
169:       NBMIN = 2
170:       NX = 1
171:       IWS = M
172:       IF( NB.GT.1 .AND. NB.LT.M ) THEN
173: *
174: *        Determine when to cross over from blocked to unblocked code.
175: *
176:          NX = MAX( 0, ILAENV( 3, 'SGERQF', ' ', M, N, -1, -1 ) )
177:          IF( NX.LT.M ) THEN
178: *
179: *           Determine if workspace is large enough for blocked code.
180: *
181:             LDWORK = M
182:             IWS = LDWORK*NB
183:             IF( LWORK.LT.IWS ) THEN
184: *
185: *              Not enough workspace to use optimal NB:  reduce NB and
186: *              determine the minimum value of NB.
187: *
188:                NB = LWORK / LDWORK
189:                NBMIN = MAX( 2, ILAENV( 2, 'SGERQF', ' ', M, N, -1,
190:      $                 -1 ) )
191:             END IF
192:          END IF
193:       END IF
194: *
195:       IF( NB.GE.NBMIN .AND. NB.LT.M .AND. NX.LT.M ) THEN
196: *
197: *        Use blocked code initially.
198: *        The last kk rows are handled by the block method.
199: *
200:          M1 = MIN( M+1, N )
201:          KI = ( ( M-NX-1 ) / NB )*NB
202:          KK = MIN( M, KI+NB )
203: *
204:          DO 20 I = M - KK + KI + 1, M - KK + 1, -NB
205:             IB = MIN( M-I+1, NB )
206: *
207: *           Compute the TZ factorization of the current block
208: *           A(i:i+ib-1,i:n)
209: *
210:             CALL SLATRZ( IB, N-I+1, N-M, A( I, I ), LDA, TAU( I ),
211:      $                   WORK )
212:             IF( I.GT.1 ) THEN
213: *
214: *              Form the triangular factor of the block reflector
215: *              H = H(i+ib-1) . . . H(i+1) H(i)
216: *
217:                CALL SLARZT( 'Backward', 'Rowwise', N-M, IB, A( I, M1 ),
218:      $                      LDA, TAU( I ), WORK, LDWORK )
219: *
220: *              Apply H to A(1:i-1,i:n) from the right
221: *
222:                CALL SLARZB( 'Right', 'No transpose', 'Backward',
223:      $                      'Rowwise', I-1, N-I+1, IB, N-M, A( I, M1 ),
224:      $                      LDA, WORK, LDWORK, A( 1, I ), LDA,
225:      $                      WORK( IB+1 ), LDWORK )
226:             END IF
227:    20    CONTINUE
228:          MU = I + NB - 1
229:       ELSE
230:          MU = M
231:       END IF
232: *
233: *     Use unblocked code to factor the last or only block
234: *
235:       IF( MU.GT.0 )
236:      $   CALL SLATRZ( MU, N, N-M, A, LDA, TAU, WORK )
237: *
238:       WORK( 1 ) = LWKOPT
239: *
240:       RETURN
241: *
242: *     End of STZRZF
243: *
244:       END
245: