Annotation of rpl/lapack/lapack/zlalsd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
! 2: $ RANK, WORK, RWORK, IWORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER UPLO
! 11: INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
! 12: DOUBLE PRECISION RCOND
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER IWORK( * )
! 16: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
! 17: COMPLEX*16 B( LDB, * ), WORK( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * ZLALSD uses the singular value decomposition of A to solve the least
! 24: * squares problem of finding X to minimize the Euclidean norm of each
! 25: * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
! 26: * are N-by-NRHS. The solution X overwrites B.
! 27: *
! 28: * The singular values of A smaller than RCOND times the largest
! 29: * singular value are treated as zero in solving the least squares
! 30: * problem; in this case a minimum norm solution is returned.
! 31: * The actual singular values are returned in D in ascending order.
! 32: *
! 33: * This code makes very mild assumptions about floating point
! 34: * arithmetic. It will work on machines with a guard digit in
! 35: * add/subtract, or on those binary machines without guard digits
! 36: * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
! 37: * It could conceivably fail on hexadecimal or decimal machines
! 38: * without guard digits, but we know of none.
! 39: *
! 40: * Arguments
! 41: * =========
! 42: *
! 43: * UPLO (input) CHARACTER*1
! 44: * = 'U': D and E define an upper bidiagonal matrix.
! 45: * = 'L': D and E define a lower bidiagonal matrix.
! 46: *
! 47: * SMLSIZ (input) INTEGER
! 48: * The maximum size of the subproblems at the bottom of the
! 49: * computation tree.
! 50: *
! 51: * N (input) INTEGER
! 52: * The dimension of the bidiagonal matrix. N >= 0.
! 53: *
! 54: * NRHS (input) INTEGER
! 55: * The number of columns of B. NRHS must be at least 1.
! 56: *
! 57: * D (input/output) DOUBLE PRECISION array, dimension (N)
! 58: * On entry D contains the main diagonal of the bidiagonal
! 59: * matrix. On exit, if INFO = 0, D contains its singular values.
! 60: *
! 61: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
! 62: * Contains the super-diagonal entries of the bidiagonal matrix.
! 63: * On exit, E has been destroyed.
! 64: *
! 65: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
! 66: * On input, B contains the right hand sides of the least
! 67: * squares problem. On output, B contains the solution X.
! 68: *
! 69: * LDB (input) INTEGER
! 70: * The leading dimension of B in the calling subprogram.
! 71: * LDB must be at least max(1,N).
! 72: *
! 73: * RCOND (input) DOUBLE PRECISION
! 74: * The singular values of A less than or equal to RCOND times
! 75: * the largest singular value are treated as zero in solving
! 76: * the least squares problem. If RCOND is negative,
! 77: * machine precision is used instead.
! 78: * For example, if diag(S)*X=B were the least squares problem,
! 79: * where diag(S) is a diagonal matrix of singular values, the
! 80: * solution would be X(i) = B(i) / S(i) if S(i) is greater than
! 81: * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
! 82: * RCOND*max(S).
! 83: *
! 84: * RANK (output) INTEGER
! 85: * The number of singular values of A greater than RCOND times
! 86: * the largest singular value.
! 87: *
! 88: * WORK (workspace) COMPLEX*16 array, dimension at least
! 89: * (N * NRHS).
! 90: *
! 91: * RWORK (workspace) DOUBLE PRECISION array, dimension at least
! 92: * (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2),
! 93: * where
! 94: * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
! 95: *
! 96: * IWORK (workspace) INTEGER array, dimension at least
! 97: * (3*N*NLVL + 11*N).
! 98: *
! 99: * INFO (output) INTEGER
! 100: * = 0: successful exit.
! 101: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 102: * > 0: The algorithm failed to compute an singular value while
! 103: * working on the submatrix lying in rows and columns
! 104: * INFO/(N+1) through MOD(INFO,N+1).
! 105: *
! 106: * Further Details
! 107: * ===============
! 108: *
! 109: * Based on contributions by
! 110: * Ming Gu and Ren-Cang Li, Computer Science Division, University of
! 111: * California at Berkeley, USA
! 112: * Osni Marques, LBNL/NERSC, USA
! 113: *
! 114: * =====================================================================
! 115: *
! 116: * .. Parameters ..
! 117: DOUBLE PRECISION ZERO, ONE, TWO
! 118: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
! 119: COMPLEX*16 CZERO
! 120: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) )
! 121: * ..
! 122: * .. Local Scalars ..
! 123: INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
! 124: $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
! 125: $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
! 126: $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
! 127: $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
! 128: $ U, VT, Z
! 129: DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
! 130: * ..
! 131: * .. External Functions ..
! 132: INTEGER IDAMAX
! 133: DOUBLE PRECISION DLAMCH, DLANST
! 134: EXTERNAL IDAMAX, DLAMCH, DLANST
! 135: * ..
! 136: * .. External Subroutines ..
! 137: EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
! 138: $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
! 139: $ ZLASCL, ZLASET
! 140: * ..
! 141: * .. Intrinsic Functions ..
! 142: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
! 143: * ..
! 144: * .. Executable Statements ..
! 145: *
! 146: * Test the input parameters.
! 147: *
! 148: INFO = 0
! 149: *
! 150: IF( N.LT.0 ) THEN
! 151: INFO = -3
! 152: ELSE IF( NRHS.LT.1 ) THEN
! 153: INFO = -4
! 154: ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
! 155: INFO = -8
! 156: END IF
! 157: IF( INFO.NE.0 ) THEN
! 158: CALL XERBLA( 'ZLALSD', -INFO )
! 159: RETURN
! 160: END IF
! 161: *
! 162: EPS = DLAMCH( 'Epsilon' )
! 163: *
! 164: * Set up the tolerance.
! 165: *
! 166: IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
! 167: RCND = EPS
! 168: ELSE
! 169: RCND = RCOND
! 170: END IF
! 171: *
! 172: RANK = 0
! 173: *
! 174: * Quick return if possible.
! 175: *
! 176: IF( N.EQ.0 ) THEN
! 177: RETURN
! 178: ELSE IF( N.EQ.1 ) THEN
! 179: IF( D( 1 ).EQ.ZERO ) THEN
! 180: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
! 181: ELSE
! 182: RANK = 1
! 183: CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
! 184: D( 1 ) = ABS( D( 1 ) )
! 185: END IF
! 186: RETURN
! 187: END IF
! 188: *
! 189: * Rotate the matrix if it is lower bidiagonal.
! 190: *
! 191: IF( UPLO.EQ.'L' ) THEN
! 192: DO 10 I = 1, N - 1
! 193: CALL DLARTG( D( I ), E( I ), CS, SN, R )
! 194: D( I ) = R
! 195: E( I ) = SN*D( I+1 )
! 196: D( I+1 ) = CS*D( I+1 )
! 197: IF( NRHS.EQ.1 ) THEN
! 198: CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
! 199: ELSE
! 200: RWORK( I*2-1 ) = CS
! 201: RWORK( I*2 ) = SN
! 202: END IF
! 203: 10 CONTINUE
! 204: IF( NRHS.GT.1 ) THEN
! 205: DO 30 I = 1, NRHS
! 206: DO 20 J = 1, N - 1
! 207: CS = RWORK( J*2-1 )
! 208: SN = RWORK( J*2 )
! 209: CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
! 210: 20 CONTINUE
! 211: 30 CONTINUE
! 212: END IF
! 213: END IF
! 214: *
! 215: * Scale.
! 216: *
! 217: NM1 = N - 1
! 218: ORGNRM = DLANST( 'M', N, D, E )
! 219: IF( ORGNRM.EQ.ZERO ) THEN
! 220: CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
! 221: RETURN
! 222: END IF
! 223: *
! 224: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
! 225: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
! 226: *
! 227: * If N is smaller than the minimum divide size SMLSIZ, then solve
! 228: * the problem with another solver.
! 229: *
! 230: IF( N.LE.SMLSIZ ) THEN
! 231: IRWU = 1
! 232: IRWVT = IRWU + N*N
! 233: IRWWRK = IRWVT + N*N
! 234: IRWRB = IRWWRK
! 235: IRWIB = IRWRB + N*NRHS
! 236: IRWB = IRWIB + N*NRHS
! 237: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
! 238: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
! 239: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
! 240: $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
! 241: $ RWORK( IRWWRK ), INFO )
! 242: IF( INFO.NE.0 ) THEN
! 243: RETURN
! 244: END IF
! 245: *
! 246: * In the real version, B is passed to DLASDQ and multiplied
! 247: * internally by Q'. Here B is complex and that product is
! 248: * computed below in two steps (real and imaginary parts).
! 249: *
! 250: J = IRWB - 1
! 251: DO 50 JCOL = 1, NRHS
! 252: DO 40 JROW = 1, N
! 253: J = J + 1
! 254: RWORK( J ) = DBLE( B( JROW, JCOL ) )
! 255: 40 CONTINUE
! 256: 50 CONTINUE
! 257: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
! 258: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
! 259: J = IRWB - 1
! 260: DO 70 JCOL = 1, NRHS
! 261: DO 60 JROW = 1, N
! 262: J = J + 1
! 263: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
! 264: 60 CONTINUE
! 265: 70 CONTINUE
! 266: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
! 267: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
! 268: JREAL = IRWRB - 1
! 269: JIMAG = IRWIB - 1
! 270: DO 90 JCOL = 1, NRHS
! 271: DO 80 JROW = 1, N
! 272: JREAL = JREAL + 1
! 273: JIMAG = JIMAG + 1
! 274: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
! 275: $ RWORK( JIMAG ) )
! 276: 80 CONTINUE
! 277: 90 CONTINUE
! 278: *
! 279: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
! 280: DO 100 I = 1, N
! 281: IF( D( I ).LE.TOL ) THEN
! 282: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
! 283: ELSE
! 284: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
! 285: $ LDB, INFO )
! 286: RANK = RANK + 1
! 287: END IF
! 288: 100 CONTINUE
! 289: *
! 290: * Since B is complex, the following call to DGEMM is performed
! 291: * in two steps (real and imaginary parts). That is for V * B
! 292: * (in the real version of the code V' is stored in WORK).
! 293: *
! 294: * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
! 295: * $ WORK( NWORK ), N )
! 296: *
! 297: J = IRWB - 1
! 298: DO 120 JCOL = 1, NRHS
! 299: DO 110 JROW = 1, N
! 300: J = J + 1
! 301: RWORK( J ) = DBLE( B( JROW, JCOL ) )
! 302: 110 CONTINUE
! 303: 120 CONTINUE
! 304: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
! 305: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
! 306: J = IRWB - 1
! 307: DO 140 JCOL = 1, NRHS
! 308: DO 130 JROW = 1, N
! 309: J = J + 1
! 310: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
! 311: 130 CONTINUE
! 312: 140 CONTINUE
! 313: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
! 314: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
! 315: JREAL = IRWRB - 1
! 316: JIMAG = IRWIB - 1
! 317: DO 160 JCOL = 1, NRHS
! 318: DO 150 JROW = 1, N
! 319: JREAL = JREAL + 1
! 320: JIMAG = JIMAG + 1
! 321: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
! 322: $ RWORK( JIMAG ) )
! 323: 150 CONTINUE
! 324: 160 CONTINUE
! 325: *
! 326: * Unscale.
! 327: *
! 328: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
! 329: CALL DLASRT( 'D', N, D, INFO )
! 330: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
! 331: *
! 332: RETURN
! 333: END IF
! 334: *
! 335: * Book-keeping and setting up some constants.
! 336: *
! 337: NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
! 338: *
! 339: SMLSZP = SMLSIZ + 1
! 340: *
! 341: U = 1
! 342: VT = 1 + SMLSIZ*N
! 343: DIFL = VT + SMLSZP*N
! 344: DIFR = DIFL + NLVL*N
! 345: Z = DIFR + NLVL*N*2
! 346: C = Z + NLVL*N
! 347: S = C + N
! 348: POLES = S + N
! 349: GIVNUM = POLES + 2*NLVL*N
! 350: NRWORK = GIVNUM + 2*NLVL*N
! 351: BX = 1
! 352: *
! 353: IRWRB = NRWORK
! 354: IRWIB = IRWRB + SMLSIZ*NRHS
! 355: IRWB = IRWIB + SMLSIZ*NRHS
! 356: *
! 357: SIZEI = 1 + N
! 358: K = SIZEI + N
! 359: GIVPTR = K + N
! 360: PERM = GIVPTR + N
! 361: GIVCOL = PERM + NLVL*N
! 362: IWK = GIVCOL + NLVL*N*2
! 363: *
! 364: ST = 1
! 365: SQRE = 0
! 366: ICMPQ1 = 1
! 367: ICMPQ2 = 0
! 368: NSUB = 0
! 369: *
! 370: DO 170 I = 1, N
! 371: IF( ABS( D( I ) ).LT.EPS ) THEN
! 372: D( I ) = SIGN( EPS, D( I ) )
! 373: END IF
! 374: 170 CONTINUE
! 375: *
! 376: DO 240 I = 1, NM1
! 377: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
! 378: NSUB = NSUB + 1
! 379: IWORK( NSUB ) = ST
! 380: *
! 381: * Subproblem found. First determine its size and then
! 382: * apply divide and conquer on it.
! 383: *
! 384: IF( I.LT.NM1 ) THEN
! 385: *
! 386: * A subproblem with E(I) small for I < NM1.
! 387: *
! 388: NSIZE = I - ST + 1
! 389: IWORK( SIZEI+NSUB-1 ) = NSIZE
! 390: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
! 391: *
! 392: * A subproblem with E(NM1) not too small but I = NM1.
! 393: *
! 394: NSIZE = N - ST + 1
! 395: IWORK( SIZEI+NSUB-1 ) = NSIZE
! 396: ELSE
! 397: *
! 398: * A subproblem with E(NM1) small. This implies an
! 399: * 1-by-1 subproblem at D(N), which is not solved
! 400: * explicitly.
! 401: *
! 402: NSIZE = I - ST + 1
! 403: IWORK( SIZEI+NSUB-1 ) = NSIZE
! 404: NSUB = NSUB + 1
! 405: IWORK( NSUB ) = N
! 406: IWORK( SIZEI+NSUB-1 ) = 1
! 407: CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
! 408: END IF
! 409: ST1 = ST - 1
! 410: IF( NSIZE.EQ.1 ) THEN
! 411: *
! 412: * This is a 1-by-1 subproblem and is not solved
! 413: * explicitly.
! 414: *
! 415: CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
! 416: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
! 417: *
! 418: * This is a small subproblem and is solved by DLASDQ.
! 419: *
! 420: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
! 421: $ RWORK( VT+ST1 ), N )
! 422: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
! 423: $ RWORK( U+ST1 ), N )
! 424: CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
! 425: $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
! 426: $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
! 427: $ INFO )
! 428: IF( INFO.NE.0 ) THEN
! 429: RETURN
! 430: END IF
! 431: *
! 432: * In the real version, B is passed to DLASDQ and multiplied
! 433: * internally by Q'. Here B is complex and that product is
! 434: * computed below in two steps (real and imaginary parts).
! 435: *
! 436: J = IRWB - 1
! 437: DO 190 JCOL = 1, NRHS
! 438: DO 180 JROW = ST, ST + NSIZE - 1
! 439: J = J + 1
! 440: RWORK( J ) = DBLE( B( JROW, JCOL ) )
! 441: 180 CONTINUE
! 442: 190 CONTINUE
! 443: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
! 444: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
! 445: $ ZERO, RWORK( IRWRB ), NSIZE )
! 446: J = IRWB - 1
! 447: DO 210 JCOL = 1, NRHS
! 448: DO 200 JROW = ST, ST + NSIZE - 1
! 449: J = J + 1
! 450: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
! 451: 200 CONTINUE
! 452: 210 CONTINUE
! 453: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
! 454: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
! 455: $ ZERO, RWORK( IRWIB ), NSIZE )
! 456: JREAL = IRWRB - 1
! 457: JIMAG = IRWIB - 1
! 458: DO 230 JCOL = 1, NRHS
! 459: DO 220 JROW = ST, ST + NSIZE - 1
! 460: JREAL = JREAL + 1
! 461: JIMAG = JIMAG + 1
! 462: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
! 463: $ RWORK( JIMAG ) )
! 464: 220 CONTINUE
! 465: 230 CONTINUE
! 466: *
! 467: CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
! 468: $ WORK( BX+ST1 ), N )
! 469: ELSE
! 470: *
! 471: * A large problem. Solve it using divide and conquer.
! 472: *
! 473: CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
! 474: $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
! 475: $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
! 476: $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
! 477: $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
! 478: $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
! 479: $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
! 480: $ RWORK( S+ST1 ), RWORK( NRWORK ),
! 481: $ IWORK( IWK ), INFO )
! 482: IF( INFO.NE.0 ) THEN
! 483: RETURN
! 484: END IF
! 485: BXST = BX + ST1
! 486: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
! 487: $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
! 488: $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
! 489: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
! 490: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
! 491: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
! 492: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
! 493: $ RWORK( C+ST1 ), RWORK( S+ST1 ),
! 494: $ RWORK( NRWORK ), IWORK( IWK ), INFO )
! 495: IF( INFO.NE.0 ) THEN
! 496: RETURN
! 497: END IF
! 498: END IF
! 499: ST = I + 1
! 500: END IF
! 501: 240 CONTINUE
! 502: *
! 503: * Apply the singular values and treat the tiny ones as zero.
! 504: *
! 505: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
! 506: *
! 507: DO 250 I = 1, N
! 508: *
! 509: * Some of the elements in D can be negative because 1-by-1
! 510: * subproblems were not solved explicitly.
! 511: *
! 512: IF( ABS( D( I ) ).LE.TOL ) THEN
! 513: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
! 514: ELSE
! 515: RANK = RANK + 1
! 516: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
! 517: $ WORK( BX+I-1 ), N, INFO )
! 518: END IF
! 519: D( I ) = ABS( D( I ) )
! 520: 250 CONTINUE
! 521: *
! 522: * Now apply back the right singular vectors.
! 523: *
! 524: ICMPQ2 = 1
! 525: DO 320 I = 1, NSUB
! 526: ST = IWORK( I )
! 527: ST1 = ST - 1
! 528: NSIZE = IWORK( SIZEI+I-1 )
! 529: BXST = BX + ST1
! 530: IF( NSIZE.EQ.1 ) THEN
! 531: CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
! 532: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
! 533: *
! 534: * Since B and BX are complex, the following call to DGEMM
! 535: * is performed in two steps (real and imaginary parts).
! 536: *
! 537: * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
! 538: * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
! 539: * $ B( ST, 1 ), LDB )
! 540: *
! 541: J = BXST - N - 1
! 542: JREAL = IRWB - 1
! 543: DO 270 JCOL = 1, NRHS
! 544: J = J + N
! 545: DO 260 JROW = 1, NSIZE
! 546: JREAL = JREAL + 1
! 547: RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
! 548: 260 CONTINUE
! 549: 270 CONTINUE
! 550: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
! 551: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
! 552: $ RWORK( IRWRB ), NSIZE )
! 553: J = BXST - N - 1
! 554: JIMAG = IRWB - 1
! 555: DO 290 JCOL = 1, NRHS
! 556: J = J + N
! 557: DO 280 JROW = 1, NSIZE
! 558: JIMAG = JIMAG + 1
! 559: RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
! 560: 280 CONTINUE
! 561: 290 CONTINUE
! 562: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
! 563: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
! 564: $ RWORK( IRWIB ), NSIZE )
! 565: JREAL = IRWRB - 1
! 566: JIMAG = IRWIB - 1
! 567: DO 310 JCOL = 1, NRHS
! 568: DO 300 JROW = ST, ST + NSIZE - 1
! 569: JREAL = JREAL + 1
! 570: JIMAG = JIMAG + 1
! 571: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
! 572: $ RWORK( JIMAG ) )
! 573: 300 CONTINUE
! 574: 310 CONTINUE
! 575: ELSE
! 576: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
! 577: $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
! 578: $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
! 579: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
! 580: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
! 581: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
! 582: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
! 583: $ RWORK( C+ST1 ), RWORK( S+ST1 ),
! 584: $ RWORK( NRWORK ), IWORK( IWK ), INFO )
! 585: IF( INFO.NE.0 ) THEN
! 586: RETURN
! 587: END IF
! 588: END IF
! 589: 320 CONTINUE
! 590: *
! 591: * Unscale and sort the singular values.
! 592: *
! 593: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
! 594: CALL DLASRT( 'D', N, D, INFO )
! 595: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
! 596: *
! 597: RETURN
! 598: *
! 599: * End of ZLALSD
! 600: *
! 601: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlalsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack