LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sormql.f
Go to the documentation of this file.
1*> \brief \b SORMQL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download SORMQL + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sormql.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sormql.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sormql.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
20* WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER SIDE, TRANS
24* INTEGER INFO, K, LDA, LDC, LWORK, M, N
25* ..
26* .. Array Arguments ..
27* REAL A( LDA, * ), C( LDC, * ), TAU( * ),
28* $ WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SORMQL overwrites the general real M-by-N matrix C with
38*>
39*> SIDE = 'L' SIDE = 'R'
40*> TRANS = 'N': Q * C C * Q
41*> TRANS = 'T': Q**T * C C * Q**T
42*>
43*> where Q is a real orthogonal matrix defined as the product of k
44*> elementary reflectors
45*>
46*> Q = H(k) . . . H(2) H(1)
47*>
48*> as returned by SGEQLF. Q is of order M if SIDE = 'L' and of order N
49*> if SIDE = 'R'.
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] SIDE
56*> \verbatim
57*> SIDE is CHARACTER*1
58*> = 'L': apply Q or Q**T from the Left;
59*> = 'R': apply Q or Q**T from the Right.
60*> \endverbatim
61*>
62*> \param[in] TRANS
63*> \verbatim
64*> TRANS is CHARACTER*1
65*> = 'N': No transpose, apply Q;
66*> = 'T': Transpose, apply Q**T.
67*> \endverbatim
68*>
69*> \param[in] M
70*> \verbatim
71*> M is INTEGER
72*> The number of rows of the matrix C. M >= 0.
73*> \endverbatim
74*>
75*> \param[in] N
76*> \verbatim
77*> N is INTEGER
78*> The number of columns of the matrix C. N >= 0.
79*> \endverbatim
80*>
81*> \param[in] K
82*> \verbatim
83*> K is INTEGER
84*> The number of elementary reflectors whose product defines
85*> the matrix Q.
86*> If SIDE = 'L', M >= K >= 0;
87*> if SIDE = 'R', N >= K >= 0.
88*> \endverbatim
89*>
90*> \param[in] A
91*> \verbatim
92*> A is REAL array, dimension (LDA,K)
93*> The i-th column must contain the vector which defines the
94*> elementary reflector H(i), for i = 1,2,...,k, as returned by
95*> SGEQLF in the last k columns of its array argument A.
96*> \endverbatim
97*>
98*> \param[in] LDA
99*> \verbatim
100*> LDA is INTEGER
101*> The leading dimension of the array A.
102*> If SIDE = 'L', LDA >= max(1,M);
103*> if SIDE = 'R', LDA >= max(1,N).
104*> \endverbatim
105*>
106*> \param[in] TAU
107*> \verbatim
108*> TAU is REAL array, dimension (K)
109*> TAU(i) must contain the scalar factor of the elementary
110*> reflector H(i), as returned by SGEQLF.
111*> \endverbatim
112*>
113*> \param[in,out] C
114*> \verbatim
115*> C is REAL array, dimension (LDC,N)
116*> On entry, the M-by-N matrix C.
117*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
118*> \endverbatim
119*>
120*> \param[in] LDC
121*> \verbatim
122*> LDC is INTEGER
123*> The leading dimension of the array C. LDC >= max(1,M).
124*> \endverbatim
125*>
126*> \param[out] WORK
127*> \verbatim
128*> WORK is REAL array, dimension (MAX(1,LWORK))
129*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
130*> \endverbatim
131*>
132*> \param[in] LWORK
133*> \verbatim
134*> LWORK is INTEGER
135*> The dimension of the array WORK.
136*> If SIDE = 'L', LWORK >= max(1,N);
137*> if SIDE = 'R', LWORK >= max(1,M).
138*> For good performance, LWORK should generally be larger.
139*>
140*> If LWORK = -1, then a workspace query is assumed; the routine
141*> only calculates the optimal size of the WORK array, returns
142*> this value as the first entry of the WORK array, and no error
143*> message related to LWORK is issued by XERBLA.
144*> \endverbatim
145*>
146*> \param[out] INFO
147*> \verbatim
148*> INFO is INTEGER
149*> = 0: successful exit
150*> < 0: if INFO = -i, the i-th argument had an illegal value
151*> \endverbatim
152*
153* Authors:
154* ========
155*
156*> \author Univ. of Tennessee
157*> \author Univ. of California Berkeley
158*> \author Univ. of Colorado Denver
159*> \author NAG Ltd.
160*
161*> \ingroup unmql
162*
163* =====================================================================
164 SUBROUTINE sormql( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
165 $ WORK, LWORK, INFO )
166*
167* -- LAPACK computational routine --
168* -- LAPACK is a software package provided by Univ. of Tennessee, --
169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*
171* .. Scalar Arguments ..
172 CHARACTER SIDE, TRANS
173 INTEGER INFO, K, LDA, LDC, LWORK, M, N
174* ..
175* .. Array Arguments ..
176 REAL A( LDA, * ), C( LDC, * ), TAU( * ),
177 $ work( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 INTEGER NBMAX, LDT, TSIZE
184 parameter( nbmax = 64, ldt = nbmax+1,
185 $ tsize = ldt*nbmax )
186* ..
187* .. Local Scalars ..
188 LOGICAL LEFT, LQUERY, NOTRAN
189 INTEGER I, I1, I2, I3, IB, IINFO, IWT, LDWORK, LWKOPT,
190 $ mi, nb, nbmin, ni, nq, nw
191* ..
192* .. External Functions ..
193 LOGICAL LSAME
194 INTEGER ILAENV
195 REAL SROUNDUP_LWORK
196 EXTERNAL lsame, ilaenv, sroundup_lwork
197* ..
198* .. External Subroutines ..
199 EXTERNAL slarfb, slarft, sorm2l, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC max, min
203* ..
204* .. Executable Statements ..
205*
206* Test the input arguments
207*
208 info = 0
209 left = lsame( side, 'L' )
210 notran = lsame( trans, 'N' )
211 lquery = ( lwork.EQ.-1 )
212*
213* NQ is the order of Q and NW is the minimum dimension of WORK
214*
215 IF( left ) THEN
216 nq = m
217 nw = max( 1, n )
218 ELSE
219 nq = n
220 nw = max( 1, m )
221 END IF
222 IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
223 info = -1
224 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
225 info = -2
226 ELSE IF( m.LT.0 ) THEN
227 info = -3
228 ELSE IF( n.LT.0 ) THEN
229 info = -4
230 ELSE IF( k.LT.0 .OR. k.GT.nq ) THEN
231 info = -5
232 ELSE IF( lda.LT.max( 1, nq ) ) THEN
233 info = -7
234 ELSE IF( ldc.LT.max( 1, m ) ) THEN
235 info = -10
236 ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
237 info = -12
238 END IF
239*
240 IF( info.EQ.0 ) THEN
241*
242* Compute the workspace requirements
243*
244 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
245 lwkopt = 1
246 ELSE
247 nb = min( nbmax, ilaenv( 1, 'SORMQL', side // trans, m,
248 $ n,
249 $ k, -1 ) )
250 lwkopt = nw*nb + tsize
251 END IF
252 work( 1 ) = sroundup_lwork(lwkopt)
253 END IF
254*
255 IF( info.NE.0 ) THEN
256 CALL xerbla( 'SORMQL', -info )
257 RETURN
258 ELSE IF( lquery ) THEN
259 RETURN
260 END IF
261*
262* Quick return if possible
263*
264 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
265 RETURN
266 END IF
267*
268 nbmin = 2
269 ldwork = nw
270 IF( nb.GT.1 .AND. nb.LT.k ) THEN
271 IF( lwork.LT.lwkopt ) THEN
272 nb = (lwork-tsize) / ldwork
273 nbmin = max( 2, ilaenv( 2, 'SORMQL', side // trans, m, n,
274 $ k,
275 $ -1 ) )
276 END IF
277 END IF
278*
279 IF( nb.LT.nbmin .OR. nb.GE.k ) THEN
280*
281* Use unblocked code
282*
283 CALL sorm2l( side, trans, m, n, k, a, lda, tau, c, ldc,
284 $ work,
285 $ iinfo )
286 ELSE
287*
288* Use blocked code
289*
290 iwt = 1 + nw*nb
291 IF( ( left .AND. notran ) .OR.
292 $ ( .NOT.left .AND. .NOT.notran ) ) THEN
293 i1 = 1
294 i2 = k
295 i3 = nb
296 ELSE
297 i1 = ( ( k-1 ) / nb )*nb + 1
298 i2 = 1
299 i3 = -nb
300 END IF
301*
302 IF( left ) THEN
303 ni = n
304 ELSE
305 mi = m
306 END IF
307*
308 DO 10 i = i1, i2, i3
309 ib = min( nb, k-i+1 )
310*
311* Form the triangular factor of the block reflector
312* H = H(i+ib-1) . . . H(i+1) H(i)
313*
314 CALL slarft( 'Backward', 'Columnwise', nq-k+i+ib-1, ib,
315 $ a( 1, i ), lda, tau( i ), work( iwt ), ldt )
316 IF( left ) THEN
317*
318* H or H**T is applied to C(1:m-k+i+ib-1,1:n)
319*
320 mi = m - k + i + ib - 1
321 ELSE
322*
323* H or H**T is applied to C(1:m,1:n-k+i+ib-1)
324*
325 ni = n - k + i + ib - 1
326 END IF
327*
328* Apply H or H**T
329*
330 CALL slarfb( side, trans, 'Backward', 'Columnwise', mi,
331 $ ni,
332 $ ib, a( 1, i ), lda, work( iwt ), ldt, c, ldc,
333 $ work, ldwork )
334 10 CONTINUE
335 END IF
336 work( 1 ) = sroundup_lwork(lwkopt)
337 RETURN
338*
339* End of SORMQL
340*
341 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine slarfb(side, trans, direct, storev, m, n, k, v, ldv, t, ldt, c, ldc, work, ldwork)
SLARFB applies a block reflector or its transpose to a general rectangular matrix.
Definition slarfb.f:195
recursive subroutine slarft(direct, storev, n, k, v, ldv, tau, t, ldt)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition slarft.f:162
subroutine sorm2l(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
SORM2L multiplies a general matrix by the orthogonal matrix from a QL factorization determined by sge...
Definition sorm2l.f:157
subroutine sormql(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQL
Definition sormql.f:166