001:       SUBROUTINE SORMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
002:      $                   LDC, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          SIDE, TRANS, VECT
010:       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               A( LDA, * ), C( LDC, * ), TAU( * ),
014:      $                   WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  If VECT = 'Q', SORMBR overwrites the general real M-by-N matrix C
021: *  with
022: *                  SIDE = 'L'     SIDE = 'R'
023: *  TRANS = 'N':      Q * C          C * Q
024: *  TRANS = 'T':      Q**T * C       C * Q**T
025: *
026: *  If VECT = 'P', SORMBR overwrites the general real M-by-N matrix C
027: *  with
028: *                  SIDE = 'L'     SIDE = 'R'
029: *  TRANS = 'N':      P * C          C * P
030: *  TRANS = 'T':      P**T * C       C * P**T
031: *
032: *  Here Q and P**T are the orthogonal matrices determined by SGEBRD when
033: *  reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and
034: *  P**T are defined as products of elementary reflectors H(i) and G(i)
035: *  respectively.
036: *
037: *  Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
038: *  order of the orthogonal matrix Q or P**T that is applied.
039: *
040: *  If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
041: *  if nq >= k, Q = H(1) H(2) . . . H(k);
042: *  if nq < k, Q = H(1) H(2) . . . H(nq-1).
043: *
044: *  If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
045: *  if k < nq, P = G(1) G(2) . . . G(k);
046: *  if k >= nq, P = G(1) G(2) . . . G(nq-1).
047: *
048: *  Arguments
049: *  =========
050: *
051: *  VECT    (input) CHARACTER*1
052: *          = 'Q': apply Q or Q**T;
053: *          = 'P': apply P or P**T.
054: *
055: *  SIDE    (input) CHARACTER*1
056: *          = 'L': apply Q, Q**T, P or P**T from the Left;
057: *          = 'R': apply Q, Q**T, P or P**T from the Right.
058: *
059: *  TRANS   (input) CHARACTER*1
060: *          = 'N':  No transpose, apply Q  or P;
061: *          = 'T':  Transpose, apply Q**T or P**T.
062: *
063: *  M       (input) INTEGER
064: *          The number of rows of the matrix C. M >= 0.
065: *
066: *  N       (input) INTEGER
067: *          The number of columns of the matrix C. N >= 0.
068: *
069: *  K       (input) INTEGER
070: *          If VECT = 'Q', the number of columns in the original
071: *          matrix reduced by SGEBRD.
072: *          If VECT = 'P', the number of rows in the original
073: *          matrix reduced by SGEBRD.
074: *          K >= 0.
075: *
076: *  A       (input) REAL array, dimension
077: *                                (LDA,min(nq,K)) if VECT = 'Q'
078: *                                (LDA,nq)        if VECT = 'P'
079: *          The vectors which define the elementary reflectors H(i) and
080: *          G(i), whose products determine the matrices Q and P, as
081: *          returned by SGEBRD.
082: *
083: *  LDA     (input) INTEGER
084: *          The leading dimension of the array A.
085: *          If VECT = 'Q', LDA >= max(1,nq);
086: *          if VECT = 'P', LDA >= max(1,min(nq,K)).
087: *
088: *  TAU     (input) REAL array, dimension (min(nq,K))
089: *          TAU(i) must contain the scalar factor of the elementary
090: *          reflector H(i) or G(i) which determines Q or P, as returned
091: *          by SGEBRD in the array argument TAUQ or TAUP.
092: *
093: *  C       (input/output) REAL array, dimension (LDC,N)
094: *          On entry, the M-by-N matrix C.
095: *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q
096: *          or P*C or P**T*C or C*P or C*P**T.
097: *
098: *  LDC     (input) INTEGER
099: *          The leading dimension of the array C. LDC >= max(1,M).
100: *
101: *  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
102: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103: *
104: *  LWORK   (input) INTEGER
105: *          The dimension of the array WORK.
106: *          If SIDE = 'L', LWORK >= max(1,N);
107: *          if SIDE = 'R', LWORK >= max(1,M).
108: *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
109: *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
110: *          blocksize.
111: *
112: *          If LWORK = -1, then a workspace query is assumed; the routine
113: *          only calculates the optimal size of the WORK array, returns
114: *          this value as the first entry of the WORK array, and no error
115: *          message related to LWORK is issued by XERBLA.
116: *
117: *  INFO    (output) INTEGER
118: *          = 0:  successful exit
119: *          < 0:  if INFO = -i, the i-th argument had an illegal value
120: *
121: *  =====================================================================
122: *
123: *     .. Local Scalars ..
124:       LOGICAL            APPLYQ, LEFT, LQUERY, NOTRAN
125:       CHARACTER          TRANST
126:       INTEGER            I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW
127: *     ..
128: *     .. External Functions ..
129:       LOGICAL            LSAME
130:       INTEGER            ILAENV
131:       EXTERNAL           ILAENV, LSAME
132: *     ..
133: *     .. External Subroutines ..
134:       EXTERNAL           SORMLQ, SORMQR, XERBLA
135: *     ..
136: *     .. Intrinsic Functions ..
137:       INTRINSIC          MAX, MIN
138: *     ..
139: *     .. Executable Statements ..
140: *
141: *     Test the input arguments
142: *
143:       INFO = 0
144:       APPLYQ = LSAME( VECT, 'Q' )
145:       LEFT = LSAME( SIDE, 'L' )
146:       NOTRAN = LSAME( TRANS, 'N' )
147:       LQUERY = ( LWORK.EQ.-1 )
148: *
149: *     NQ is the order of Q or P and NW is the minimum dimension of WORK
150: *
151:       IF( LEFT ) THEN
152:          NQ = M
153:          NW = N
154:       ELSE
155:          NQ = N
156:          NW = M
157:       END IF
158:       IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
159:          INFO = -1
160:       ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
161:          INFO = -2
162:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
163:          INFO = -3
164:       ELSE IF( M.LT.0 ) THEN
165:          INFO = -4
166:       ELSE IF( N.LT.0 ) THEN
167:          INFO = -5
168:       ELSE IF( K.LT.0 ) THEN
169:          INFO = -6
170:       ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR.
171:      $         ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) )
172:      $          THEN
173:          INFO = -8
174:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
175:          INFO = -11
176:       ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
177:          INFO = -13
178:       END IF
179: *
180:       IF( INFO.EQ.0 ) THEN
181:          IF( APPLYQ ) THEN
182:             IF( LEFT ) THEN
183:                NB = ILAENV( 1, 'SORMQR', SIDE // TRANS, M-1, N, M-1,
184:      $                      -1 )
185:             ELSE
186:                NB = ILAENV( 1, 'SORMQR', SIDE // TRANS, M, N-1, N-1,
187:      $                      -1 )
188:             END IF   
189:          ELSE
190:             IF( LEFT ) THEN
191:                NB = ILAENV( 1, 'SORMLQ', SIDE // TRANS, M-1, N, M-1,
192:      $                      -1 ) 
193:             ELSE
194:                NB = ILAENV( 1, 'SORMLQ', SIDE // TRANS, M, N-1, N-1,
195:      $                      -1 )
196:             END IF
197:          END IF
198:          LWKOPT = MAX( 1, NW )*NB
199:          WORK( 1 ) = LWKOPT 
200:       END IF
201: *
202:       IF( INFO.NE.0 ) THEN
203:          CALL XERBLA( 'SORMBR', -INFO )
204:          RETURN
205:       ELSE IF( LQUERY ) THEN
206:          RETURN
207:       END IF
208: *
209: *     Quick return if possible
210: *
211:       WORK( 1 ) = 1
212:       IF( M.EQ.0 .OR. N.EQ.0 )
213:      $   RETURN
214: *
215:       IF( APPLYQ ) THEN
216: *
217: *        Apply Q
218: *
219:          IF( NQ.GE.K ) THEN
220: *
221: *           Q was determined by a call to SGEBRD with nq >= k
222: *
223:             CALL SORMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
224:      $                   WORK, LWORK, IINFO )
225:          ELSE IF( NQ.GT.1 ) THEN
226: *
227: *           Q was determined by a call to SGEBRD with nq < k
228: *
229:             IF( LEFT ) THEN
230:                MI = M - 1
231:                NI = N
232:                I1 = 2
233:                I2 = 1
234:             ELSE
235:                MI = M
236:                NI = N - 1
237:                I1 = 1
238:                I2 = 2
239:             END IF
240:             CALL SORMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU,
241:      $                   C( I1, I2 ), LDC, WORK, LWORK, IINFO )
242:          END IF
243:       ELSE
244: *
245: *        Apply P
246: *
247:          IF( NOTRAN ) THEN
248:             TRANST = 'T'
249:          ELSE
250:             TRANST = 'N'
251:          END IF
252:          IF( NQ.GT.K ) THEN
253: *
254: *           P was determined by a call to SGEBRD with nq > k
255: *
256:             CALL SORMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC,
257:      $                   WORK, LWORK, IINFO )
258:          ELSE IF( NQ.GT.1 ) THEN
259: *
260: *           P was determined by a call to SGEBRD with nq <= k
261: *
262:             IF( LEFT ) THEN
263:                MI = M - 1
264:                NI = N
265:                I1 = 2
266:                I2 = 1
267:             ELSE
268:                MI = M
269:                NI = N - 1
270:                I1 = 1
271:                I2 = 2
272:             END IF
273:             CALL SORMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA,
274:      $                   TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO )
275:          END IF
276:       END IF
277:       WORK( 1 ) = LWKOPT
278:       RETURN
279: *
280: *     End of SORMBR
281: *
282:       END
283: