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