2
SUBROUTINE DORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
5
* -- LAPACK routine (version 2.0) --
6
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
7
* Courant Institute, Argonne National Lab, and Rice University
10
* .. Scalar Arguments ..
12
INTEGER INFO, K, LDA, LDC, LWORK, M, N
14
* .. Array Arguments ..
15
DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ),
22
* DORMQL overwrites the general real M-by-N matrix C with
24
* SIDE = 'L' SIDE = 'R'
25
* TRANS = 'N': Q * C C * Q
26
* TRANS = 'T': Q**T * C C * Q**T
28
* where Q is a real orthogonal matrix defined as the product of k
29
* elementary reflectors
31
* Q = H(k) . . . H(2) H(1)
33
* as returned by DGEQLF. Q is of order M if SIDE = 'L' and of order N
39
* SIDE (input) CHARACTER*1
40
* = 'L': apply Q or Q**T from the Left;
41
* = 'R': apply Q or Q**T from the Right.
43
* TRANS (input) CHARACTER*1
44
* = 'N': No transpose, apply Q;
45
* = 'T': Transpose, apply Q**T.
48
* The number of rows of the matrix C. M >= 0.
51
* The number of columns of the matrix C. N >= 0.
54
* The number of elementary reflectors whose product defines
56
* If SIDE = 'L', M >= K >= 0;
57
* if SIDE = 'R', N >= K >= 0.
59
* A (input) DOUBLE PRECISION array, dimension (LDA,K)
60
* The i-th column must contain the vector which defines the
61
* elementary reflector H(i), for i = 1,2,...,k, as returned by
62
* DGEQLF in the last k columns of its array argument A.
63
* A is modified by the routine but restored on exit.
66
* The leading dimension of the array A.
67
* If SIDE = 'L', LDA >= max(1,M);
68
* if SIDE = 'R', LDA >= max(1,N).
70
* TAU (input) DOUBLE PRECISION array, dimension (K)
71
* TAU(i) must contain the scalar factor of the elementary
72
* reflector H(i), as returned by DGEQLF.
74
* C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
75
* On entry, the M-by-N matrix C.
76
* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
79
* The leading dimension of the array C. LDC >= max(1,M).
81
* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
82
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
84
* LWORK (input) INTEGER
85
* The dimension of the array WORK.
86
* If SIDE = 'L', LWORK >= max(1,N);
87
* if SIDE = 'R', LWORK >= max(1,M).
88
* For optimum performance LWORK >= N*NB if SIDE = 'L', and
89
* LWORK >= M*NB if SIDE = 'R', where NB is the optimal
92
* INFO (output) INTEGER
93
* = 0: successful exit
94
* < 0: if INFO = -i, the i-th argument had an illegal value
96
* =====================================================================
100
PARAMETER ( NBMAX = 64, LDT = NBMAX+1 )
102
* .. Local Scalars ..
104
INTEGER I, I1, I2, I3, IB, IINFO, IWS, LDWORK, MI, NB,
108
DOUBLE PRECISION T( LDT, NBMAX )
110
* .. External Functions ..
113
EXTERNAL LSAME, ILAENV
115
* .. External Subroutines ..
116
EXTERNAL DLARFB, DLARFT, DORM2L, XERBLA
118
* .. Intrinsic Functions ..
121
* .. Executable Statements ..
123
* Test the input arguments
126
LEFT = LSAME( SIDE, 'L' )
127
NOTRAN = LSAME( TRANS, 'N' )
129
* NQ is the order of Q and NW is the minimum dimension of WORK
138
IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
140
ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
142
ELSE IF( M.LT.0 ) THEN
144
ELSE IF( N.LT.0 ) THEN
146
ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
148
ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN
150
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
152
ELSE IF( LWORK.LT.MAX( 1, NW ) ) THEN
156
CALL XERBLA( 'DORMQL', -INFO )
160
* Quick return if possible
162
IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN
167
* Determine the block size. NB may be at most NBMAX, where NBMAX
168
* is used to define the local array T.
170
NB = MIN( NBMAX, ILAENV( 1, 'DORMQL', SIDE // TRANS, M, N, K,
174
IF( NB.GT.1 .AND. NB.LT.K ) THEN
176
IF( LWORK.LT.IWS ) THEN
178
NBMIN = MAX( 2, ILAENV( 2, 'DORMQL', SIDE // TRANS, M, N, K,
185
IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
189
CALL DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK,
195
IF( ( LEFT .AND. NOTRAN ) .OR.
196
$ ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) THEN
201
I1 = ( ( K-1 ) / NB )*NB + 1
213
IB = MIN( NB, K-I+1 )
215
* Form the triangular factor of the block reflector
216
* H = H(i+ib-1) . . . H(i+1) H(i)
218
CALL DLARFT( 'Backward', 'Columnwise', NQ-K+I+IB-1, IB,
219
$ A( 1, I ), LDA, TAU( I ), T, LDT )
222
* H or H' is applied to C(1:m-k+i+ib-1,1:n)
224
MI = M - K + I + IB - 1
227
* H or H' is applied to C(1:m,1:n-k+i+ib-1)
229
NI = N - K + I + IB - 1
234
CALL DLARFB( SIDE, TRANS, 'Backward', 'Columnwise', MI, NI,
235
$ IB, A( 1, I ), LDA, T, LDT, C, LDC, WORK,
246
SUBROUTINE DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
249
* -- LAPACK routine (version 2.0) --
250
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
251
* Courant Institute, Argonne National Lab, and Rice University
254
* .. Scalar Arguments ..
255
CHARACTER SIDE, TRANS
256
INTEGER INFO, K, LDA, LDC, M, N
258
* .. Array Arguments ..
259
DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
265
* DORM2L overwrites the general real m by n matrix C with
267
* Q * C if SIDE = 'L' and TRANS = 'N', or
269
* Q'* C if SIDE = 'L' and TRANS = 'T', or
271
* C * Q if SIDE = 'R' and TRANS = 'N', or
273
* C * Q' if SIDE = 'R' and TRANS = 'T',
275
* where Q is a real orthogonal matrix defined as the product of k
276
* elementary reflectors
278
* Q = H(k) . . . H(2) H(1)
280
* as returned by DGEQLF. Q is of order m if SIDE = 'L' and of order n
286
* SIDE (input) CHARACTER*1
287
* = 'L': apply Q or Q' from the Left
288
* = 'R': apply Q or Q' from the Right
290
* TRANS (input) CHARACTER*1
291
* = 'N': apply Q (No transpose)
292
* = 'T': apply Q' (Transpose)
295
* The number of rows of the matrix C. M >= 0.
298
* The number of columns of the matrix C. N >= 0.
301
* The number of elementary reflectors whose product defines
303
* If SIDE = 'L', M >= K >= 0;
304
* if SIDE = 'R', N >= K >= 0.
306
* A (input) DOUBLE PRECISION array, dimension (LDA,K)
307
* The i-th column must contain the vector which defines the
308
* elementary reflector H(i), for i = 1,2,...,k, as returned by
309
* DGEQLF in the last k columns of its array argument A.
310
* A is modified by the routine but restored on exit.
312
* LDA (input) INTEGER
313
* The leading dimension of the array A.
314
* If SIDE = 'L', LDA >= max(1,M);
315
* if SIDE = 'R', LDA >= max(1,N).
317
* TAU (input) DOUBLE PRECISION array, dimension (K)
318
* TAU(i) must contain the scalar factor of the elementary
319
* reflector H(i), as returned by DGEQLF.
321
* C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
322
* On entry, the m by n matrix C.
323
* On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q.
325
* LDC (input) INTEGER
326
* The leading dimension of the array C. LDC >= max(1,M).
328
* WORK (workspace) DOUBLE PRECISION array, dimension
332
* INFO (output) INTEGER
333
* = 0: successful exit
334
* < 0: if INFO = -i, the i-th argument had an illegal value
336
* =====================================================================
340
PARAMETER ( ONE = 1.0D+0 )
342
* .. Local Scalars ..
344
INTEGER I, I1, I2, I3, MI, NI, NQ
347
* .. External Functions ..
351
* .. External Subroutines ..
352
EXTERNAL DLARF, XERBLA
354
* .. Intrinsic Functions ..
357
* .. Executable Statements ..
359
* Test the input arguments
362
LEFT = LSAME( SIDE, 'L' )
363
NOTRAN = LSAME( TRANS, 'N' )
365
* NQ is the order of Q
372
IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
374
ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
376
ELSE IF( M.LT.0 ) THEN
378
ELSE IF( N.LT.0 ) THEN
380
ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
382
ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN
384
ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
388
CALL XERBLA( 'DORM2L', -INFO )
392
* Quick return if possible
394
IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 )
397
IF( ( LEFT .AND. NOTRAN ) .OR. ( .NOT.LEFT .AND. .NOT.NOTRAN ) )
417
* H(i) is applied to C(1:m-k+i,1:n)
422
* H(i) is applied to C(1:m,1:n-k+i)
431
CALL DLARF( SIDE, MI, NI, A( 1, I ), 1, TAU( I ), C, LDC,
added
removed
routines/lapack/dormql.f
