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
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