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