LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
clals0.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine clals0 (ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO)
 CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. Used by sgelsd. More...
 

Function/Subroutine Documentation

subroutine clals0 ( integer  ICOMPQ,
integer  NL,
integer  NR,
integer  SQRE,
integer  NRHS,
complex, dimension( ldb, * )  B,
integer  LDB,
complex, dimension( ldbx, * )  BX,
integer  LDBX,
integer, dimension( * )  PERM,
integer  GIVPTR,
integer, dimension( ldgcol, * )  GIVCOL,
integer  LDGCOL,
real, dimension( ldgnum, * )  GIVNUM,
integer  LDGNUM,
real, dimension( ldgnum, * )  POLES,
real, dimension( * )  DIFL,
real, dimension( ldgnum, * )  DIFR,
real, dimension( * )  Z,
integer  K,
real  C,
real  S,
real, dimension( * )  RWORK,
integer  INFO 
)

CLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. Used by sgelsd.

Download CLALS0 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CLALS0 applies back the multiplying factors of either the left or the
 right singular vector matrix of a diagonal matrix appended by a row
 to the right hand side matrix B in solving the least squares problem
 using the divide-and-conquer SVD approach.

 For the left singular vector matrix, three types of orthogonal
 matrices are involved:

 (1L) Givens rotations: the number of such rotations is GIVPTR; the
      pairs of columns/rows they were applied to are stored in GIVCOL;
      and the C- and S-values of these rotations are stored in GIVNUM.

 (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
      row, and for J=2:N, PERM(J)-th row of B is to be moved to the
      J-th row.

 (3L) The left singular vector matrix of the remaining matrix.

 For the right singular vector matrix, four types of orthogonal
 matrices are involved:

 (1R) The right singular vector matrix of the remaining matrix.

 (2R) If SQRE = 1, one extra Givens rotation to generate the right
      null space.

 (3R) The inverse transformation of (2L).

 (4R) The inverse transformation of (1L).
Parameters
[in]ICOMPQ
          ICOMPQ is INTEGER
         Specifies whether singular vectors are to be computed in
         factored form:
         = 0: Left singular vector matrix.
         = 1: Right singular vector matrix.
[in]NL
          NL is INTEGER
         The row dimension of the upper block. NL >= 1.
[in]NR
          NR is INTEGER
         The row dimension of the lower block. NR >= 1.
[in]SQRE
          SQRE is INTEGER
         = 0: the lower block is an NR-by-NR square matrix.
         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

         The bidiagonal matrix has row dimension N = NL + NR + 1,
         and column dimension M = N + SQRE.
[in]NRHS
          NRHS is INTEGER
         The number of columns of B and BX. NRHS must be at least 1.
[in,out]B
          B is COMPLEX array, dimension ( LDB, NRHS )
         On input, B contains the right hand sides of the least
         squares problem in rows 1 through M. On output, B contains
         the solution X in rows 1 through N.
[in]LDB
          LDB is INTEGER
         The leading dimension of B. LDB must be at least
         max(1,MAX( M, N ) ).
[out]BX
          BX is COMPLEX array, dimension ( LDBX, NRHS )
[in]LDBX
          LDBX is INTEGER
         The leading dimension of BX.
[in]PERM
          PERM is INTEGER array, dimension ( N )
         The permutations (from deflation and sorting) applied
         to the two blocks.
[in]GIVPTR
          GIVPTR is INTEGER
         The number of Givens rotations which took place in this
         subproblem.
[in]GIVCOL
          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
         Each pair of numbers indicates a pair of rows/columns
         involved in a Givens rotation.
[in]LDGCOL
          LDGCOL is INTEGER
         The leading dimension of GIVCOL, must be at least N.
[in]GIVNUM
          GIVNUM is REAL array, dimension ( LDGNUM, 2 )
         Each number indicates the C or S value used in the
         corresponding Givens rotation.
[in]LDGNUM
          LDGNUM is INTEGER
         The leading dimension of arrays DIFR, POLES and
         GIVNUM, must be at least K.
[in]POLES
          POLES is REAL array, dimension ( LDGNUM, 2 )
         On entry, POLES(1:K, 1) contains the new singular
         values obtained from solving the secular equation, and
         POLES(1:K, 2) is an array containing the poles in the secular
         equation.
[in]DIFL
          DIFL is REAL array, dimension ( K ).
         On entry, DIFL(I) is the distance between I-th updated
         (undeflated) singular value and the I-th (undeflated) old
         singular value.
[in]DIFR
          DIFR is REAL array, dimension ( LDGNUM, 2 ).
         On entry, DIFR(I, 1) contains the distances between I-th
         updated (undeflated) singular value and the I+1-th
         (undeflated) old singular value. And DIFR(I, 2) is the
         normalizing factor for the I-th right singular vector.
[in]Z
          Z is REAL array, dimension ( K )
         Contain the components of the deflation-adjusted updating row
         vector.
[in]K
          K is INTEGER
         Contains the dimension of the non-deflated matrix,
         This is the order of the related secular equation. 1 <= K <=N.
[in]C
          C is REAL
         C contains garbage if SQRE =0 and the C-value of a Givens
         rotation related to the right null space if SQRE = 1.
[in]S
          S is REAL
         S contains garbage if SQRE =0 and the S-value of a Givens
         rotation related to the right null space if SQRE = 1.
[out]RWORK
          RWORK is REAL array, dimension
         ( K*(1+NRHS) + 2*NRHS )
[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.
Date
September 2012
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 272 of file clals0.f.

272 *
273 * -- LAPACK computational routine (version 3.4.2) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * September 2012
277 *
278 * .. Scalar Arguments ..
279  INTEGER givptr, icompq, info, k, ldb, ldbx, ldgcol,
280  $ ldgnum, nl, nr, nrhs, sqre
281  REAL c, s
282 * ..
283 * .. Array Arguments ..
284  INTEGER givcol( ldgcol, * ), perm( * )
285  REAL difl( * ), difr( ldgnum, * ),
286  $ givnum( ldgnum, * ), poles( ldgnum, * ),
287  $ rwork( * ), z( * )
288  COMPLEX b( ldb, * ), bx( ldbx, * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  REAL one, zero, negone
295  parameter( one = 1.0e0, zero = 0.0e0, negone = -1.0e0 )
296 * ..
297 * .. Local Scalars ..
298  INTEGER i, j, jcol, jrow, m, n, nlp1
299  REAL diflj, difrj, dj, dsigj, dsigjp, temp
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL ccopy, clacpy, clascl, csrot, csscal, sgemv,
303  $ xerbla
304 * ..
305 * .. External Functions ..
306  REAL slamc3, snrm2
307  EXTERNAL slamc3, snrm2
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC aimag, cmplx, max, real
311 * ..
312 * .. Executable Statements ..
313 *
314 * Test the input parameters.
315 *
316  info = 0
317 *
318  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
319  info = -1
320  ELSE IF( nl.LT.1 ) THEN
321  info = -2
322  ELSE IF( nr.LT.1 ) THEN
323  info = -3
324  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
325  info = -4
326  END IF
327 *
328  n = nl + nr + 1
329 *
330  IF( nrhs.LT.1 ) THEN
331  info = -5
332  ELSE IF( ldb.LT.n ) THEN
333  info = -7
334  ELSE IF( ldbx.LT.n ) THEN
335  info = -9
336  ELSE IF( givptr.LT.0 ) THEN
337  info = -11
338  ELSE IF( ldgcol.LT.n ) THEN
339  info = -13
340  ELSE IF( ldgnum.LT.n ) THEN
341  info = -15
342  ELSE IF( k.LT.1 ) THEN
343  info = -20
344  END IF
345  IF( info.NE.0 ) THEN
346  CALL xerbla( 'CLALS0', -info )
347  RETURN
348  END IF
349 *
350  m = n + sqre
351  nlp1 = nl + 1
352 *
353  IF( icompq.EQ.0 ) THEN
354 *
355 * Apply back orthogonal transformations from the left.
356 *
357 * Step (1L): apply back the Givens rotations performed.
358 *
359  DO 10 i = 1, givptr
360  CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
361  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
362  $ givnum( i, 1 ) )
363  10 CONTINUE
364 *
365 * Step (2L): permute rows of B.
366 *
367  CALL ccopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
368  DO 20 i = 2, n
369  CALL ccopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
370  20 CONTINUE
371 *
372 * Step (3L): apply the inverse of the left singular vector
373 * matrix to BX.
374 *
375  IF( k.EQ.1 ) THEN
376  CALL ccopy( nrhs, bx, ldbx, b, ldb )
377  IF( z( 1 ).LT.zero ) THEN
378  CALL csscal( nrhs, negone, b, ldb )
379  END IF
380  ELSE
381  DO 100 j = 1, k
382  diflj = difl( j )
383  dj = poles( j, 1 )
384  dsigj = -poles( j, 2 )
385  IF( j.LT.k ) THEN
386  difrj = -difr( j, 1 )
387  dsigjp = -poles( j+1, 2 )
388  END IF
389  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
390  $ THEN
391  rwork( j ) = zero
392  ELSE
393  rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
394  $ ( poles( j, 2 )+dj )
395  END IF
396  DO 30 i = 1, j - 1
397  IF( ( z( i ).EQ.zero ) .OR.
398  $ ( poles( i, 2 ).EQ.zero ) ) THEN
399  rwork( i ) = zero
400  ELSE
401  rwork( i ) = poles( i, 2 )*z( i ) /
402  $ ( slamc3( poles( i, 2 ), dsigj )-
403  $ diflj ) / ( poles( i, 2 )+dj )
404  END IF
405  30 CONTINUE
406  DO 40 i = j + 1, k
407  IF( ( z( i ).EQ.zero ) .OR.
408  $ ( poles( i, 2 ).EQ.zero ) ) THEN
409  rwork( i ) = zero
410  ELSE
411  rwork( i ) = poles( i, 2 )*z( i ) /
412  $ ( slamc3( poles( i, 2 ), dsigjp )+
413  $ difrj ) / ( poles( i, 2 )+dj )
414  END IF
415  40 CONTINUE
416  rwork( 1 ) = negone
417  temp = snrm2( k, rwork, 1 )
418 *
419 * Since B and BX are complex, the following call to SGEMV
420 * is performed in two steps (real and imaginary parts).
421 *
422 * CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
423 * $ B( J, 1 ), LDB )
424 *
425  i = k + nrhs*2
426  DO 60 jcol = 1, nrhs
427  DO 50 jrow = 1, k
428  i = i + 1
429  rwork( i ) = REAL( BX( JROW, JCOL ) )
430  50 CONTINUE
431  60 CONTINUE
432  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
433  $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
434  i = k + nrhs*2
435  DO 80 jcol = 1, nrhs
436  DO 70 jrow = 1, k
437  i = i + 1
438  rwork( i ) = aimag( bx( jrow, jcol ) )
439  70 CONTINUE
440  80 CONTINUE
441  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
442  $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
443  DO 90 jcol = 1, nrhs
444  b( j, jcol ) = cmplx( rwork( jcol+k ),
445  $ rwork( jcol+k+nrhs ) )
446  90 CONTINUE
447  CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
448  $ ldb, info )
449  100 CONTINUE
450  END IF
451 *
452 * Move the deflated rows of BX to B also.
453 *
454  IF( k.LT.max( m, n ) )
455  $ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
456  $ b( k+1, 1 ), ldb )
457  ELSE
458 *
459 * Apply back the right orthogonal transformations.
460 *
461 * Step (1R): apply back the new right singular vector matrix
462 * to B.
463 *
464  IF( k.EQ.1 ) THEN
465  CALL ccopy( nrhs, b, ldb, bx, ldbx )
466  ELSE
467  DO 180 j = 1, k
468  dsigj = poles( j, 2 )
469  IF( z( j ).EQ.zero ) THEN
470  rwork( j ) = zero
471  ELSE
472  rwork( j ) = -z( j ) / difl( j ) /
473  $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
474  END IF
475  DO 110 i = 1, j - 1
476  IF( z( j ).EQ.zero ) THEN
477  rwork( i ) = zero
478  ELSE
479  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
480  $ 2 ) )-difr( i, 1 ) ) /
481  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
482  END IF
483  110 CONTINUE
484  DO 120 i = j + 1, k
485  IF( z( j ).EQ.zero ) THEN
486  rwork( i ) = zero
487  ELSE
488  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
489  $ 2 ) )-difl( i ) ) /
490  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
491  END IF
492  120 CONTINUE
493 *
494 * Since B and BX are complex, the following call to SGEMV
495 * is performed in two steps (real and imaginary parts).
496 *
497 * CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
498 * $ BX( J, 1 ), LDBX )
499 *
500  i = k + nrhs*2
501  DO 140 jcol = 1, nrhs
502  DO 130 jrow = 1, k
503  i = i + 1
504  rwork( i ) = REAL( B( JROW, JCOL ) )
505  130 CONTINUE
506  140 CONTINUE
507  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
508  $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
509  i = k + nrhs*2
510  DO 160 jcol = 1, nrhs
511  DO 150 jrow = 1, k
512  i = i + 1
513  rwork( i ) = aimag( b( jrow, jcol ) )
514  150 CONTINUE
515  160 CONTINUE
516  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
517  $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
518  DO 170 jcol = 1, nrhs
519  bx( j, jcol ) = cmplx( rwork( jcol+k ),
520  $ rwork( jcol+k+nrhs ) )
521  170 CONTINUE
522  180 CONTINUE
523  END IF
524 *
525 * Step (2R): if SQRE = 1, apply back the rotation that is
526 * related to the right null space of the subproblem.
527 *
528  IF( sqre.EQ.1 ) THEN
529  CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
530  CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
531  END IF
532  IF( k.LT.max( m, n ) )
533  $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
534  $ bx( k+1, 1 ), ldbx )
535 *
536 * Step (3R): permute rows of B.
537 *
538  CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
539  IF( sqre.EQ.1 ) THEN
540  CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
541  END IF
542  DO 190 i = 2, n
543  CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
544  190 CONTINUE
545 *
546 * Step (4R): apply back the Givens rotations performed.
547 *
548  DO 200 i = givptr, 1, -1
549  CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
550  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
551  $ -givnum( i, 1 ) )
552  200 CONTINUE
553  END IF
554 *
555  RETURN
556 *
557 * End of CLALS0
558 *
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:100
real function slamc3(A, B)
SLAMC3
Definition: slamch.f:172
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:141
real function snrm2(N, X, INCX)
SNRM2
Definition: snrm2.f:56
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: