diff options
Diffstat (limited to 'SRC/zlals0.f')
| -rw-r--r-- | SRC/zlals0.f | 433 |
1 files changed, 433 insertions, 0 deletions
diff --git a/SRC/zlals0.f b/SRC/zlals0.f new file mode 100644 index 00000000..9d419612 --- /dev/null +++ b/SRC/zlals0.f @@ -0,0 +1,433 @@ + SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, + $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, + $ POLES, DIFL, DIFR, Z, K, C, S, RWORK, 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 + DOUBLE PRECISION C, S +* .. +* .. Array Arguments .. + INTEGER GIVCOL( LDGCOL, * ), PERM( * ) + DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ), + $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), + $ RWORK( * ), Z( * ) + COMPLEX*16 B( LDB, * ), BX( LDBX, * ) +* .. +* +* 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). +* +* 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) 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. +* +* LDB (input) INTEGER +* The leading dimension of B. LDB must be at least +* max(1,MAX( M, N ) ). +* +* BX (workspace) COMPLEX*16 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) DOUBLE PRECISION 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) 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. +* +* DIFL (input) 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. +* +* DIFR (input) 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. +* +* Z (input) DOUBLE PRECISION 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) 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. +* +* S (input) 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. +* +* RWORK (workspace) DOUBLE PRECISION array, dimension +* ( K*(1+NRHS) + 2*NRHS ) +* +* 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 .. + DOUBLE PRECISION ONE, ZERO, NEGONE + PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, JCOL, JROW, M, N, NLP1 + DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, XERBLA, ZCOPY, ZDROT, ZDSCAL, ZLACPY, + $ ZLASCL +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMC3, DNRM2 + EXTERNAL DLAMC3, DNRM2 +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, DCMPLX, DIMAG, 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( 'ZLALS0', -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 ZDROT( 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 ZCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) + DO 20 I = 2, N + CALL ZCOPY( 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 ZCOPY( NRHS, BX, LDBX, B, LDB ) + IF( Z( 1 ).LT.ZERO ) THEN + CALL ZDSCAL( NRHS, NEGONE, B, LDB ) + END IF + ELSE + DO 100 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 + RWORK( J ) = ZERO + ELSE + RWORK( 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 + RWORK( I ) = ZERO + ELSE + RWORK( I ) = POLES( I, 2 )*Z( I ) / + $ ( DLAMC3( 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 + RWORK( I ) = ZERO + ELSE + RWORK( I ) = POLES( I, 2 )*Z( I ) / + $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+ + $ DIFRJ ) / ( POLES( I, 2 )+DJ ) + END IF + 40 CONTINUE + RWORK( 1 ) = NEGONE + TEMP = DNRM2( K, RWORK, 1 ) +* +* Since B and BX are complex, the following call to DGEMV +* is performed in two steps (real and imaginary parts). +* +* CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO, +* $ B( J, 1 ), LDB ) +* + I = K + NRHS*2 + DO 60 JCOL = 1, NRHS + DO 50 JROW = 1, K + I = I + 1 + RWORK( I ) = DBLE( BX( JROW, JCOL ) ) + 50 CONTINUE + 60 CONTINUE + CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K, + $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 ) + I = K + NRHS*2 + DO 80 JCOL = 1, NRHS + DO 70 JROW = 1, K + I = I + 1 + RWORK( I ) = DIMAG( BX( JROW, JCOL ) ) + 70 CONTINUE + 80 CONTINUE + CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K, + $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 ) + DO 90 JCOL = 1, NRHS + B( J, JCOL ) = DCMPLX( RWORK( JCOL+K ), + $ RWORK( JCOL+K+NRHS ) ) + 90 CONTINUE + CALL ZLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ), + $ LDB, INFO ) + 100 CONTINUE + END IF +* +* Move the deflated rows of BX to B also. +* + IF( K.LT.MAX( M, N ) ) + $ CALL ZLACPY( '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 ZCOPY( NRHS, B, LDB, BX, LDBX ) + ELSE + DO 180 J = 1, K + DSIGJ = POLES( J, 2 ) + IF( Z( J ).EQ.ZERO ) THEN + RWORK( J ) = ZERO + ELSE + RWORK( J ) = -Z( J ) / DIFL( J ) / + $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 ) + END IF + DO 110 I = 1, J - 1 + IF( Z( J ).EQ.ZERO ) THEN + RWORK( I ) = ZERO + ELSE + RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, + $ 2 ) )-DIFR( I, 1 ) ) / + $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) + END IF + 110 CONTINUE + DO 120 I = J + 1, K + IF( Z( J ).EQ.ZERO ) THEN + RWORK( I ) = ZERO + ELSE + RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I, + $ 2 ) )-DIFL( I ) ) / + $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) + END IF + 120 CONTINUE +* +* Since B and BX are complex, the following call to DGEMV +* is performed in two steps (real and imaginary parts). +* +* CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO, +* $ BX( J, 1 ), LDBX ) +* + I = K + NRHS*2 + DO 140 JCOL = 1, NRHS + DO 130 JROW = 1, K + I = I + 1 + RWORK( I ) = DBLE( B( JROW, JCOL ) ) + 130 CONTINUE + 140 CONTINUE + CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K, + $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 ) + I = K + NRHS*2 + DO 160 JCOL = 1, NRHS + DO 150 JROW = 1, K + I = I + 1 + RWORK( I ) = DIMAG( B( JROW, JCOL ) ) + 150 CONTINUE + 160 CONTINUE + CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K, + $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 ) + DO 170 JCOL = 1, NRHS + BX( J, JCOL ) = DCMPLX( RWORK( JCOL+K ), + $ RWORK( JCOL+K+NRHS ) ) + 170 CONTINUE + 180 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 ZCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) + CALL ZDROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) + END IF + IF( K.LT.MAX( M, N ) ) + $ CALL ZLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ), + $ LDBX ) +* +* Step (3R): permute rows of B. +* + CALL ZCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) + IF( SQRE.EQ.1 ) THEN + CALL ZCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) + END IF + DO 190 I = 2, N + CALL ZCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) + 190 CONTINUE +* +* Step (4R): apply back the Givens rotations performed. +* + DO 200 I = GIVPTR, 1, -1 + CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, + $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), + $ -GIVNUM( I, 1 ) ) + 200 CONTINUE + END IF +* + RETURN +* +* End of ZLALS0 +* + END |