![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.2.2.
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.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: * June 2010 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 + 93: * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ), 94: * where 95: * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) 96: * 97: * IWORK (workspace) INTEGER array, dimension at least 98: * (3*N*NLVL + 11*N). 99: * 100: * INFO (output) INTEGER 101: * = 0: successful exit. 102: * < 0: if INFO = -i, the i-th argument had an illegal value. 103: * > 0: The algorithm failed to compute a singular value while 104: * working on the submatrix lying in rows and columns 105: * INFO/(N+1) through MOD(INFO,N+1). 106: * 107: * Further Details 108: * =============== 109: * 110: * Based on contributions by 111: * Ming Gu and Ren-Cang Li, Computer Science Division, University of 112: * California at Berkeley, USA 113: * Osni Marques, LBNL/NERSC, USA 114: * 115: * ===================================================================== 116: * 117: * .. Parameters .. 118: DOUBLE PRECISION ZERO, ONE, TWO 119: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 ) 120: COMPLEX*16 CZERO 121: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) ) 122: * .. 123: * .. Local Scalars .. 124: INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM, 125: $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB, 126: $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG, 127: $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB, 128: $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1, 129: $ U, VT, Z 130: DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL 131: * .. 132: * .. External Functions .. 133: INTEGER IDAMAX 134: DOUBLE PRECISION DLAMCH, DLANST 135: EXTERNAL IDAMAX, DLAMCH, DLANST 136: * .. 137: * .. External Subroutines .. 138: EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET, 139: $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA, 140: $ ZLASCL, ZLASET 141: * .. 142: * .. Intrinsic Functions .. 143: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN 144: * .. 145: * .. Executable Statements .. 146: * 147: * Test the input parameters. 148: * 149: INFO = 0 150: * 151: IF( N.LT.0 ) THEN 152: INFO = -3 153: ELSE IF( NRHS.LT.1 ) THEN 154: INFO = -4 155: ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN 156: INFO = -8 157: END IF 158: IF( INFO.NE.0 ) THEN 159: CALL XERBLA( 'ZLALSD', -INFO ) 160: RETURN 161: END IF 162: * 163: EPS = DLAMCH( 'Epsilon' ) 164: * 165: * Set up the tolerance. 166: * 167: IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN 168: RCND = EPS 169: ELSE 170: RCND = RCOND 171: END IF 172: * 173: RANK = 0 174: * 175: * Quick return if possible. 176: * 177: IF( N.EQ.0 ) THEN 178: RETURN 179: ELSE IF( N.EQ.1 ) THEN 180: IF( D( 1 ).EQ.ZERO ) THEN 181: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB ) 182: ELSE 183: RANK = 1 184: CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO ) 185: D( 1 ) = ABS( D( 1 ) ) 186: END IF 187: RETURN 188: END IF 189: * 190: * Rotate the matrix if it is lower bidiagonal. 191: * 192: IF( UPLO.EQ.'L' ) THEN 193: DO 10 I = 1, N - 1 194: CALL DLARTG( D( I ), E( I ), CS, SN, R ) 195: D( I ) = R 196: E( I ) = SN*D( I+1 ) 197: D( I+1 ) = CS*D( I+1 ) 198: IF( NRHS.EQ.1 ) THEN 199: CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN ) 200: ELSE 201: RWORK( I*2-1 ) = CS 202: RWORK( I*2 ) = SN 203: END IF 204: 10 CONTINUE 205: IF( NRHS.GT.1 ) THEN 206: DO 30 I = 1, NRHS 207: DO 20 J = 1, N - 1 208: CS = RWORK( J*2-1 ) 209: SN = RWORK( J*2 ) 210: CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN ) 211: 20 CONTINUE 212: 30 CONTINUE 213: END IF 214: END IF 215: * 216: * Scale. 217: * 218: NM1 = N - 1 219: ORGNRM = DLANST( 'M', N, D, E ) 220: IF( ORGNRM.EQ.ZERO ) THEN 221: CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB ) 222: RETURN 223: END IF 224: * 225: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 226: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO ) 227: * 228: * If N is smaller than the minimum divide size SMLSIZ, then solve 229: * the problem with another solver. 230: * 231: IF( N.LE.SMLSIZ ) THEN 232: IRWU = 1 233: IRWVT = IRWU + N*N 234: IRWWRK = IRWVT + N*N 235: IRWRB = IRWWRK 236: IRWIB = IRWRB + N*NRHS 237: IRWB = IRWIB + N*NRHS 238: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N ) 239: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N ) 240: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N, 241: $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1, 242: $ RWORK( IRWWRK ), INFO ) 243: IF( INFO.NE.0 ) THEN 244: RETURN 245: END IF 246: * 247: * In the real version, B is passed to DLASDQ and multiplied 248: * internally by Q'. Here B is complex and that product is 249: * computed below in two steps (real and imaginary parts). 250: * 251: J = IRWB - 1 252: DO 50 JCOL = 1, NRHS 253: DO 40 JROW = 1, N 254: J = J + 1 255: RWORK( J ) = DBLE( B( JROW, JCOL ) ) 256: 40 CONTINUE 257: 50 CONTINUE 258: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 259: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 260: J = IRWB - 1 261: DO 70 JCOL = 1, NRHS 262: DO 60 JROW = 1, N 263: J = J + 1 264: RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 265: 60 CONTINUE 266: 70 CONTINUE 267: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N, 268: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 269: JREAL = IRWRB - 1 270: JIMAG = IRWIB - 1 271: DO 90 JCOL = 1, NRHS 272: DO 80 JROW = 1, N 273: JREAL = JREAL + 1 274: JIMAG = JIMAG + 1 275: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 276: $ RWORK( JIMAG ) ) 277: 80 CONTINUE 278: 90 CONTINUE 279: * 280: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 281: DO 100 I = 1, N 282: IF( D( I ).LE.TOL ) THEN 283: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB ) 284: ELSE 285: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ), 286: $ LDB, INFO ) 287: RANK = RANK + 1 288: END IF 289: 100 CONTINUE 290: * 291: * Since B is complex, the following call to DGEMM is performed 292: * in two steps (real and imaginary parts). That is for V * B 293: * (in the real version of the code V' is stored in WORK). 294: * 295: * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO, 296: * $ WORK( NWORK ), N ) 297: * 298: J = IRWB - 1 299: DO 120 JCOL = 1, NRHS 300: DO 110 JROW = 1, N 301: J = J + 1 302: RWORK( J ) = DBLE( B( JROW, JCOL ) ) 303: 110 CONTINUE 304: 120 CONTINUE 305: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 306: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N ) 307: J = IRWB - 1 308: DO 140 JCOL = 1, NRHS 309: DO 130 JROW = 1, N 310: J = J + 1 311: RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 312: 130 CONTINUE 313: 140 CONTINUE 314: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N, 315: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N ) 316: JREAL = IRWRB - 1 317: JIMAG = IRWIB - 1 318: DO 160 JCOL = 1, NRHS 319: DO 150 JROW = 1, N 320: JREAL = JREAL + 1 321: JIMAG = JIMAG + 1 322: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 323: $ RWORK( JIMAG ) ) 324: 150 CONTINUE 325: 160 CONTINUE 326: * 327: * Unscale. 328: * 329: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 330: CALL DLASRT( 'D', N, D, INFO ) 331: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 332: * 333: RETURN 334: END IF 335: * 336: * Book-keeping and setting up some constants. 337: * 338: NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1 339: * 340: SMLSZP = SMLSIZ + 1 341: * 342: U = 1 343: VT = 1 + SMLSIZ*N 344: DIFL = VT + SMLSZP*N 345: DIFR = DIFL + NLVL*N 346: Z = DIFR + NLVL*N*2 347: C = Z + NLVL*N 348: S = C + N 349: POLES = S + N 350: GIVNUM = POLES + 2*NLVL*N 351: NRWORK = GIVNUM + 2*NLVL*N 352: BX = 1 353: * 354: IRWRB = NRWORK 355: IRWIB = IRWRB + SMLSIZ*NRHS 356: IRWB = IRWIB + SMLSIZ*NRHS 357: * 358: SIZEI = 1 + N 359: K = SIZEI + N 360: GIVPTR = K + N 361: PERM = GIVPTR + N 362: GIVCOL = PERM + NLVL*N 363: IWK = GIVCOL + NLVL*N*2 364: * 365: ST = 1 366: SQRE = 0 367: ICMPQ1 = 1 368: ICMPQ2 = 0 369: NSUB = 0 370: * 371: DO 170 I = 1, N 372: IF( ABS( D( I ) ).LT.EPS ) THEN 373: D( I ) = SIGN( EPS, D( I ) ) 374: END IF 375: 170 CONTINUE 376: * 377: DO 240 I = 1, NM1 378: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN 379: NSUB = NSUB + 1 380: IWORK( NSUB ) = ST 381: * 382: * Subproblem found. First determine its size and then 383: * apply divide and conquer on it. 384: * 385: IF( I.LT.NM1 ) THEN 386: * 387: * A subproblem with E(I) small for I < NM1. 388: * 389: NSIZE = I - ST + 1 390: IWORK( SIZEI+NSUB-1 ) = NSIZE 391: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN 392: * 393: * A subproblem with E(NM1) not too small but I = NM1. 394: * 395: NSIZE = N - ST + 1 396: IWORK( SIZEI+NSUB-1 ) = NSIZE 397: ELSE 398: * 399: * A subproblem with E(NM1) small. This implies an 400: * 1-by-1 subproblem at D(N), which is not solved 401: * explicitly. 402: * 403: NSIZE = I - ST + 1 404: IWORK( SIZEI+NSUB-1 ) = NSIZE 405: NSUB = NSUB + 1 406: IWORK( NSUB ) = N 407: IWORK( SIZEI+NSUB-1 ) = 1 408: CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N ) 409: END IF 410: ST1 = ST - 1 411: IF( NSIZE.EQ.1 ) THEN 412: * 413: * This is a 1-by-1 subproblem and is not solved 414: * explicitly. 415: * 416: CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N ) 417: ELSE IF( NSIZE.LE.SMLSIZ ) THEN 418: * 419: * This is a small subproblem and is solved by DLASDQ. 420: * 421: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 422: $ RWORK( VT+ST1 ), N ) 423: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE, 424: $ RWORK( U+ST1 ), N ) 425: CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ), 426: $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ), 427: $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ), 428: $ INFO ) 429: IF( INFO.NE.0 ) THEN 430: RETURN 431: END IF 432: * 433: * In the real version, B is passed to DLASDQ and multiplied 434: * internally by Q'. Here B is complex and that product is 435: * computed below in two steps (real and imaginary parts). 436: * 437: J = IRWB - 1 438: DO 190 JCOL = 1, NRHS 439: DO 180 JROW = ST, ST + NSIZE - 1 440: J = J + 1 441: RWORK( J ) = DBLE( B( JROW, JCOL ) ) 442: 180 CONTINUE 443: 190 CONTINUE 444: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 445: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 446: $ ZERO, RWORK( IRWRB ), NSIZE ) 447: J = IRWB - 1 448: DO 210 JCOL = 1, NRHS 449: DO 200 JROW = ST, ST + NSIZE - 1 450: J = J + 1 451: RWORK( J ) = DIMAG( B( JROW, JCOL ) ) 452: 200 CONTINUE 453: 210 CONTINUE 454: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 455: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE, 456: $ ZERO, RWORK( IRWIB ), NSIZE ) 457: JREAL = IRWRB - 1 458: JIMAG = IRWIB - 1 459: DO 230 JCOL = 1, NRHS 460: DO 220 JROW = ST, ST + NSIZE - 1 461: JREAL = JREAL + 1 462: JIMAG = JIMAG + 1 463: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 464: $ RWORK( JIMAG ) ) 465: 220 CONTINUE 466: 230 CONTINUE 467: * 468: CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB, 469: $ WORK( BX+ST1 ), N ) 470: ELSE 471: * 472: * A large problem. Solve it using divide and conquer. 473: * 474: CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ), 475: $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ), 476: $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ), 477: $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ), 478: $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ), 479: $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ), 480: $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ), 481: $ RWORK( S+ST1 ), RWORK( NRWORK ), 482: $ IWORK( IWK ), INFO ) 483: IF( INFO.NE.0 ) THEN 484: RETURN 485: END IF 486: BXST = BX + ST1 487: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ), 488: $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N, 489: $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 490: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 491: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 492: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 493: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 494: $ RWORK( C+ST1 ), RWORK( S+ST1 ), 495: $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 496: IF( INFO.NE.0 ) THEN 497: RETURN 498: END IF 499: END IF 500: ST = I + 1 501: END IF 502: 240 CONTINUE 503: * 504: * Apply the singular values and treat the tiny ones as zero. 505: * 506: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) ) 507: * 508: DO 250 I = 1, N 509: * 510: * Some of the elements in D can be negative because 1-by-1 511: * subproblems were not solved explicitly. 512: * 513: IF( ABS( D( I ) ).LE.TOL ) THEN 514: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N ) 515: ELSE 516: RANK = RANK + 1 517: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, 518: $ WORK( BX+I-1 ), N, INFO ) 519: END IF 520: D( I ) = ABS( D( I ) ) 521: 250 CONTINUE 522: * 523: * Now apply back the right singular vectors. 524: * 525: ICMPQ2 = 1 526: DO 320 I = 1, NSUB 527: ST = IWORK( I ) 528: ST1 = ST - 1 529: NSIZE = IWORK( SIZEI+I-1 ) 530: BXST = BX + ST1 531: IF( NSIZE.EQ.1 ) THEN 532: CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB ) 533: ELSE IF( NSIZE.LE.SMLSIZ ) THEN 534: * 535: * Since B and BX are complex, the following call to DGEMM 536: * is performed in two steps (real and imaginary parts). 537: * 538: * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 539: * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO, 540: * $ B( ST, 1 ), LDB ) 541: * 542: J = BXST - N - 1 543: JREAL = IRWB - 1 544: DO 270 JCOL = 1, NRHS 545: J = J + N 546: DO 260 JROW = 1, NSIZE 547: JREAL = JREAL + 1 548: RWORK( JREAL ) = DBLE( WORK( J+JROW ) ) 549: 260 CONTINUE 550: 270 CONTINUE 551: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 552: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 553: $ RWORK( IRWRB ), NSIZE ) 554: J = BXST - N - 1 555: JIMAG = IRWB - 1 556: DO 290 JCOL = 1, NRHS 557: J = J + N 558: DO 280 JROW = 1, NSIZE 559: JIMAG = JIMAG + 1 560: RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) ) 561: 280 CONTINUE 562: 290 CONTINUE 563: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE, 564: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO, 565: $ RWORK( IRWIB ), NSIZE ) 566: JREAL = IRWRB - 1 567: JIMAG = IRWIB - 1 568: DO 310 JCOL = 1, NRHS 569: DO 300 JROW = ST, ST + NSIZE - 1 570: JREAL = JREAL + 1 571: JIMAG = JIMAG + 1 572: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ), 573: $ RWORK( JIMAG ) ) 574: 300 CONTINUE 575: 310 CONTINUE 576: ELSE 577: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N, 578: $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N, 579: $ RWORK( VT+ST1 ), IWORK( K+ST1 ), 580: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ), 581: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ), 582: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N, 583: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ), 584: $ RWORK( C+ST1 ), RWORK( S+ST1 ), 585: $ RWORK( NRWORK ), IWORK( IWK ), INFO ) 586: IF( INFO.NE.0 ) THEN 587: RETURN 588: END IF 589: END IF 590: 320 CONTINUE 591: * 592: * Unscale and sort the singular values. 593: * 594: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 595: CALL DLASRT( 'D', N, D, INFO ) 596: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO ) 597: * 598: RETURN 599: * 600: * End of ZLALSD 601: * 602: END