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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine zlals0 ( integer  ICOMPQ,
integer  NL,
integer  NR,
integer  SQRE,
integer  NRHS,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, 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( * )  RWORK,
integer  INFO 
)

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

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

Purpose:
 ZLALS0 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*16 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*16 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]RWORK
          RWORK is DOUBLE PRECISION 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 zlals0.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  DOUBLE PRECISION c, s
282 * ..
283 * .. Array Arguments ..
284  INTEGER givcol( ldgcol, * ), perm( * )
285  DOUBLE PRECISION difl( * ), difr( ldgnum, * ),
286  $ givnum( ldgnum, * ), poles( ldgnum, * ),
287  $ rwork( * ), z( * )
288  COMPLEX*16 b( ldb, * ), bx( ldbx, * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  DOUBLE PRECISION one, zero, negone
295  parameter( one = 1.0d0, zero = 0.0d0, negone = -1.0d0 )
296 * ..
297 * .. Local Scalars ..
298  INTEGER i, j, jcol, jrow, m, n, nlp1
299  DOUBLE PRECISION diflj, difrj, dj, dsigj, dsigjp, temp
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL dgemv, xerbla, zcopy, zdrot, zdscal, zlacpy,
303  $ zlascl
304 * ..
305 * .. External Functions ..
306  DOUBLE PRECISION dlamc3, dnrm2
307  EXTERNAL dlamc3, dnrm2
308 * ..
309 * .. Intrinsic Functions ..
310  INTRINSIC dble, dcmplx, dimag, max
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( 'ZLALS0', -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 zdrot( 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 zcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
368  DO 20 i = 2, n
369  CALL zcopy( 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 zcopy( nrhs, bx, ldbx, b, ldb )
377  IF( z( 1 ).LT.zero ) THEN
378  CALL zdscal( 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  $ ( dlamc3( 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  $ ( dlamc3( poles( i, 2 ), dsigjp )+
413  $ difrj ) / ( poles( i, 2 )+dj )
414  END IF
415  40 CONTINUE
416  rwork( 1 ) = negone
417  temp = dnrm2( k, rwork, 1 )
418 *
419 * Since B and BX are complex, the following call to DGEMV
420 * is performed in two steps (real and imaginary parts).
421 *
422 * CALL DGEMV( '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 ) = dble( bx( jrow, jcol ) )
430  50 CONTINUE
431  60 CONTINUE
432  CALL dgemv( '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 ) = dimag( bx( jrow, jcol ) )
439  70 CONTINUE
440  80 CONTINUE
441  CALL dgemv( '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 ) = dcmplx( rwork( jcol+k ),
445  $ rwork( jcol+k+nrhs ) )
446  90 CONTINUE
447  CALL zlascl( '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 zlacpy( '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 zcopy( 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 ) / ( dlamc3( 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 ) / ( dlamc3( 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 DGEMV
495 * is performed in two steps (real and imaginary parts).
496 *
497 * CALL DGEMV( '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 ) = dble( b( jrow, jcol ) )
505  130 CONTINUE
506  140 CONTINUE
507  CALL dgemv( '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 ) = dimag( b( jrow, jcol ) )
514  150 CONTINUE
515  160 CONTINUE
516  CALL dgemv( '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 ) = dcmplx( 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 zcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
530  CALL zdrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
531  END IF
532  IF( k.LT.max( m, n ) )
533  $ CALL zlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
534  $ ldbx )
535 *
536 * Step (3R): permute rows of B.
537 *
538  CALL zcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
539  IF( sqre.EQ.1 ) THEN
540  CALL zcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
541  END IF
542  DO 190 i = 2, n
543  CALL zcopy( 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 zdrot( 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 ZLALS0
558 *
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:54
subroutine zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:141
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:52
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:173
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158

Here is the call graph for this function:

Here is the caller graph for this function: