SUBROUTINE CLALS0( 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 (instrumented to count ops, version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * December 22, 1999 * * .. 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 DIFL( * ), DIFR( LDGNUM, * ), $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), $ RWORK( * ), Z( * ) COMPLEX B( LDB, * ), BX( LDBX, * ) * .. * .. Common block to return operation count .. COMMON / LATIME / OPS, ITCNT * .. * .. Scalars in Common .. REAL ITCNT, OPS * .. * * 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). * * 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 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 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. * * RWORK (workspace) REAL 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. * * ===================================================================== * * .. Parameters .. REAL ONE, ZERO, NEGONE PARAMETER ( ONE = 1.0E0, ZERO = 0.0E0, NEGONE = -1.0E0 ) * .. * .. Local Scalars .. INTEGER I, J, JCOL, JROW, M, N, NLP1 REAL DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP * .. * .. External Subroutines .. EXTERNAL CSROT, CSSCAL, CCOPY, $ CLACPY, CLASCL, SGEMV, XERBLA * .. * .. External Functions .. REAL SLAMC3, SNRM2, SOPBL2 EXTERNAL SLAMC3, SNRM2, SOPBL2 * .. * .. Intrinsic Functions .. INTRINSIC CMPLX, REAL, AIMAG * .. * .. 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( 'CLALS0', -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. * OPS = OPS + REAL( 12*NRHS*GIVPTR ) DO 10 I = 1, GIVPTR CALL CSROT( 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 CCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) DO 20 I = 2, N CALL CCOPY( 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 CCOPY( NRHS, BX, LDBX, B, LDB ) IF( Z( 1 ).LT.ZERO ) THEN OPS = OPS + REAL( 6*NRHS ) CALL CSSCAL( 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 OPS = OPS + REAL( 4 ) 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 OPS = OPS + REAL( 6 ) RWORK( 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 RWORK( I ) = ZERO ELSE OPS = OPS + REAL( 6 ) RWORK( I ) = POLES( I, 2 )*Z( I ) / $ ( SLAMC3( POLES( I, 2 ), DSIGJP )+ $ DIFRJ ) / ( POLES( I, 2 )+DJ ) END IF 40 CONTINUE OPS = OPS + REAL( 2*K ) RWORK( 1 ) = NEGONE TEMP = SNRM2( K, RWORK, 1 ) * * Since B and BX are complex, the following call to SGEMV * is performed in two steps (real and imaginary parts). * * CALL SGEMV( '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 ) = REAL( BX( JROW, JCOL ) ) 50 CONTINUE 60 CONTINUE OPS = OPS + SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) CALL SGEMV( '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 ) = AIMAG( BX( JROW, JCOL ) ) 70 CONTINUE 80 CONTINUE OPS = OPS + SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) CALL SGEMV( '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 ) = CMPLX( RWORK( JCOL+K ), $ RWORK( JCOL+K+NRHS ) ) 90 CONTINUE OPS = OPS + REAL( 6*NRHS ) CALL CLASCL( '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 CLACPY( '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 CCOPY( 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 OPS = OPS + REAL( 4 ) 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 OPS = OPS + REAL( 6 ) RWORK( I ) = Z( J ) / ( SLAMC3( 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 OPS = OPS + REAL( 6 ) RWORK( I ) = Z( J ) / ( SLAMC3( 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 SGEMV * is performed in two steps (real and imaginary parts). * * CALL SGEMV( '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 ) = REAL( B( JROW, JCOL ) ) 130 CONTINUE 140 CONTINUE OPS = OPS + SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) CALL SGEMV( '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 ) = AIMAG( B( JROW, JCOL ) ) 150 CONTINUE 160 CONTINUE OPS = OPS + SOPBL2( 'SGEMV ', K, NRHS, 0, 0 ) CALL SGEMV( '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 ) = CMPLX( 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 OPS = OPS + REAL( 12*NRHS ) CALL CCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) CALL CSROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) END IF IF( K.LT.MAX( M, N ) ) $ CALL CLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, $ BX( K+1, 1 ), LDBX ) * * Step (3R): permute rows of B. * CALL CCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) IF( SQRE.EQ.1 ) THEN CALL CCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) END IF DO 190 I = 2, N CALL CCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) 190 CONTINUE * * Step (4R): apply back the Givens rotations performed. * OPS = OPS + REAL( 12*NRHS*GIVPTR ) DO 200 I = GIVPTR, 1, -1 CALL CSROT( 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 CLALS0 * END