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

◆ 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.

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

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.
!>
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 clals0.f.

269*
270* -- LAPACK computational routine --
271* -- LAPACK is a software package provided by Univ. of Tennessee, --
272* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
273*
274* .. Scalar Arguments ..
275 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
276 $ LDGNUM, NL, NR, NRHS, SQRE
277 REAL C, S
278* ..
279* .. Array Arguments ..
280 INTEGER GIVCOL( LDGCOL, * ), PERM( * )
281 REAL DIFL( * ), DIFR( LDGNUM, * ),
282 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
283 $ RWORK( * ), Z( * )
284 COMPLEX B( LDB, * ), BX( LDBX, * )
285* ..
286*
287* =====================================================================
288*
289* .. Parameters ..
290 REAL ONE, ZERO, NEGONE
291 parameter( one = 1.0e0, zero = 0.0e0, negone = -1.0e0 )
292* ..
293* .. Local Scalars ..
294 INTEGER I, J, JCOL, JROW, M, N, NLP1
295 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
296* ..
297* .. External Subroutines ..
298 EXTERNAL ccopy, clacpy, clascl, csrot, csscal,
299 $ 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 ),
364 $ ldbx )
365 20 CONTINUE
366*
367* Step (3L): apply the inverse of the left singular vector
368* matrix to BX.
369*
370 IF( k.EQ.1 ) THEN
371 CALL ccopy( nrhs, bx, ldbx, b, ldb )
372 IF( z( 1 ).LT.zero ) THEN
373 CALL csscal( nrhs, negone, b, ldb )
374 END IF
375 ELSE
376 DO 100 j = 1, k
377 diflj = difl( j )
378 dj = poles( j, 1 )
379 dsigj = -poles( j, 2 )
380 IF( j.LT.k ) THEN
381 difrj = -difr( j, 1 )
382 dsigjp = -poles( j+1, 2 )
383 END IF
384 IF( ( z( j ).EQ.zero ) .OR. ( poles( j, 2 ).EQ.zero ) )
385 $ THEN
386 rwork( j ) = zero
387 ELSE
388 rwork( j ) = -poles( j, 2 )*z( j ) / diflj /
389 $ ( poles( j, 2 )+dj )
390 END IF
391 DO 30 i = 1, j - 1
392 IF( ( z( i ).EQ.zero ) .OR.
393 $ ( poles( i, 2 ).EQ.zero ) ) THEN
394 rwork( i ) = zero
395 ELSE
396*
397* Use calls to the subroutine SLAMC3 to enforce the
398* parentheses (x+y)+z. The goal is to prevent
399* optimizing compilers from doing x+(y+z).
400*
401 rwork( i ) = poles( i, 2 )*z( i ) /
402 $ ( slamc3( 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 $ ( slamc3( poles( i, 2 ), dsigjp )+
413 $ difrj ) / ( poles( i, 2 )+dj )
414 END IF
415 40 CONTINUE
416 rwork( 1 ) = negone
417 temp = snrm2( k, rwork, 1 )
418*
419* Since B and BX are complex, the following call to SGEMV
420* is performed in two steps (real and imaginary parts).
421*
422* CALL SGEMV( '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 ) = real( bx( jrow, jcol ) )
430 50 CONTINUE
431 60 CONTINUE
432 CALL sgemv( '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 ) = aimag( bx( jrow, jcol ) )
439 70 CONTINUE
440 80 CONTINUE
441 CALL sgemv( '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 ) = cmplx( rwork( jcol+k ),
445 $ rwork( jcol+k+nrhs ) )
446 90 CONTINUE
447 CALL clascl( '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 clacpy( '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 ccopy( 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*
480* Use calls to the subroutine SLAMC3 to enforce the
481* parentheses (x+y)+z. The goal is to prevent optimizing
482* compilers from doing x+(y+z).
483*
484 rwork( i ) = z( j ) / ( slamc3( dsigj,
485 $ -poles( i+1,
486 $ 2 ) )-difr( i, 1 ) ) /
487 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
488 END IF
489 110 CONTINUE
490 DO 120 i = j + 1, k
491 IF( z( j ).EQ.zero ) THEN
492 rwork( i ) = zero
493 ELSE
494 rwork( i ) = z( j ) / ( slamc3( dsigj,
495 $ -poles( i,
496 $ 2 ) )-difl( i ) ) /
497 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
498 END IF
499 120 CONTINUE
500*
501* Since B and BX are complex, the following call to SGEMV
502* is performed in two steps (real and imaginary parts).
503*
504* CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
505* $ BX( J, 1 ), LDBX )
506*
507 i = k + nrhs*2
508 DO 140 jcol = 1, nrhs
509 DO 130 jrow = 1, k
510 i = i + 1
511 rwork( i ) = real( b( jrow, jcol ) )
512 130 CONTINUE
513 140 CONTINUE
514 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
515 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
516 i = k + nrhs*2
517 DO 160 jcol = 1, nrhs
518 DO 150 jrow = 1, k
519 i = i + 1
520 rwork( i ) = aimag( b( jrow, jcol ) )
521 150 CONTINUE
522 160 CONTINUE
523 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
524 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
525 DO 170 jcol = 1, nrhs
526 bx( j, jcol ) = cmplx( rwork( jcol+k ),
527 $ rwork( jcol+k+nrhs ) )
528 170 CONTINUE
529 180 CONTINUE
530 END IF
531*
532* Step (2R): if SQRE = 1, apply back the rotation that is
533* related to the right null space of the subproblem.
534*
535 IF( sqre.EQ.1 ) THEN
536 CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
537 CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c,
538 $ s )
539 END IF
540 IF( k.LT.max( m, n ) )
541 $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
542 $ bx( k+1, 1 ), ldbx )
543*
544* Step (3R): permute rows of B.
545*
546 CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
547 IF( sqre.EQ.1 ) THEN
548 CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
549 END IF
550 DO 190 i = 2, n
551 CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ),
552 $ ldb )
553 190 CONTINUE
554*
555* Step (4R): apply back the Givens rotations performed.
556*
557 DO 200 i = givptr, 1, -1
558 CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
559 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
560 $ -givnum( i, 1 ) )
561 200 CONTINUE
562 END IF
563*
564 RETURN
565*
566* End of CLALS0
567*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
slamc3
real function slamc3(a, b)
SLAMC3
Definition slamch.f:171
clascl
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:142
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
csrot
subroutine csrot(n, cx, incx, cy, incy, c, s)
CSROT
Definition csrot.f:98
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
Here is the call graph for this function:
Here is the caller graph for this function: