LAPACK 3.12.1
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 263 of file dlals0.f.

267*
268* -- LAPACK computational routine --
269* -- LAPACK is a software package provided by Univ. of Tennessee, --
270* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271*
272* .. Scalar Arguments ..
273 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
274 $ LDGNUM, NL, NR, NRHS, SQRE
275 DOUBLE PRECISION C, S
276* ..
277* .. Array Arguments ..
278 INTEGER GIVCOL( LDGCOL, * ), PERM( * )
279 DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ),
280 $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
281 $ POLES( LDGNUM, * ), WORK( * ), Z( * )
282* ..
283*
284* =====================================================================
285*
286* .. Parameters ..
287 DOUBLE PRECISION ONE, ZERO, NEGONE
288 parameter( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
289* ..
290* .. Local Scalars ..
291 INTEGER I, J, M, N, NLP1
292 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
293* ..
294* .. External Subroutines ..
295 EXTERNAL dcopy, dgemv, dlacpy, dlascl, drot,
296 $ 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 ),
361 $ ldbx )
362 20 CONTINUE
363*
364* Step (3L): apply the inverse of the left singular vector
365* matrix to BX.
366*
367 IF( k.EQ.1 ) THEN
368 CALL dcopy( nrhs, bx, ldbx, b, ldb )
369 IF( z( 1 ).LT.zero ) THEN
370 CALL dscal( nrhs, negone, b, ldb )
371 END IF
372 ELSE
373 DO 50 j = 1, k
374 diflj = difl( j )
375 dj = poles( j, 1 )
376 dsigj = -poles( j, 2 )
377 IF( j.LT.k ) THEN
378 difrj = -difr( j, 1 )
379 dsigjp = -poles( j+1, 2 )
380 END IF
381 IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
382 $ THEN
383 work( j ) = zero
384 ELSE
385 work( j ) = -poles( j, 2 )*z( j ) / diflj /
386 $ ( poles( j, 2 )+dj )
387 END IF
388 DO 30 i = 1, j - 1
389 IF( ( z( i ).EQ.zero ) .OR.
390 $ ( poles( i, 2 ).EQ.zero ) ) THEN
391 work( i ) = zero
392 ELSE
393*
394* Use calls to the subroutine DLAMC3 to enforce the
395* parentheses (x+y)+z. The goal is to prevent
396* optimizing compilers from doing x+(y+z).
397*
398 work( i ) = poles( i, 2 )*z( i ) /
399 $ ( dlamc3( poles( i, 2 ), dsigj )-
400 $ diflj ) / ( poles( i, 2 )+dj )
401 END IF
402 30 CONTINUE
403 DO 40 i = j + 1, k
404 IF( ( z( i ).EQ.zero ) .OR.
405 $ ( poles( i, 2 ).EQ.zero ) ) THEN
406 work( i ) = zero
407 ELSE
408 work( i ) = poles( i, 2 )*z( i ) /
409 $ ( dlamc3( poles( i, 2 ), dsigjp )+
410 $ difrj ) / ( poles( i, 2 )+dj )
411 END IF
412 40 CONTINUE
413 work( 1 ) = negone
414 temp = dnrm2( k, work, 1 )
415 CALL dgemv( 'T', k, nrhs, one, bx, ldbx, work, 1,
416 $ zero,
417 $ b( j, 1 ), ldb )
418 CALL dlascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
419 $ ldb, info )
420 50 CONTINUE
421 END IF
422*
423* Move the deflated rows of BX to B also.
424*
425 IF( k.LT.max( m, n ) )
426 $ CALL dlacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
427 $ b( k+1, 1 ), ldb )
428 ELSE
429*
430* Apply back the right orthogonal transformations.
431*
432* Step (1R): apply back the new right singular vector matrix
433* to B.
434*
435 IF( k.EQ.1 ) THEN
436 CALL dcopy( nrhs, b, ldb, bx, ldbx )
437 ELSE
438 DO 80 j = 1, k
439 dsigj = poles( j, 2 )
440 IF( z( j ).EQ.zero ) THEN
441 work( j ) = zero
442 ELSE
443 work( j ) = -z( j ) / difl( j ) /
444 $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
445 END IF
446 DO 60 i = 1, j - 1
447 IF( z( j ).EQ.zero ) THEN
448 work( i ) = zero
449 ELSE
450*
451* Use calls to the subroutine DLAMC3 to enforce the
452* parentheses (x+y)+z. The goal is to prevent
453* optimizing compilers from doing x+(y+z).
454*
455 work( i ) = z( j ) / ( dlamc3( dsigj,
456 $ -poles( i+1,
457 $ 2 ) )-difr( i, 1 ) ) /
458 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
459 END IF
460 60 CONTINUE
461 DO 70 i = j + 1, k
462 IF( z( j ).EQ.zero ) THEN
463 work( i ) = zero
464 ELSE
465 work( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
466 $ 2 ) )-difl( i ) ) /
467 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
468 END IF
469 70 CONTINUE
470 CALL dgemv( 'T', k, nrhs, one, b, ldb, work, 1, zero,
471 $ bx( j, 1 ), ldbx )
472 80 CONTINUE
473 END IF
474*
475* Step (2R): if SQRE = 1, apply back the rotation that is
476* related to the right null space of the subproblem.
477*
478 IF( sqre.EQ.1 ) THEN
479 CALL dcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
480 CALL drot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c,
481 $ s )
482 END IF
483 IF( k.LT.max( m, n ) )
484 $ CALL dlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1,
485 $ 1 ),
486 $ ldbx )
487*
488* Step (3R): permute rows of B.
489*
490 CALL dcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
491 IF( sqre.EQ.1 ) THEN
492 CALL dcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
493 END IF
494 DO 90 i = 2, n
495 CALL dcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ),
496 $ ldb )
497 90 CONTINUE
498*
499* Step (4R): apply back the Givens rotations performed.
500*
501 DO 100 i = givptr, 1, -1
502 CALL drot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
503 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
504 $ -givnum( i, 1 ) )
505 100 CONTINUE
506 END IF
507*
508 RETURN
509*
510* End of DLALS0
511*
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:101
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:142
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: