LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slamtsqr()

subroutine slamtsqr ( character  side,
character  trans,
integer  m,
integer  n,
integer  k,
integer  mb,
integer  nb,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( ldt, * )  t,
integer  ldt,
real, dimension(ldc, * )  c,
integer  ldc,
real, dimension( * )  work,
integer  lwork,
integer  info 
)

SLAMTSQR

Purpose:
      SLAMTSQR overwrites the general real M-by-N matrix C with


                 SIDE = 'L'     SIDE = 'R'
 TRANS = 'N':      Q * C          C * Q
 TRANS = 'T':      Q**T * C       C * Q**T
      where Q is a real orthogonal matrix defined as the product
      of blocked elementary reflectors computed by tall skinny
      QR factorization (SLATSQR)
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.
[in]TRANS
          TRANS is CHARACTER*1
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >=0.
[in]N
          N is INTEGER
          The number of columns of the matrix C. N >= 0.
[in]K
          K is INTEGER
          The number of elementary reflectors whose product defines
          the matrix Q. M >= K >= 0;
[in]MB
          MB is INTEGER
          The block size to be used in the blocked QR.
          MB > N. (must be the same as SLATSQR)
[in]NB
          NB is INTEGER
          The column block size to be used in the blocked QR.
          N >= NB >= 1.
[in]A
          A is REAL array, dimension (LDA,K)
          The i-th column must contain the vector which defines the
          blockedelementary reflector H(i), for i = 1,2,...,k, as
          returned by SLATSQR in the first k columns of
          its array argument A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.
          If SIDE = 'L', LDA >= max(1,M);
          if SIDE = 'R', LDA >= max(1,N).
[in]T
          T is REAL array, dimension
          ( N * Number of blocks(CEIL(M-K/MB-K)),
          The blocked upper triangular block reflectors stored in compact form
          as a sequence of upper triangular blocks.  See below
          for further details.
[in]LDT
          LDT is INTEGER
          The leading dimension of the array T.  LDT >= NB.
[in,out]C
          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).
[out]WORK
         (workspace) REAL array, dimension (MAX(1,LWORK))
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.

          If SIDE = 'L', LWORK >= max(1,N)*NB;
          if SIDE = 'R', LWORK >= max(1,MB)*NB.
          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
 Tall-Skinny QR (TSQR) performs QR by a sequence of orthogonal transformations,
 representing Q as a product of other orthogonal matrices
   Q = Q(1) * Q(2) * . . . * Q(k)
 where each Q(i) zeros out subdiagonal entries of a block of MB rows of A:
   Q(1) zeros out the subdiagonal entries of rows 1:MB of A
   Q(2) zeros out the bottom MB-N rows of rows [1:N,MB+1:2*MB-N] of A
   Q(3) zeros out the bottom MB-N rows of rows [1:N,2*MB-N+1:3*MB-2*N] of A
   . . .

 Q(1) is computed by GEQRT, which represents Q(1) by Householder vectors
 stored under the diagonal of rows 1:MB of A, and by upper triangular
 block reflectors, stored in array T(1:LDT,1:N).
 For more information see Further Details in GEQRT.

 Q(i) for i>1 is computed by TPQRT, which represents Q(i) by Householder vectors
 stored in rows [(i-1)*(MB-N)+N+1:i*(MB-N)+N] of A, and by upper triangular
 block reflectors, stored in array T(1:LDT,(i-1)*N+1:i*N).
 The last Q(k) may use fewer rows.
 For more information see Further Details in TPQRT.

 For more details of the overall algorithm, see the description of
 Sequential TSQR in Section 2.2 of [1].

 [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,”
     J. Demmel, L. Grigori, M. Hoemmen, J. Langou,
     SIAM J. Sci. Comput, vol. 34, no. 1, 2012

Definition at line 197 of file slamtsqr.f.

199*
200* -- LAPACK computational routine --
201* -- LAPACK is a software package provided by Univ. of Tennessee, --
202* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203*
204* .. Scalar Arguments ..
205 CHARACTER SIDE, TRANS
206 INTEGER INFO, LDA, M, N, K, MB, NB, LDT, LWORK, LDC
207* ..
208* .. Array Arguments ..
209 REAL A( LDA, * ), WORK( * ), C(LDC, * ),
210 $ T( LDT, * )
211* ..
212*
213* =====================================================================
214*
215* ..
216* .. Local Scalars ..
217 LOGICAL LEFT, RIGHT, TRAN, NOTRAN, LQUERY
218 INTEGER I, II, KK, LW, CTR, Q
219* ..
220* .. External Functions ..
221 LOGICAL LSAME
222 EXTERNAL lsame
223* .. External Subroutines ..
224 EXTERNAL sgemqrt, stpmqrt, xerbla
225* ..
226* .. Executable Statements ..
227*
228* Test the input arguments
229*
230 lquery = lwork.LT.0
231 notran = lsame( trans, 'N' )
232 tran = lsame( trans, 'T' )
233 left = lsame( side, 'L' )
234 right = lsame( side, 'R' )
235 IF (left) THEN
236 lw = n * nb
237 q = m
238 ELSE
239 lw = mb * nb
240 q = n
241 END IF
242*
243 info = 0
244 IF( .NOT.left .AND. .NOT.right ) THEN
245 info = -1
246 ELSE IF( .NOT.tran .AND. .NOT.notran ) THEN
247 info = -2
248 ELSE IF( m.LT.k ) THEN
249 info = -3
250 ELSE IF( n.LT.0 ) THEN
251 info = -4
252 ELSE IF( k.LT.0 ) THEN
253 info = -5
254 ELSE IF( k.LT.nb .OR. nb.LT.1 ) THEN
255 info = -7
256 ELSE IF( lda.LT.max( 1, q ) ) THEN
257 info = -9
258 ELSE IF( ldt.LT.max( 1, nb) ) THEN
259 info = -11
260 ELSE IF( ldc.LT.max( 1, m ) ) THEN
261 info = -13
262 ELSE IF(( lwork.LT.max(1,lw)).AND.(.NOT.lquery)) THEN
263 info = -15
264 END IF
265*
266* Determine the block size if it is tall skinny or short and wide
267*
268 IF( info.EQ.0) THEN
269 work(1) = lw
270 END IF
271*
272 IF( info.NE.0 ) THEN
273 CALL xerbla( 'SLAMTSQR', -info )
274 RETURN
275 ELSE IF (lquery) THEN
276 RETURN
277 END IF
278*
279* Quick return if possible
280*
281 IF( min(m,n,k).EQ.0 ) THEN
282 RETURN
283 END IF
284*
285 IF((mb.LE.k).OR.(mb.GE.max(m,n,k))) THEN
286 CALL sgemqrt( side, trans, m, n, k, nb, a, lda,
287 $ t, ldt, c, ldc, work, info)
288 RETURN
289 END IF
290*
291 IF(left.AND.notran) THEN
292*
293* Multiply Q to the last block of C
294*
295 kk = mod((m-k),(mb-k))
296 ctr = (m-k)/(mb-k)
297 IF (kk.GT.0) THEN
298 ii=m-kk+1
299 CALL stpmqrt('L','N',kk , n, k, 0, nb, a(ii,1), lda,
300 $ t(1,ctr*k+1),ldt , c(1,1), ldc,
301 $ c(ii,1), ldc, work, info )
302 ELSE
303 ii=m+1
304 END IF
305*
306 DO i=ii-(mb-k),mb+1,-(mb-k)
307*
308* Multiply Q to the current block of C (I:I+MB,1:N)
309*
310 ctr = ctr - 1
311 CALL stpmqrt('L','N',mb-k , n, k, 0,nb, a(i,1), lda,
312 $ t(1, ctr * k + 1), ldt, c(1,1), ldc,
313 $ c(i,1), ldc, work, info )
314*
315 END DO
316*
317* Multiply Q to the first block of C (1:MB,1:N)
318*
319 CALL sgemqrt('L','N',mb , n, k, nb, a(1,1), lda, t
320 $ ,ldt ,c(1,1), ldc, work, info )
321*
322 ELSE IF (left.AND.tran) THEN
323*
324* Multiply Q to the first block of C
325*
326 kk = mod((m-k),(mb-k))
327 ii=m-kk+1
328 ctr = 1
329 CALL sgemqrt('L','T',mb , n, k, nb, a(1,1), lda, t
330 $ ,ldt ,c(1,1), ldc, work, info )
331*
332 DO i=mb+1,ii-mb+k,(mb-k)
333*
334* Multiply Q to the current block of C (I:I+MB,1:N)
335*
336 CALL stpmqrt('L','T',mb-k , n, k, 0,nb, a(i,1), lda,
337 $ t(1,ctr * k + 1),ldt, c(1,1), ldc,
338 $ c(i,1), ldc, work, info )
339 ctr = ctr + 1
340*
341 END DO
342 IF(ii.LE.m) THEN
343*
344* Multiply Q to the last block of C
345*
346 CALL stpmqrt('L','T',kk , n, k, 0,nb, a(ii,1), lda,
347 $ t(1, ctr * k + 1), ldt, c(1,1), ldc,
348 $ c(ii,1), ldc, work, info )
349*
350 END IF
351*
352 ELSE IF(right.AND.tran) THEN
353*
354* Multiply Q to the last block of C
355*
356 kk = mod((n-k),(mb-k))
357 ctr = (n-k)/(mb-k)
358 IF (kk.GT.0) THEN
359 ii=n-kk+1
360 CALL stpmqrt('R','T',m , kk, k, 0, nb, a(ii,1), lda,
361 $ t(1, ctr * k + 1), ldt, c(1,1), ldc,
362 $ c(1,ii), ldc, work, info )
363 ELSE
364 ii=n+1
365 END IF
366*
367 DO i=ii-(mb-k),mb+1,-(mb-k)
368*
369* Multiply Q to the current block of C (1:M,I:I+MB)
370*
371 ctr = ctr - 1
372 CALL stpmqrt('R','T',m , mb-k, k, 0,nb, a(i,1), lda,
373 $ t(1, ctr * k + 1), ldt, c(1,1), ldc,
374 $ c(1,i), ldc, work, info )
375*
376 END DO
377*
378* Multiply Q to the first block of C (1:M,1:MB)
379*
380 CALL sgemqrt('R','T',m , mb, k, nb, a(1,1), lda, t
381 $ ,ldt ,c(1,1), ldc, work, info )
382*
383 ELSE IF (right.AND.notran) THEN
384*
385* Multiply Q to the first block of C
386*
387 kk = mod((n-k),(mb-k))
388 ii=n-kk+1
389 ctr = 1
390 CALL sgemqrt('R','N', m, mb , k, nb, a(1,1), lda, t
391 $ ,ldt ,c(1,1), ldc, work, info )
392*
393 DO i=mb+1,ii-mb+k,(mb-k)
394*
395* Multiply Q to the current block of C (1:M,I:I+MB)
396*
397 CALL stpmqrt('R','N', m, mb-k, k, 0,nb, a(i,1), lda,
398 $ t(1, ctr * k + 1),ldt, c(1,1), ldc,
399 $ c(1,i), ldc, work, info )
400 ctr = ctr + 1
401*
402 END DO
403 IF(ii.LE.n) THEN
404*
405* Multiply Q to the last block of C
406*
407 CALL stpmqrt('R','N', m, kk , k, 0,nb, a(ii,1), lda,
408 $ t(1, ctr * k + 1),ldt, c(1,1), ldc,
409 $ c(1,ii), ldc, work, info )
410*
411 END IF
412*
413 END IF
414*
415 work(1) = lw
416 RETURN
417*
418* End of SLAMTSQR
419*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
SGEMQRT
Definition sgemqrt.f:168
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine stpmqrt(side, trans, m, n, k, l, nb, v, ldv, t, ldt, a, lda, b, ldb, work, info)
STPMQRT
Definition stpmqrt.f:216
Here is the call graph for this function:
Here is the caller graph for this function: