378 lines
13 KiB
FortranFixed
378 lines
13 KiB
FortranFixed
SUBROUTINE SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
|
|
$ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
|
|
$ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.1) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
|
* November 2006
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
|
|
$ LDGNUM, NL, NR, NRHS, SQRE
|
|
REAL C, S
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER GIVCOL( LDGCOL, * ), PERM( * )
|
|
REAL B( LDB, * ), BX( LDBX, * ), DIFL( * ),
|
|
$ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ),
|
|
$ POLES( LDGNUM, * ), WORK( * ), Z( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* SLALS0 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).
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* ICOMPQ (input) INTEGER
|
|
* Specifies whether singular vectors are to be computed in
|
|
* factored form:
|
|
* = 0: Left singular vector matrix.
|
|
* = 1: Right singular vector matrix.
|
|
*
|
|
* NL (input) INTEGER
|
|
* The row dimension of the upper block. NL >= 1.
|
|
*
|
|
* NR (input) INTEGER
|
|
* The row dimension of the lower block. NR >= 1.
|
|
*
|
|
* SQRE (input) 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.
|
|
*
|
|
* NRHS (input) INTEGER
|
|
* The number of columns of B and BX. NRHS must be at least 1.
|
|
*
|
|
* B (input/output) REAL 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.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of B. LDB must be at least
|
|
* max(1,MAX( M, N ) ).
|
|
*
|
|
* BX (workspace) REAL array, dimension ( LDBX, NRHS )
|
|
*
|
|
* LDBX (input) INTEGER
|
|
* The leading dimension of BX.
|
|
*
|
|
* PERM (input) INTEGER array, dimension ( N )
|
|
* The permutations (from deflation and sorting) applied
|
|
* to the two blocks.
|
|
*
|
|
* GIVPTR (input) INTEGER
|
|
* The number of Givens rotations which took place in this
|
|
* subproblem.
|
|
*
|
|
* GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
|
|
* Each pair of numbers indicates a pair of rows/columns
|
|
* involved in a Givens rotation.
|
|
*
|
|
* LDGCOL (input) INTEGER
|
|
* The leading dimension of GIVCOL, must be at least N.
|
|
*
|
|
* GIVNUM (input) REAL array, dimension ( LDGNUM, 2 )
|
|
* Each number indicates the C or S value used in the
|
|
* corresponding Givens rotation.
|
|
*
|
|
* LDGNUM (input) INTEGER
|
|
* The leading dimension of arrays DIFR, POLES and
|
|
* GIVNUM, must be at least K.
|
|
*
|
|
* POLES (input) 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.
|
|
*
|
|
* DIFL (input) 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.
|
|
*
|
|
* DIFR (input) 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.
|
|
*
|
|
* Z (input) REAL array, dimension ( K )
|
|
* Contain the components of the deflation-adjusted updating row
|
|
* vector.
|
|
*
|
|
* K (input) INTEGER
|
|
* Contains the dimension of the non-deflated matrix,
|
|
* This is the order of the related secular equation. 1 <= K <=N.
|
|
*
|
|
* C (input) REAL
|
|
* C contains garbage if SQRE =0 and the C-value of a Givens
|
|
* rotation related to the right null space if SQRE = 1.
|
|
*
|
|
* S (input) REAL
|
|
* S contains garbage if SQRE =0 and the S-value of a Givens
|
|
* rotation related to the right null space if SQRE = 1.
|
|
*
|
|
* WORK (workspace) REAL array, dimension ( K )
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit.
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value.
|
|
*
|
|
* Further Details
|
|
* ===============
|
|
*
|
|
* Based on contributions by
|
|
* Ming Gu and Ren-Cang Li, Computer Science Division, University of
|
|
* California at Berkeley, USA
|
|
* Osni Marques, LBNL/NERSC, USA
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
REAL ONE, ZERO, NEGONE
|
|
PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER I, J, M, N, NLP1
|
|
REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL SCOPY, SGEMV, SLACPY, SLASCL, SROT, SSCAL,
|
|
$ XERBLA
|
|
* ..
|
|
* .. External Functions ..
|
|
REAL SLAMC3, SNRM2
|
|
EXTERNAL SLAMC3, SNRM2
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
*
|
|
IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( NL.LT.1 ) THEN
|
|
INFO = -2
|
|
ELSE IF( NR.LT.1 ) THEN
|
|
INFO = -3
|
|
ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
*
|
|
N = NL + NR + 1
|
|
*
|
|
IF( NRHS.LT.1 ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDB.LT.N ) THEN
|
|
INFO = -7
|
|
ELSE IF( LDBX.LT.N ) THEN
|
|
INFO = -9
|
|
ELSE IF( GIVPTR.LT.0 ) THEN
|
|
INFO = -11
|
|
ELSE IF( LDGCOL.LT.N ) THEN
|
|
INFO = -13
|
|
ELSE IF( LDGNUM.LT.N ) THEN
|
|
INFO = -15
|
|
ELSE IF( K.LT.1 ) THEN
|
|
INFO = -20
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'SLALS0', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
M = N + SQRE
|
|
NLP1 = NL + 1
|
|
*
|
|
IF( ICOMPQ.EQ.0 ) THEN
|
|
*
|
|
* Apply back orthogonal transformations from the left.
|
|
*
|
|
* Step (1L): apply back the Givens rotations performed.
|
|
*
|
|
DO 10 I = 1, GIVPTR
|
|
CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
|
|
$ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
|
|
$ GIVNUM( I, 1 ) )
|
|
10 CONTINUE
|
|
*
|
|
* Step (2L): permute rows of B.
|
|
*
|
|
CALL SCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
|
|
DO 20 I = 2, N
|
|
CALL SCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
|
|
20 CONTINUE
|
|
*
|
|
* Step (3L): apply the inverse of the left singular vector
|
|
* matrix to BX.
|
|
*
|
|
IF( K.EQ.1 ) THEN
|
|
CALL SCOPY( NRHS, BX, LDBX, B, LDB )
|
|
IF( Z( 1 ).LT.ZERO ) THEN
|
|
CALL SSCAL( NRHS, NEGONE, B, LDB )
|
|
END IF
|
|
ELSE
|
|
DO 50 J = 1, K
|
|
DIFLJ = DIFL( J )
|
|
DJ = POLES( J, 1 )
|
|
DSIGJ = -POLES( J, 2 )
|
|
IF( J.LT.K ) THEN
|
|
DIFRJ = -DIFR( J, 1 )
|
|
DSIGJP = -POLES( J+1, 2 )
|
|
END IF
|
|
IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
|
|
$ THEN
|
|
WORK( J ) = ZERO
|
|
ELSE
|
|
WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
|
|
$ ( POLES( J, 2 )+DJ )
|
|
END IF
|
|
DO 30 I = 1, J - 1
|
|
IF( ( Z( I ).EQ.ZERO ) .OR.
|
|
$ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
|
|
WORK( I ) = ZERO
|
|
ELSE
|
|
WORK( I ) = POLES( I, 2 )*Z( I ) /
|
|
$ ( SLAMC3( POLES( I, 2 ), DSIGJ )-
|
|
$ DIFLJ ) / ( POLES( I, 2 )+DJ )
|
|
END IF
|
|
30 CONTINUE
|
|
DO 40 I = J + 1, K
|
|
IF( ( Z( I ).EQ.ZERO ) .OR.
|
|
$ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
|
|
WORK( I ) = ZERO
|
|
ELSE
|
|
WORK( I ) = POLES( I, 2 )*Z( I ) /
|
|
$ ( SLAMC3( POLES( I, 2 ), DSIGJP )+
|
|
$ DIFRJ ) / ( POLES( I, 2 )+DJ )
|
|
END IF
|
|
40 CONTINUE
|
|
WORK( 1 ) = NEGONE
|
|
TEMP = SNRM2( K, WORK, 1 )
|
|
CALL SGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
|
|
$ B( J, 1 ), LDB )
|
|
CALL SLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
|
|
$ LDB, INFO )
|
|
50 CONTINUE
|
|
END IF
|
|
*
|
|
* Move the deflated rows of BX to B also.
|
|
*
|
|
IF( K.LT.MAX( M, N ) )
|
|
$ CALL SLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
|
|
$ B( K+1, 1 ), LDB )
|
|
ELSE
|
|
*
|
|
* Apply back the right orthogonal transformations.
|
|
*
|
|
* Step (1R): apply back the new right singular vector matrix
|
|
* to B.
|
|
*
|
|
IF( K.EQ.1 ) THEN
|
|
CALL SCOPY( NRHS, B, LDB, BX, LDBX )
|
|
ELSE
|
|
DO 80 J = 1, K
|
|
DSIGJ = POLES( J, 2 )
|
|
IF( Z( J ).EQ.ZERO ) THEN
|
|
WORK( J ) = ZERO
|
|
ELSE
|
|
WORK( J ) = -Z( J ) / DIFL( J ) /
|
|
$ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
|
|
END IF
|
|
DO 60 I = 1, J - 1
|
|
IF( Z( J ).EQ.ZERO ) THEN
|
|
WORK( I ) = ZERO
|
|
ELSE
|
|
WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1,
|
|
$ 2 ) )-DIFR( I, 1 ) ) /
|
|
$ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
|
|
END IF
|
|
60 CONTINUE
|
|
DO 70 I = J + 1, K
|
|
IF( Z( J ).EQ.ZERO ) THEN
|
|
WORK( I ) = ZERO
|
|
ELSE
|
|
WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I,
|
|
$ 2 ) )-DIFL( I ) ) /
|
|
$ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
|
|
END IF
|
|
70 CONTINUE
|
|
CALL SGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
|
|
$ BX( J, 1 ), LDBX )
|
|
80 CONTINUE
|
|
END IF
|
|
*
|
|
* Step (2R): if SQRE = 1, apply back the rotation that is
|
|
* related to the right null space of the subproblem.
|
|
*
|
|
IF( SQRE.EQ.1 ) THEN
|
|
CALL SCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
|
|
CALL SROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
|
|
END IF
|
|
IF( K.LT.MAX( M, N ) )
|
|
$ CALL SLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
|
|
$ LDBX )
|
|
*
|
|
* Step (3R): permute rows of B.
|
|
*
|
|
CALL SCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
|
|
IF( SQRE.EQ.1 ) THEN
|
|
CALL SCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
|
|
END IF
|
|
DO 90 I = 2, N
|
|
CALL SCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
|
|
90 CONTINUE
|
|
*
|
|
* Step (4R): apply back the Givens rotations performed.
|
|
*
|
|
DO 100 I = GIVPTR, 1, -1
|
|
CALL SROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
|
|
$ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
|
|
$ -GIVNUM( I, 1 ) )
|
|
100 CONTINUE
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of SLALS0
|
|
*
|
|
END
|