 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ clals0()

 subroutine clals0 ( integer ICOMPQ, integer NL, integer NR, integer SQRE, integer NRHS, complex, dimension( ldb, * ) B, integer LDB, complex, dimension( ldbx, * ) BX, integer LDBX, integer, dimension( * ) PERM, integer GIVPTR, integer, dimension( ldgcol, * ) GIVCOL, integer LDGCOL, real, dimension( ldgnum, * ) GIVNUM, integer LDGNUM, real, dimension( ldgnum, * ) POLES, real, dimension( * ) DIFL, real, dimension( ldgnum, * ) DIFR, real, dimension( * ) Z, integer K, real C, real S, real, dimension( * ) RWORK, integer INFO )

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

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

Definition at line 267 of file clals0.f.

270 *
271 * -- LAPACK computational routine --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 *
275 * .. Scalar Arguments ..
276  INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
277  \$ LDGNUM, NL, NR, NRHS, SQRE
278  REAL C, S
279 * ..
280 * .. Array Arguments ..
281  INTEGER GIVCOL( LDGCOL, * ), PERM( * )
282  REAL DIFL( * ), DIFR( LDGNUM, * ),
283  \$ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
284  \$ RWORK( * ), Z( * )
285  COMPLEX B( LDB, * ), BX( LDBX, * )
286 * ..
287 *
288 * =====================================================================
289 *
290 * .. Parameters ..
291  REAL ONE, ZERO, NEGONE
292  parameter( one = 1.0e0, zero = 0.0e0, negone = -1.0e0 )
293 * ..
294 * .. Local Scalars ..
295  INTEGER I, J, JCOL, JROW, M, N, NLP1
296  REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
297 * ..
298 * .. External Subroutines ..
299  EXTERNAL ccopy, clacpy, clascl, csrot, csscal, sgemv,
300  \$ xerbla
301 * ..
302 * .. External Functions ..
303  REAL SLAMC3, SNRM2
304  EXTERNAL slamc3, snrm2
305 * ..
306 * .. Intrinsic Functions ..
307  INTRINSIC aimag, cmplx, max, real
308 * ..
309 * .. Executable Statements ..
310 *
311 * Test the input parameters.
312 *
313  info = 0
314  n = nl + nr + 1
315 *
316  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
317  info = -1
318  ELSE IF( nl.LT.1 ) THEN
319  info = -2
320  ELSE IF( nr.LT.1 ) THEN
321  info = -3
322  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
323  info = -4
324  ELSE IF( nrhs.LT.1 ) THEN
325  info = -5
326  ELSE IF( ldb.LT.n ) THEN
327  info = -7
328  ELSE IF( ldbx.LT.n ) THEN
329  info = -9
330  ELSE IF( givptr.LT.0 ) THEN
331  info = -11
332  ELSE IF( ldgcol.LT.n ) THEN
333  info = -13
334  ELSE IF( ldgnum.LT.n ) THEN
335  info = -15
336  ELSE IF( k.LT.1 ) THEN
337  info = -20
338  END IF
339  IF( info.NE.0 ) THEN
340  CALL xerbla( 'CLALS0', -info )
341  RETURN
342  END IF
343 *
344  m = n + sqre
345  nlp1 = nl + 1
346 *
347  IF( icompq.EQ.0 ) THEN
348 *
349 * Apply back orthogonal transformations from the left.
350 *
351 * Step (1L): apply back the Givens rotations performed.
352 *
353  DO 10 i = 1, givptr
354  CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
355  \$ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
356  \$ givnum( i, 1 ) )
357  10 CONTINUE
358 *
359 * Step (2L): permute rows of B.
360 *
361  CALL ccopy( nrhs, b( nlp1, 1 ), ldb, bx( 1, 1 ), ldbx )
362  DO 20 i = 2, n
363  CALL ccopy( nrhs, b( perm( i ), 1 ), ldb, bx( i, 1 ), ldbx )
364  20 CONTINUE
365 *
366 * Step (3L): apply the inverse of the left singular vector
367 * matrix to BX.
368 *
369  IF( k.EQ.1 ) THEN
370  CALL ccopy( nrhs, bx, ldbx, b, ldb )
371  IF( z( 1 ).LT.zero ) THEN
372  CALL csscal( nrhs, negone, b, ldb )
373  END IF
374  ELSE
375  DO 100 j = 1, k
376  diflj = difl( j )
377  dj = poles( j, 1 )
378  dsigj = -poles( j, 2 )
379  IF( j.LT.k ) THEN
380  difrj = -difr( j, 1 )
381  dsigjp = -poles( j+1, 2 )
382  END IF
383  IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
384  \$ THEN
385  rwork( j ) = zero
386  ELSE
387  rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
388  \$ ( poles( j, 2 )+dj )
389  END IF
390  DO 30 i = 1, j - 1
391  IF( ( z( i ).EQ.zero ) .OR.
392  \$ ( poles( i, 2 ).EQ.zero ) ) THEN
393  rwork( i ) = zero
394  ELSE
395  rwork( i ) = poles( i, 2 )*z( i ) /
396  \$ ( slamc3( poles( i, 2 ), dsigj )-
397  \$ diflj ) / ( poles( i, 2 )+dj )
398  END IF
399  30 CONTINUE
400  DO 40 i = j + 1, k
401  IF( ( z( i ).EQ.zero ) .OR.
402  \$ ( poles( i, 2 ).EQ.zero ) ) THEN
403  rwork( i ) = zero
404  ELSE
405  rwork( i ) = poles( i, 2 )*z( i ) /
406  \$ ( slamc3( poles( i, 2 ), dsigjp )+
407  \$ difrj ) / ( poles( i, 2 )+dj )
408  END IF
409  40 CONTINUE
410  rwork( 1 ) = negone
411  temp = snrm2( k, rwork, 1 )
412 *
413 * Since B and BX are complex, the following call to SGEMV
414 * is performed in two steps (real and imaginary parts).
415 *
416 * CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
417 * \$ B( J, 1 ), LDB )
418 *
419  i = k + nrhs*2
420  DO 60 jcol = 1, nrhs
421  DO 50 jrow = 1, k
422  i = i + 1
423  rwork( i ) = real( bx( jrow, jcol ) )
424  50 CONTINUE
425  60 CONTINUE
426  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
427  \$ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
428  i = k + nrhs*2
429  DO 80 jcol = 1, nrhs
430  DO 70 jrow = 1, k
431  i = i + 1
432  rwork( i ) = aimag( bx( jrow, jcol ) )
433  70 CONTINUE
434  80 CONTINUE
435  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
436  \$ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
437  DO 90 jcol = 1, nrhs
438  b( j, jcol ) = cmplx( rwork( jcol+k ),
439  \$ rwork( jcol+k+nrhs ) )
440  90 CONTINUE
441  CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
442  \$ ldb, info )
443  100 CONTINUE
444  END IF
445 *
446 * Move the deflated rows of BX to B also.
447 *
448  IF( k.LT.max( m, n ) )
449  \$ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
450  \$ b( k+1, 1 ), ldb )
451  ELSE
452 *
453 * Apply back the right orthogonal transformations.
454 *
455 * Step (1R): apply back the new right singular vector matrix
456 * to B.
457 *
458  IF( k.EQ.1 ) THEN
459  CALL ccopy( nrhs, b, ldb, bx, ldbx )
460  ELSE
461  DO 180 j = 1, k
462  dsigj = poles( j, 2 )
463  IF( z( j ).EQ.zero ) THEN
464  rwork( j ) = zero
465  ELSE
466  rwork( j ) = -z( j ) / difl( j ) /
467  \$ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
468  END IF
469  DO 110 i = 1, j - 1
470  IF( z( j ).EQ.zero ) THEN
471  rwork( i ) = zero
472  ELSE
473  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
474  \$ 2 ) )-difr( i, 1 ) ) /
475  \$ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
476  END IF
477  110 CONTINUE
478  DO 120 i = j + 1, k
479  IF( z( j ).EQ.zero ) THEN
480  rwork( i ) = zero
481  ELSE
482  rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
483  \$ 2 ) )-difl( i ) ) /
484  \$ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
485  END IF
486  120 CONTINUE
487 *
488 * Since B and BX are complex, the following call to SGEMV
489 * is performed in two steps (real and imaginary parts).
490 *
491 * CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
492 * \$ BX( J, 1 ), LDBX )
493 *
494  i = k + nrhs*2
495  DO 140 jcol = 1, nrhs
496  DO 130 jrow = 1, k
497  i = i + 1
498  rwork( i ) = real( b( jrow, jcol ) )
499  130 CONTINUE
500  140 CONTINUE
501  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
502  \$ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
503  i = k + nrhs*2
504  DO 160 jcol = 1, nrhs
505  DO 150 jrow = 1, k
506  i = i + 1
507  rwork( i ) = aimag( b( jrow, jcol ) )
508  150 CONTINUE
509  160 CONTINUE
510  CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
511  \$ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
512  DO 170 jcol = 1, nrhs
513  bx( j, jcol ) = cmplx( rwork( jcol+k ),
514  \$ rwork( jcol+k+nrhs ) )
515  170 CONTINUE
516  180 CONTINUE
517  END IF
518 *
519 * Step (2R): if SQRE = 1, apply back the rotation that is
520 * related to the right null space of the subproblem.
521 *
522  IF( sqre.EQ.1 ) THEN
523  CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
524  CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
525  END IF
526  IF( k.LT.max( m, n ) )
527  \$ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
528  \$ bx( k+1, 1 ), ldbx )
529 *
530 * Step (3R): permute rows of B.
531 *
532  CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
533  IF( sqre.EQ.1 ) THEN
534  CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
535  END IF
536  DO 190 i = 2, n
537  CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
538  190 CONTINUE
539 *
540 * Step (4R): apply back the Givens rotations performed.
541 *
542  DO 200 i = givptr, 1, -1
543  CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
544  \$ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
545  \$ -givnum( i, 1 ) )
546  200 CONTINUE
547  END IF
548 *
549  RETURN
550 *
551 * End of CLALS0
552 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine csscal(N, SA, CX, INCX)
CSSCAL
Definition: csscal.f:78
subroutine csrot(N, CX, INCX, CY, INCY, C, S)
CSROT
Definition: csrot.f:98
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
real function slamc3(A, B)
SLAMC3
Definition: slamch.f:169
Here is the call graph for this function:
Here is the caller graph for this function: