LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dlals0 ( integer  ICOMPQ,
integer  NL,
integer  NR,
integer  SQRE,
integer  NRHS,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldbx, * )  BX,
integer  LDBX,
integer, dimension( * )  PERM,
integer  GIVPTR,
integer, dimension( ldgcol, * )  GIVCOL,
integer  LDGCOL,
double precision, dimension( ldgnum, * )  GIVNUM,
integer  LDGNUM,
double precision, dimension( ldgnum, * )  POLES,
double precision, dimension( * )  DIFL,
double precision, dimension( ldgnum, * )  DIFR,
double precision, dimension( * )  Z,
integer  K,
double precision  C,
double precision  S,
double precision, dimension( * )  WORK,
integer  INFO 
)

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

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

Purpose: 
 DLALS0 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
         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 DOUBLE PRECISION
         S contains garbage if SQRE =0 and the S-value of a Givens
         rotation related to the right null space if SQRE = 1.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension ( K )
[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
November 2015
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 270 of file dlals0.f.

270 *
271 * -- LAPACK computational routine (version 3.6.0) --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 * November 2015
275 *
276 * .. Scalar Arguments ..
277  INTEGER givptr, icompq, info, k, ldb, ldbx, ldgcol,
278  $ ldgnum, nl, nr, nrhs, sqre
279  DOUBLE PRECISION c, s
280 * ..
281 * .. Array Arguments ..
282  INTEGER givcol( ldgcol, * ), perm( * )
283  DOUBLE PRECISION b( ldb, * ), bx( ldbx, * ), difl( * ),
284  $ difr( ldgnum, * ), givnum( ldgnum, * ),
285  $ poles( ldgnum, * ), work( * ), z( * )
286 * ..
287 *
288 * =====================================================================
289 *
290 * .. Parameters ..
291  DOUBLE PRECISION one, zero, negone
292  parameter ( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
293 * ..
294 * .. Local Scalars ..
295  INTEGER i, j, m, n, nlp1
296  DOUBLE PRECISION diflj, difrj, dj, dsigj, dsigjp, temp
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL dcopy, dgemv, dlacpy, dlascl, drot, dscal,
300  $ xerbla
301 * ..
302 * .. External Functions ..
303  DOUBLE PRECISION dlamc3, dnrm2
304  EXTERNAL dlamc3, dnrm2
305 * ..
306 * .. Intrinsic Functions ..
307  INTRINSIC max
308 * ..
309 * .. Executable Statements ..
310 *
311 * Test the input parameters.
312 *
313  info = 0
314  n = nl + nr + 1
315 *
316  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
317  info = -1
318  ELSE IF( nl.LT.1 ) THEN
319  info = -2
320  ELSE IF( nr.LT.1 ) THEN
321  info = -3
322  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
323  info = -4
324  ELSE IF( nrhs.LT.1 ) THEN
325  info = -5
326  ELSE IF( ldb.LT.n ) THEN
327  info = -7
328  ELSE IF( ldbx.LT.n ) THEN
329  info = -9
330  ELSE IF( givptr.LT.0 ) THEN
331  info = -11
332  ELSE IF( ldgcol.LT.n ) THEN
333  info = -13
334  ELSE IF( ldgnum.LT.n ) THEN
335  info = -15
336  ELSE IF( k.LT.1 ) THEN
337  info = -20
338  END IF
339  IF( info.NE.0 ) THEN
340  CALL xerbla( 'DLALS0', -info )
341  RETURN
342  END IF
343 *
344  m = n + sqre
345  nlp1 = nl + 1
346 *
347  IF( icompq.EQ.0 ) THEN
348 *
349 * Apply back orthogonal transformations from the left.
350 *
351 * Step (1L): apply back the Givens rotations performed.
352 *
353  DO 10 i = 1, givptr
354  CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
355  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
356  $ givnum( i, 1 ) )
357  10 CONTINUE
358 *
359 * Step (2L): permute rows of B.
360 *
361  CALL dcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
362  DO 20 i = 2, n
363  CALL dcopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
364  20 CONTINUE
365 *
366 * Step (3L): apply the inverse of the left singular vector
367 * matrix to BX.
368 *
369  IF( k.EQ.1 ) THEN
370  CALL dcopy( nrhs, bx, ldbx, b, ldb )
371  IF( z( 1 ).LT.zero ) THEN
372  CALL dscal( nrhs, negone, b, ldb )
373  END IF
374  ELSE
375  DO 50 j = 1, k
376  diflj = difl( j )
377  dj = poles( j, 1 )
378  dsigj = -poles( j, 2 )
379  IF( j.LT.k ) THEN
380  difrj = -difr( j, 1 )
381  dsigjp = -poles( j+1, 2 )
382  END IF
383  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
384  $ THEN
385  work( j ) = zero
386  ELSE
387  work( j ) = -poles( j, 2 )*z( j ) / diflj /
388  $ ( poles( j, 2 )+dj )
389  END IF
390  DO 30 i = 1, j - 1
391  IF( ( z( i ).EQ.zero ) .OR.
392  $ ( poles( i, 2 ).EQ.zero ) ) THEN
393  work( i ) = zero
394  ELSE
395  work( i ) = poles( i, 2 )*z( i ) /
396  $ ( dlamc3( poles( i, 2 ), dsigj )-
397  $ diflj ) / ( poles( i, 2 )+dj )
398  END IF
399  30 CONTINUE
400  DO 40 i = j + 1, k
401  IF( ( z( i ).EQ.zero ) .OR.
402  $ ( poles( i, 2 ).EQ.zero ) ) THEN
403  work( i ) = zero
404  ELSE
405  work( i ) = poles( i, 2 )*z( i ) /
406  $ ( dlamc3( poles( i, 2 ), dsigjp )+
407  $ difrj ) / ( poles( i, 2 )+dj )
408  END IF
409  40 CONTINUE
410  work( 1 ) = negone
411  temp = dnrm2( k, work, 1 )
412  CALL dgemv( 'T', k, nrhs, one, bx, ldbx, work, 1, zero,
413  $ b( j, 1 ), ldb )
414  CALL dlascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
415  $ ldb, info )
416  50 CONTINUE
417  END IF
418 *
419 * Move the deflated rows of BX to B also.
420 *
421  IF( k.LT.max( m, n ) )
422  $ CALL dlacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
423  $ b( k+1, 1 ), ldb )
424  ELSE
425 *
426 * Apply back the right orthogonal transformations.
427 *
428 * Step (1R): apply back the new right singular vector matrix
429 * to B.
430 *
431  IF( k.EQ.1 ) THEN
432  CALL dcopy( nrhs, b, ldb, bx, ldbx )
433  ELSE
434  DO 80 j = 1, k
435  dsigj = poles( j, 2 )
436  IF( z( j ).EQ.zero ) THEN
437  work( j ) = zero
438  ELSE
439  work( j ) = -z( j ) / difl( j ) /
440  $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
441  END IF
442  DO 60 i = 1, j - 1
443  IF( z( j ).EQ.zero ) THEN
444  work( i ) = zero
445  ELSE
446  work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i+1,
447  $ 2 ) )-difr( i, 1 ) ) /
448  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
449  END IF
450  60 CONTINUE
451  DO 70 i = j + 1, k
452  IF( z( j ).EQ.zero ) THEN
453  work( i ) = zero
454  ELSE
455  work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
456  $ 2 ) )-difl( i ) ) /
457  $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
458  END IF
459  70 CONTINUE
460  CALL dgemv( 'T', k, nrhs, one, b, ldb, work, 1, zero,
461  $ bx( j, 1 ), ldbx )
462  80 CONTINUE
463  END IF
464 *
465 * Step (2R): if SQRE = 1, apply back the rotation that is
466 * related to the right null space of the subproblem.
467 *
468  IF( sqre.EQ.1 ) THEN
469  CALL dcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
470  CALL drot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
471  END IF
472  IF( k.LT.max( m, n ) )
473  $ CALL dlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
474  $ ldbx )
475 *
476 * Step (3R): permute rows of B.
477 *
478  CALL dcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
479  IF( sqre.EQ.1 ) THEN
480  CALL dcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
481  END IF
482  DO 90 i = 2, n
483  CALL dcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
484  90 CONTINUE
485 *
486 * Step (4R): apply back the Givens rotations performed.
487 *
488  DO 100 i = givptr, 1, -1
489  CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
490  $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
491  $ -givnum( i, 1 ) )
492  100 CONTINUE
493  END IF
494 *
495  RETURN
496 *
497 * End of DLALS0
498 *
dlamc3
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
dlascl
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dnrm2
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function: