LAPACK 3.12.0
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 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*
396* Use calls to the subroutine SLAMC3 to enforce the
397* parentheses (x+y)+z. The goal is to prevent
398* optimizing compilers from doing x+(y+z).
399*
400 rwork( i ) = poles( i, 2 )*z( i ) /
401 $ ( slamc3( poles( i, 2 ), dsigj )-
402 $ diflj ) / ( poles( i, 2 )+dj )
403 END IF
404 30 CONTINUE
405 DO 40 i = j + 1, k
406 IF( ( z( i ).EQ.zero ) .OR.
407 $ ( poles( i, 2 ).EQ.zero ) ) THEN
408 rwork( i ) = zero
409 ELSE
410 rwork( i ) = poles( i, 2 )*z( i ) /
411 $ ( slamc3( poles( i, 2 ), dsigjp )+
412 $ difrj ) / ( poles( i, 2 )+dj )
413 END IF
414 40 CONTINUE
415 rwork( 1 ) = negone
416 temp = snrm2( k, rwork, 1 )
417*
418* Since B and BX are complex, the following call to SGEMV
419* is performed in two steps (real and imaginary parts).
420*
421* CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
422* $ B( J, 1 ), LDB )
423*
424 i = k + nrhs*2
425 DO 60 jcol = 1, nrhs
426 DO 50 jrow = 1, k
427 i = i + 1
428 rwork( i ) = real( bx( jrow, jcol ) )
429 50 CONTINUE
430 60 CONTINUE
431 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
432 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
433 i = k + nrhs*2
434 DO 80 jcol = 1, nrhs
435 DO 70 jrow = 1, k
436 i = i + 1
437 rwork( i ) = aimag( bx( jrow, jcol ) )
438 70 CONTINUE
439 80 CONTINUE
440 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
441 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
442 DO 90 jcol = 1, nrhs
443 b( j, jcol ) = cmplx( rwork( jcol+k ),
444 $ rwork( jcol+k+nrhs ) )
445 90 CONTINUE
446 CALL clascl( 'G', 0, 0, temp, one, 1, nrhs, b( j, 1 ),
447 $ ldb, info )
448 100 CONTINUE
449 END IF
450*
451* Move the deflated rows of BX to B also.
452*
453 IF( k.LT.max( m, n ) )
454 $ CALL clacpy( 'A', n-k, nrhs, bx( k+1, 1 ), ldbx,
455 $ b( k+1, 1 ), ldb )
456 ELSE
457*
458* Apply back the right orthogonal transformations.
459*
460* Step (1R): apply back the new right singular vector matrix
461* to B.
462*
463 IF( k.EQ.1 ) THEN
464 CALL ccopy( nrhs, b, ldb, bx, ldbx )
465 ELSE
466 DO 180 j = 1, k
467 dsigj = poles( j, 2 )
468 IF( z( j ).EQ.zero ) THEN
469 rwork( j ) = zero
470 ELSE
471 rwork( j ) = -z( j ) / difl( j ) /
472 $ ( dsigj+poles( j, 1 ) ) / difr( j, 2 )
473 END IF
474 DO 110 i = 1, j - 1
475 IF( z( j ).EQ.zero ) THEN
476 rwork( i ) = zero
477 ELSE
478*
479* Use calls to the subroutine SLAMC3 to enforce the
480* parentheses (x+y)+z. The goal is to prevent optimizing
481* compilers from doing x+(y+z).
482*
483 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i+1,
484 $ 2 ) )-difr( i, 1 ) ) /
485 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
486 END IF
487 110 CONTINUE
488 DO 120 i = j + 1, k
489 IF( z( j ).EQ.zero ) THEN
490 rwork( i ) = zero
491 ELSE
492 rwork( i ) = z( j ) / ( slamc3( dsigj, -poles( i,
493 $ 2 ) )-difl( i ) ) /
494 $ ( dsigj+poles( i, 1 ) ) / difr( i, 2 )
495 END IF
496 120 CONTINUE
497*
498* Since B and BX are complex, the following call to SGEMV
499* is performed in two steps (real and imaginary parts).
500*
501* CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
502* $ BX( J, 1 ), LDBX )
503*
504 i = k + nrhs*2
505 DO 140 jcol = 1, nrhs
506 DO 130 jrow = 1, k
507 i = i + 1
508 rwork( i ) = real( b( jrow, jcol ) )
509 130 CONTINUE
510 140 CONTINUE
511 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
512 $ rwork( 1 ), 1, zero, rwork( 1+k ), 1 )
513 i = k + nrhs*2
514 DO 160 jcol = 1, nrhs
515 DO 150 jrow = 1, k
516 i = i + 1
517 rwork( i ) = aimag( b( jrow, jcol ) )
518 150 CONTINUE
519 160 CONTINUE
520 CALL sgemv( 'T', k, nrhs, one, rwork( 1+k+nrhs*2 ), k,
521 $ rwork( 1 ), 1, zero, rwork( 1+k+nrhs ), 1 )
522 DO 170 jcol = 1, nrhs
523 bx( j, jcol ) = cmplx( rwork( jcol+k ),
524 $ rwork( jcol+k+nrhs ) )
525 170 CONTINUE
526 180 CONTINUE
527 END IF
528*
529* Step (2R): if SQRE = 1, apply back the rotation that is
530* related to the right null space of the subproblem.
531*
532 IF( sqre.EQ.1 ) THEN
533 CALL ccopy( nrhs, b( m, 1 ), ldb, bx( m, 1 ), ldbx )
534 CALL csrot( nrhs, bx( 1, 1 ), ldbx, bx( m, 1 ), ldbx, c, s )
535 END IF
536 IF( k.LT.max( m, n ) )
537 $ CALL clacpy( 'A', n-k, nrhs, b( k+1, 1 ), ldb,
538 $ bx( k+1, 1 ), ldbx )
539*
540* Step (3R): permute rows of B.
541*
542 CALL ccopy( nrhs, bx( 1, 1 ), ldbx, b( nlp1, 1 ), ldb )
543 IF( sqre.EQ.1 ) THEN
544 CALL ccopy( nrhs, bx( m, 1 ), ldbx, b( m, 1 ), ldb )
545 END IF
546 DO 190 i = 2, n
547 CALL ccopy( nrhs, bx( i, 1 ), ldbx, b( perm( i ), 1 ), ldb )
548 190 CONTINUE
549*
550* Step (4R): apply back the Givens rotations performed.
551*
552 DO 200 i = givptr, 1, -1
553 CALL csrot( nrhs, b( givcol( i, 2 ), 1 ), ldb,
554 $ b( givcol( i, 1 ), 1 ), ldb, givnum( i, 2 ),
555 $ -givnum( i, 1 ) )
556 200 CONTINUE
557 END IF
558*
559 RETURN
560*
561* End of CLALS0
562*
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:103
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:143
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: