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

◆ 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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*
393* Use calls to the subroutine DLAMC3 to enforce the
394* parentheses (x+y)+z. The goal is to prevent
395* optimizing compilers from doing x+(y+z).
396*
397 work( i ) = poles( i, 2 )*z( i ) /
398 $ ( dlamc3( poles( i, 2 ), dsigj )-
399 $ diflj ) / ( poles( i, 2 )+dj )
400 END IF
401 30 CONTINUE
402 DO 40 i = j + 1, k
403 IF( ( z( i ).EQ.zero ) .OR.
404 $ ( poles( i, 2 ).EQ.zero ) ) THEN
405 work( i ) = zero
406 ELSE
407 work( i ) = poles( i, 2 )*z( i ) /
408 $ ( dlamc3( poles( i, 2 ), dsigjp )+
409 $ difrj ) / ( poles( i, 2 )+dj )
410 END IF
411 40 CONTINUE
412 work( 1 ) = negone
413 temp = dnrm2( k, work, 1 )
414 CALL dgemv( 'T', k, nrhs, one, bx, ldbx, work, 1, zero,
415 $ b( j, 1 ), ldb )
416 CALL dlascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
417 $ ldb, info )
418 50 CONTINUE
419 END IF
420*
421* Move the deflated rows of BX to B also.
422*
423 IF( k.LT.max( m, n ) )
424 $ CALL dlacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
425 $ b( k+1, 1 ), ldb )
426 ELSE
427*
428* Apply back the right orthogonal transformations.
429*
430* Step (1R): apply back the new right singular vector matrix
431* to B.
432*
433 IF( k.EQ.1 ) THEN
434 CALL dcopy( nrhs, b, ldb, bx, ldbx )
435 ELSE
436 DO 80 j = 1, k
437 dsigj = poles( j, 2 )
438 IF( z( j ).EQ.zero ) THEN
439 work( j ) = zero
440 ELSE
441 work( j ) = -z( j ) / difl( j ) /
442 $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
443 END IF
444 DO 60 i = 1, j - 1
445 IF( z( j ).EQ.zero ) THEN
446 work( i ) = zero
447 ELSE
448*
449* Use calls to the subroutine DLAMC3 to enforce the
450* parentheses (x+y)+z. The goal is to prevent
451* optimizing compilers from doing x+(y+z).
452*
453 work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i+1,
454 $ 2 ) )-difr( i, 1 ) ) /
455 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
456 END IF
457 60 CONTINUE
458 DO 70 i = j + 1, k
459 IF( z( j ).EQ.zero ) THEN
460 work( i ) = zero
461 ELSE
462 work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
463 $ 2 ) )-difl( i ) ) /
464 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
465 END IF
466 70 CONTINUE
467 CALL dgemv( 'T', k, nrhs, one, b, ldb, work, 1, zero,
468 $ bx( j, 1 ), ldbx )
469 80 CONTINUE
470 END IF
471*
472* Step (2R): if SQRE = 1, apply back the rotation that is
473* related to the right null space of the subproblem.
474*
475 IF( sqre.EQ.1 ) THEN
476 CALL dcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
477 CALL drot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
478 END IF
479 IF( k.LT.max( m, n ) )
480 $ CALL dlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
481 $ ldbx )
482*
483* Step (3R): permute rows of B.
484*
485 CALL dcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
486 IF( sqre.EQ.1 ) THEN
487 CALL dcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
488 END IF
489 DO 90 i = 2, n
490 CALL dcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
491 90 CONTINUE
492*
493* Step (4R): apply back the Givens rotations performed.
494*
495 DO 100 i = givptr, 1, -1
496 CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
497 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
498 $ -givnum( i, 1 ) )
499 100 CONTINUE
500 END IF
501*
502 RETURN
503*
504* End of DLALS0
505*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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:103
dlamc3
double precision function dlamc3(a, b)
DLAMC3
Definition dlamch.f:172
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:143
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: