 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ dlals0()

 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.```
Contributors:
Ming Gu and Ren-Cang Li, Computer Science Division, University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA

Definition at line 265 of file dlals0.f.

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