LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ zlals0()

 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.```
Date
December 2016
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.7.0) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * December 2016
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  n = nl + nr + 1
318 *
319  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
320  info = -1
321  ELSE IF( nl.LT.1 ) THEN
322  info = -2
323  ELSE IF( nr.LT.1 ) THEN
324  info = -3
325  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
326  info = -4
327  ELSE IF( nrhs.LT.1 ) THEN
328  info = -5
329  ELSE IF( ldb.LT.n ) THEN
330  info = -7
331  ELSE IF( ldbx.LT.n ) THEN
332  info = -9
333  ELSE IF( givptr.LT.0 ) THEN
334  info = -11
335  ELSE IF( ldgcol.LT.n ) THEN
336  info = -13
337  ELSE IF( ldgnum.LT.n ) THEN
338  info = -15
339  ELSE IF( k.LT.1 ) THEN
340  info = -20
341  END IF
342  IF( info.NE.0 ) THEN
343  CALL xerbla( 'ZLALS0', -info )
344  RETURN
345  END IF
346 *
347  m = n + sqre
348  nlp1 = nl + 1
349 *
350  IF( icompq.EQ.0 ) THEN
351 *
352 * Apply back orthogonal transformations from the left.
353 *
354 * Step (1L): apply back the Givens rotations performed.
355 *
356  DO 10 i = 1, givptr
357  CALL zdrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
358  \$ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
359  \$ givnum( i, 1 ) )
360  10 CONTINUE
361 *
362 * Step (2L): permute rows of B.
363 *
364  CALL zcopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
365  DO 20 i = 2, n
366  CALL zcopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
367  20 CONTINUE
368 *
369 * Step (3L): apply the inverse of the left singular vector
370 * matrix to BX.
371 *
372  IF( k.EQ.1 ) THEN
373  CALL zcopy( nrhs, bx, ldbx, b, ldb )
374  IF( z( 1 ).LT.zero ) THEN
375  CALL zdscal( nrhs, negone, b, ldb )
376  END IF
377  ELSE
378  DO 100 j = 1, k
379  diflj = difl( j )
380  dj = poles( j, 1 )
381  dsigj = -poles( j, 2 )
382  IF( j.LT.k ) THEN
383  difrj = -difr( j, 1 )
384  dsigjp = -poles( j+1, 2 )
385  END IF
386  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
387  \$ THEN
388  rwork( j ) = zero
389  ELSE
390  rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
391  \$ ( poles( j, 2 )+dj )
392  END IF
393  DO 30 i = 1, j - 1
394  IF( ( z( i ).EQ.zero ) .OR.
395  \$ ( poles( i, 2 ).EQ.zero ) ) THEN
396  rwork( i ) = zero
397  ELSE
398  rwork( 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  rwork( i ) = zero
407  ELSE
408  rwork( 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  rwork( 1 ) = negone
414  temp = dnrm2( k, rwork, 1 )
415 *
416 * Since B and BX are complex, the following call to DGEMV
417 * is performed in two steps (real and imaginary parts).
418 *
419 * CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
420 * \$ B( J, 1 ), LDB )
421 *
422  i = k + nrhs*2
423  DO 60 jcol = 1, nrhs
424  DO 50 jrow = 1, k
425  i = i + 1
426  rwork( i ) = dble( bx( jrow, jcol ) )
427  50 CONTINUE
428  60 CONTINUE
429  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
430  \$ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
431  i = k + nrhs*2
432  DO 80 jcol = 1, nrhs
433  DO 70 jrow = 1, k
434  i = i + 1
435  rwork( i ) = dimag( bx( jrow, jcol ) )
436  70 CONTINUE
437  80 CONTINUE
438  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
439  \$ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
440  DO 90 jcol = 1, nrhs
441  b( j, jcol ) = dcmplx( rwork( jcol+k ),
442  \$ rwork( jcol+k+nrhs ) )
443  90 CONTINUE
444  CALL zlascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
445  \$ ldb, info )
446  100 CONTINUE
447  END IF
448 *
449 * Move the deflated rows of BX to B also.
450 *
451  IF( k.LT.max( m, n ) )
452  \$ CALL zlacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
453  \$ b( k+1, 1 ), ldb )
454  ELSE
455 *
456 * Apply back the right orthogonal transformations.
457 *
458 * Step (1R): apply back the new right singular vector matrix
459 * to B.
460 *
461  IF( k.EQ.1 ) THEN
462  CALL zcopy( nrhs, b, ldb, bx, ldbx )
463  ELSE
464  DO 180 j = 1, k
465  dsigj = poles( j, 2 )
466  IF( z( j ).EQ.zero ) THEN
467  rwork( j ) = zero
468  ELSE
469  rwork( j ) = -z( j ) / difl( j ) /
470  \$ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
471  END IF
472  DO 110 i = 1, j - 1
473  IF( z( j ).EQ.zero ) THEN
474  rwork( i ) = zero
475  ELSE
476  rwork( i ) = z( j ) / ( dlamc3( dsigj, -poles( i+1,
477  \$ 2 ) )-difr( i, 1 ) ) /
478  \$ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
479  END IF
480  110 CONTINUE
481  DO 120 i = j + 1, k
482  IF( z( j ).EQ.zero ) THEN
483  rwork( i ) = zero
484  ELSE
485  rwork( i ) = z( j ) / ( dlamc3( dsigj, -poles( i,
486  \$ 2 ) )-difl( i ) ) /
487  \$ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
488  END IF
489  120 CONTINUE
490 *
491 * Since B and BX are complex, the following call to DGEMV
492 * is performed in two steps (real and imaginary parts).
493 *
494 * CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
495 * \$ BX( J, 1 ), LDBX )
496 *
497  i = k + nrhs*2
498  DO 140 jcol = 1, nrhs
499  DO 130 jrow = 1, k
500  i = i + 1
501  rwork( i ) = dble( b( jrow, jcol ) )
502  130 CONTINUE
503  140 CONTINUE
504  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
505  \$ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
506  i = k + nrhs*2
507  DO 160 jcol = 1, nrhs
508  DO 150 jrow = 1, k
509  i = i + 1
510  rwork( i ) = dimag( b( jrow, jcol ) )
511  150 CONTINUE
512  160 CONTINUE
513  CALL dgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
514  \$ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
515  DO 170 jcol = 1, nrhs
516  bx( j, jcol ) = dcmplx( rwork( jcol+k ),
517  \$ rwork( jcol+k+nrhs ) )
518  170 CONTINUE
519  180 CONTINUE
520  END IF
521 *
522 * Step (2R): if SQRE = 1, apply back the rotation that is
523 * related to the right null space of the subproblem.
524 *
525  IF( sqre.EQ.1 ) THEN
526  CALL zcopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
527  CALL zdrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
528  END IF
529  IF( k.LT.max( m, n ) )
530  \$ CALL zlacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb, bx( k+1, 1 ),
531  \$ ldbx )
532 *
533 * Step (3R): permute rows of B.
534 *
535  CALL zcopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
536  IF( sqre.EQ.1 ) THEN
537  CALL zcopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
538  END IF
539  DO 190 i = 2, n
540  CALL zcopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
541  190 CONTINUE
542 *
543 * Step (4R): apply back the Givens rotations performed.
544 *
545  DO 200 i = givptr, 1, -1
546  CALL zdrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
547  \$ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
548  \$ -givnum( i, 1 ) )
549  200 CONTINUE
550  END IF
551 *
552  RETURN
553 *
554 * End of ZLALS0
555 *
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
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:145
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 zdrot(N, CX, INCX, CY, INCY, C, S)
ZDROT
Definition: zdrot.f:100
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:76
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: