Annotation of rpl/lapack/lapack/zlals0.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
! 2: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
! 3: $ POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
! 4: *
! 5: * -- LAPACK routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * .. Scalar Arguments ..
! 11: INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
! 12: $ LDGNUM, NL, NR, NRHS, SQRE
! 13: DOUBLE PRECISION C, S
! 14: * ..
! 15: * .. Array Arguments ..
! 16: INTEGER GIVCOL( LDGCOL, * ), PERM( * )
! 17: DOUBLE PRECISION DIFL( * ), DIFR( LDGNUM, * ),
! 18: $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
! 19: $ RWORK( * ), Z( * )
! 20: COMPLEX*16 B( LDB, * ), BX( LDBX, * )
! 21: * ..
! 22: *
! 23: * Purpose
! 24: * =======
! 25: *
! 26: * ZLALS0 applies back the multiplying factors of either the left or the
! 27: * right singular vector matrix of a diagonal matrix appended by a row
! 28: * to the right hand side matrix B in solving the least squares problem
! 29: * using the divide-and-conquer SVD approach.
! 30: *
! 31: * For the left singular vector matrix, three types of orthogonal
! 32: * matrices are involved:
! 33: *
! 34: * (1L) Givens rotations: the number of such rotations is GIVPTR; the
! 35: * pairs of columns/rows they were applied to are stored in GIVCOL;
! 36: * and the C- and S-values of these rotations are stored in GIVNUM.
! 37: *
! 38: * (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
! 39: * row, and for J=2:N, PERM(J)-th row of B is to be moved to the
! 40: * J-th row.
! 41: *
! 42: * (3L) The left singular vector matrix of the remaining matrix.
! 43: *
! 44: * For the right singular vector matrix, four types of orthogonal
! 45: * matrices are involved:
! 46: *
! 47: * (1R) The right singular vector matrix of the remaining matrix.
! 48: *
! 49: * (2R) If SQRE = 1, one extra Givens rotation to generate the right
! 50: * null space.
! 51: *
! 52: * (3R) The inverse transformation of (2L).
! 53: *
! 54: * (4R) The inverse transformation of (1L).
! 55: *
! 56: * Arguments
! 57: * =========
! 58: *
! 59: * ICOMPQ (input) INTEGER
! 60: * Specifies whether singular vectors are to be computed in
! 61: * factored form:
! 62: * = 0: Left singular vector matrix.
! 63: * = 1: Right singular vector matrix.
! 64: *
! 65: * NL (input) INTEGER
! 66: * The row dimension of the upper block. NL >= 1.
! 67: *
! 68: * NR (input) INTEGER
! 69: * The row dimension of the lower block. NR >= 1.
! 70: *
! 71: * SQRE (input) INTEGER
! 72: * = 0: the lower block is an NR-by-NR square matrix.
! 73: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
! 74: *
! 75: * The bidiagonal matrix has row dimension N = NL + NR + 1,
! 76: * and column dimension M = N + SQRE.
! 77: *
! 78: * NRHS (input) INTEGER
! 79: * The number of columns of B and BX. NRHS must be at least 1.
! 80: *
! 81: * B (input/output) COMPLEX*16 array, dimension ( LDB, NRHS )
! 82: * On input, B contains the right hand sides of the least
! 83: * squares problem in rows 1 through M. On output, B contains
! 84: * the solution X in rows 1 through N.
! 85: *
! 86: * LDB (input) INTEGER
! 87: * The leading dimension of B. LDB must be at least
! 88: * max(1,MAX( M, N ) ).
! 89: *
! 90: * BX (workspace) COMPLEX*16 array, dimension ( LDBX, NRHS )
! 91: *
! 92: * LDBX (input) INTEGER
! 93: * The leading dimension of BX.
! 94: *
! 95: * PERM (input) INTEGER array, dimension ( N )
! 96: * The permutations (from deflation and sorting) applied
! 97: * to the two blocks.
! 98: *
! 99: * GIVPTR (input) INTEGER
! 100: * The number of Givens rotations which took place in this
! 101: * subproblem.
! 102: *
! 103: * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
! 104: * Each pair of numbers indicates a pair of rows/columns
! 105: * involved in a Givens rotation.
! 106: *
! 107: * LDGCOL (input) INTEGER
! 108: * The leading dimension of GIVCOL, must be at least N.
! 109: *
! 110: * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 111: * Each number indicates the C or S value used in the
! 112: * corresponding Givens rotation.
! 113: *
! 114: * LDGNUM (input) INTEGER
! 115: * The leading dimension of arrays DIFR, POLES and
! 116: * GIVNUM, must be at least K.
! 117: *
! 118: * POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
! 119: * On entry, POLES(1:K, 1) contains the new singular
! 120: * values obtained from solving the secular equation, and
! 121: * POLES(1:K, 2) is an array containing the poles in the secular
! 122: * equation.
! 123: *
! 124: * DIFL (input) DOUBLE PRECISION array, dimension ( K ).
! 125: * On entry, DIFL(I) is the distance between I-th updated
! 126: * (undeflated) singular value and the I-th (undeflated) old
! 127: * singular value.
! 128: *
! 129: * DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
! 130: * On entry, DIFR(I, 1) contains the distances between I-th
! 131: * updated (undeflated) singular value and the I+1-th
! 132: * (undeflated) old singular value. And DIFR(I, 2) is the
! 133: * normalizing factor for the I-th right singular vector.
! 134: *
! 135: * Z (input) DOUBLE PRECISION array, dimension ( K )
! 136: * Contain the components of the deflation-adjusted updating row
! 137: * vector.
! 138: *
! 139: * K (input) INTEGER
! 140: * Contains the dimension of the non-deflated matrix,
! 141: * This is the order of the related secular equation. 1 <= K <=N.
! 142: *
! 143: * C (input) DOUBLE PRECISION
! 144: * C contains garbage if SQRE =0 and the C-value of a Givens
! 145: * rotation related to the right null space if SQRE = 1.
! 146: *
! 147: * S (input) DOUBLE PRECISION
! 148: * S contains garbage if SQRE =0 and the S-value of a Givens
! 149: * rotation related to the right null space if SQRE = 1.
! 150: *
! 151: * RWORK (workspace) DOUBLE PRECISION array, dimension
! 152: * ( K*(1+NRHS) + 2*NRHS )
! 153: *
! 154: * INFO (output) INTEGER
! 155: * = 0: successful exit.
! 156: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 157: *
! 158: * Further Details
! 159: * ===============
! 160: *
! 161: * Based on contributions by
! 162: * Ming Gu and Ren-Cang Li, Computer Science Division, University of
! 163: * California at Berkeley, USA
! 164: * Osni Marques, LBNL/NERSC, USA
! 165: *
! 166: * =====================================================================
! 167: *
! 168: * .. Parameters ..
! 169: DOUBLE PRECISION ONE, ZERO, NEGONE
! 170: PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
! 171: * ..
! 172: * .. Local Scalars ..
! 173: INTEGER I, J, JCOL, JROW, M, N, NLP1
! 174: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
! 175: * ..
! 176: * .. External Subroutines ..
! 177: EXTERNAL DGEMV, XERBLA, ZCOPY, ZDROT, ZDSCAL, ZLACPY,
! 178: $ ZLASCL
! 179: * ..
! 180: * .. External Functions ..
! 181: DOUBLE PRECISION DLAMC3, DNRM2
! 182: EXTERNAL DLAMC3, DNRM2
! 183: * ..
! 184: * .. Intrinsic Functions ..
! 185: INTRINSIC DBLE, DCMPLX, DIMAG, MAX
! 186: * ..
! 187: * .. Executable Statements ..
! 188: *
! 189: * Test the input parameters.
! 190: *
! 191: INFO = 0
! 192: *
! 193: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
! 194: INFO = -1
! 195: ELSE IF( NL.LT.1 ) THEN
! 196: INFO = -2
! 197: ELSE IF( NR.LT.1 ) THEN
! 198: INFO = -3
! 199: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
! 200: INFO = -4
! 201: END IF
! 202: *
! 203: N = NL + NR + 1
! 204: *
! 205: IF( NRHS.LT.1 ) THEN
! 206: INFO = -5
! 207: ELSE IF( LDB.LT.N ) THEN
! 208: INFO = -7
! 209: ELSE IF( LDBX.LT.N ) THEN
! 210: INFO = -9
! 211: ELSE IF( GIVPTR.LT.0 ) THEN
! 212: INFO = -11
! 213: ELSE IF( LDGCOL.LT.N ) THEN
! 214: INFO = -13
! 215: ELSE IF( LDGNUM.LT.N ) THEN
! 216: INFO = -15
! 217: ELSE IF( K.LT.1 ) THEN
! 218: INFO = -20
! 219: END IF
! 220: IF( INFO.NE.0 ) THEN
! 221: CALL XERBLA( 'ZLALS0', -INFO )
! 222: RETURN
! 223: END IF
! 224: *
! 225: M = N + SQRE
! 226: NLP1 = NL + 1
! 227: *
! 228: IF( ICOMPQ.EQ.0 ) THEN
! 229: *
! 230: * Apply back orthogonal transformations from the left.
! 231: *
! 232: * Step (1L): apply back the Givens rotations performed.
! 233: *
! 234: DO 10 I = 1, GIVPTR
! 235: CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
! 236: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
! 237: $ GIVNUM( I, 1 ) )
! 238: 10 CONTINUE
! 239: *
! 240: * Step (2L): permute rows of B.
! 241: *
! 242: CALL ZCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
! 243: DO 20 I = 2, N
! 244: CALL ZCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
! 245: 20 CONTINUE
! 246: *
! 247: * Step (3L): apply the inverse of the left singular vector
! 248: * matrix to BX.
! 249: *
! 250: IF( K.EQ.1 ) THEN
! 251: CALL ZCOPY( NRHS, BX, LDBX, B, LDB )
! 252: IF( Z( 1 ).LT.ZERO ) THEN
! 253: CALL ZDSCAL( NRHS, NEGONE, B, LDB )
! 254: END IF
! 255: ELSE
! 256: DO 100 J = 1, K
! 257: DIFLJ = DIFL( J )
! 258: DJ = POLES( J, 1 )
! 259: DSIGJ = -POLES( J, 2 )
! 260: IF( J.LT.K ) THEN
! 261: DIFRJ = -DIFR( J, 1 )
! 262: DSIGJP = -POLES( J+1, 2 )
! 263: END IF
! 264: IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
! 265: $ THEN
! 266: RWORK( J ) = ZERO
! 267: ELSE
! 268: RWORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
! 269: $ ( POLES( J, 2 )+DJ )
! 270: END IF
! 271: DO 30 I = 1, J - 1
! 272: IF( ( Z( I ).EQ.ZERO ) .OR.
! 273: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
! 274: RWORK( I ) = ZERO
! 275: ELSE
! 276: RWORK( I ) = POLES( I, 2 )*Z( I ) /
! 277: $ ( DLAMC3( POLES( I, 2 ), DSIGJ )-
! 278: $ DIFLJ ) / ( POLES( I, 2 )+DJ )
! 279: END IF
! 280: 30 CONTINUE
! 281: DO 40 I = J + 1, K
! 282: IF( ( Z( I ).EQ.ZERO ) .OR.
! 283: $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN
! 284: RWORK( I ) = ZERO
! 285: ELSE
! 286: RWORK( I ) = POLES( I, 2 )*Z( I ) /
! 287: $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+
! 288: $ DIFRJ ) / ( POLES( I, 2 )+DJ )
! 289: END IF
! 290: 40 CONTINUE
! 291: RWORK( 1 ) = NEGONE
! 292: TEMP = DNRM2( K, RWORK, 1 )
! 293: *
! 294: * Since B and BX are complex, the following call to DGEMV
! 295: * is performed in two steps (real and imaginary parts).
! 296: *
! 297: * CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
! 298: * $ B( J, 1 ), LDB )
! 299: *
! 300: I = K + NRHS*2
! 301: DO 60 JCOL = 1, NRHS
! 302: DO 50 JROW = 1, K
! 303: I = I + 1
! 304: RWORK( I ) = DBLE( BX( JROW, JCOL ) )
! 305: 50 CONTINUE
! 306: 60 CONTINUE
! 307: CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
! 308: $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
! 309: I = K + NRHS*2
! 310: DO 80 JCOL = 1, NRHS
! 311: DO 70 JROW = 1, K
! 312: I = I + 1
! 313: RWORK( I ) = DIMAG( BX( JROW, JCOL ) )
! 314: 70 CONTINUE
! 315: 80 CONTINUE
! 316: CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
! 317: $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
! 318: DO 90 JCOL = 1, NRHS
! 319: B( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
! 320: $ RWORK( JCOL+K+NRHS ) )
! 321: 90 CONTINUE
! 322: CALL ZLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
! 323: $ LDB, INFO )
! 324: 100 CONTINUE
! 325: END IF
! 326: *
! 327: * Move the deflated rows of BX to B also.
! 328: *
! 329: IF( K.LT.MAX( M, N ) )
! 330: $ CALL ZLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
! 331: $ B( K+1, 1 ), LDB )
! 332: ELSE
! 333: *
! 334: * Apply back the right orthogonal transformations.
! 335: *
! 336: * Step (1R): apply back the new right singular vector matrix
! 337: * to B.
! 338: *
! 339: IF( K.EQ.1 ) THEN
! 340: CALL ZCOPY( NRHS, B, LDB, BX, LDBX )
! 341: ELSE
! 342: DO 180 J = 1, K
! 343: DSIGJ = POLES( J, 2 )
! 344: IF( Z( J ).EQ.ZERO ) THEN
! 345: RWORK( J ) = ZERO
! 346: ELSE
! 347: RWORK( J ) = -Z( J ) / DIFL( J ) /
! 348: $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
! 349: END IF
! 350: DO 110 I = 1, J - 1
! 351: IF( Z( J ).EQ.ZERO ) THEN
! 352: RWORK( I ) = ZERO
! 353: ELSE
! 354: RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
! 355: $ 2 ) )-DIFR( I, 1 ) ) /
! 356: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
! 357: END IF
! 358: 110 CONTINUE
! 359: DO 120 I = J + 1, K
! 360: IF( Z( J ).EQ.ZERO ) THEN
! 361: RWORK( I ) = ZERO
! 362: ELSE
! 363: RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
! 364: $ 2 ) )-DIFL( I ) ) /
! 365: $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
! 366: END IF
! 367: 120 CONTINUE
! 368: *
! 369: * Since B and BX are complex, the following call to DGEMV
! 370: * is performed in two steps (real and imaginary parts).
! 371: *
! 372: * CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
! 373: * $ BX( J, 1 ), LDBX )
! 374: *
! 375: I = K + NRHS*2
! 376: DO 140 JCOL = 1, NRHS
! 377: DO 130 JROW = 1, K
! 378: I = I + 1
! 379: RWORK( I ) = DBLE( B( JROW, JCOL ) )
! 380: 130 CONTINUE
! 381: 140 CONTINUE
! 382: CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
! 383: $ RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
! 384: I = K + NRHS*2
! 385: DO 160 JCOL = 1, NRHS
! 386: DO 150 JROW = 1, K
! 387: I = I + 1
! 388: RWORK( I ) = DIMAG( B( JROW, JCOL ) )
! 389: 150 CONTINUE
! 390: 160 CONTINUE
! 391: CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
! 392: $ RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
! 393: DO 170 JCOL = 1, NRHS
! 394: BX( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
! 395: $ RWORK( JCOL+K+NRHS ) )
! 396: 170 CONTINUE
! 397: 180 CONTINUE
! 398: END IF
! 399: *
! 400: * Step (2R): if SQRE = 1, apply back the rotation that is
! 401: * related to the right null space of the subproblem.
! 402: *
! 403: IF( SQRE.EQ.1 ) THEN
! 404: CALL ZCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
! 405: CALL ZDROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
! 406: END IF
! 407: IF( K.LT.MAX( M, N ) )
! 408: $ CALL ZLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
! 409: $ LDBX )
! 410: *
! 411: * Step (3R): permute rows of B.
! 412: *
! 413: CALL ZCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
! 414: IF( SQRE.EQ.1 ) THEN
! 415: CALL ZCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
! 416: END IF
! 417: DO 190 I = 2, N
! 418: CALL ZCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
! 419: 190 CONTINUE
! 420: *
! 421: * Step (4R): apply back the Givens rotations performed.
! 422: *
! 423: DO 200 I = GIVPTR, 1, -1
! 424: CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
! 425: $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
! 426: $ -GIVNUM( I, 1 ) )
! 427: 200 CONTINUE
! 428: END IF
! 429: *
! 430: RETURN
! 431: *
! 432: * End of ZLALS0
! 433: *
! 434: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlals0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack