Annotation of rpl/lapack/lapack/zlalsd.f, revision 1.21
1.13 bertrand 1: *> \brief \b ZLALSD uses the singular value decomposition of A to solve the least squares problem.
1.10 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.10 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 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">
1.10 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.10 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
22: * RANK, WORK, RWORK, IWORK, INFO )
1.17 bertrand 23: *
1.10 bertrand 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: * ..
1.17 bertrand 34: *
1.10 bertrand 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
1.19 bertrand 139: *> WORK is COMPLEX*16 array, dimension (N * NRHS)
1.10 bertrand 140: *> \endverbatim
141: *>
142: *> \param[out] RWORK
143: *> \verbatim
144: *> RWORK is DOUBLE PRECISION array, dimension at least
145: *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
146: *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
147: *> where
148: *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
149: *> \endverbatim
150: *>
151: *> \param[out] IWORK
152: *> \verbatim
153: *> IWORK is INTEGER array, dimension at least
154: *> (3*N*NLVL + 11*N).
155: *> \endverbatim
156: *>
157: *> \param[out] INFO
158: *> \verbatim
159: *> INFO is INTEGER
160: *> = 0: successful exit.
161: *> < 0: if INFO = -i, the i-th argument had an illegal value.
162: *> > 0: The algorithm failed to compute a singular value while
163: *> working on the submatrix lying in rows and columns
164: *> INFO/(N+1) through MOD(INFO,N+1).
165: *> \endverbatim
166: *
167: * Authors:
168: * ========
169: *
1.17 bertrand 170: *> \author Univ. of Tennessee
171: *> \author Univ. of California Berkeley
172: *> \author Univ. of Colorado Denver
173: *> \author NAG Ltd.
1.10 bertrand 174: *
175: *> \ingroup complex16OTHERcomputational
176: *
177: *> \par Contributors:
178: * ==================
179: *>
180: *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
181: *> California at Berkeley, USA \n
182: *> Osni Marques, LBNL/NERSC, USA \n
183: *
184: * =====================================================================
1.1 bertrand 185: SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
186: $ RANK, WORK, RWORK, IWORK, INFO )
187: *
1.21 ! bertrand 188: * -- LAPACK computational routine --
1.1 bertrand 189: * -- LAPACK is a software package provided by Univ. of Tennessee, --
190: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191: *
192: * .. Scalar Arguments ..
193: CHARACTER UPLO
194: INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
195: DOUBLE PRECISION RCOND
196: * ..
197: * .. Array Arguments ..
198: INTEGER IWORK( * )
199: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
200: COMPLEX*16 B( LDB, * ), WORK( * )
201: * ..
202: *
203: * =====================================================================
204: *
205: * .. Parameters ..
206: DOUBLE PRECISION ZERO, ONE, TWO
207: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
208: COMPLEX*16 CZERO
209: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) )
210: * ..
211: * .. Local Scalars ..
212: INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
213: $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
214: $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
215: $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
216: $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
217: $ U, VT, Z
218: DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
219: * ..
220: * .. External Functions ..
221: INTEGER IDAMAX
222: DOUBLE PRECISION DLAMCH, DLANST
223: EXTERNAL IDAMAX, DLAMCH, DLANST
224: * ..
225: * .. External Subroutines ..
226: EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
227: $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
228: $ ZLASCL, ZLASET
229: * ..
230: * .. Intrinsic Functions ..
231: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
232: * ..
233: * .. Executable Statements ..
234: *
235: * Test the input parameters.
236: *
237: INFO = 0
238: *
239: IF( N.LT.0 ) THEN
240: INFO = -3
241: ELSE IF( NRHS.LT.1 ) THEN
242: INFO = -4
243: ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
244: INFO = -8
245: END IF
246: IF( INFO.NE.0 ) THEN
247: CALL XERBLA( 'ZLALSD', -INFO )
248: RETURN
249: END IF
250: *
251: EPS = DLAMCH( 'Epsilon' )
252: *
253: * Set up the tolerance.
254: *
255: IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
256: RCND = EPS
257: ELSE
258: RCND = RCOND
259: END IF
260: *
261: RANK = 0
262: *
263: * Quick return if possible.
264: *
265: IF( N.EQ.0 ) THEN
266: RETURN
267: ELSE IF( N.EQ.1 ) THEN
268: IF( D( 1 ).EQ.ZERO ) THEN
269: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
270: ELSE
271: RANK = 1
272: CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
273: D( 1 ) = ABS( D( 1 ) )
274: END IF
275: RETURN
276: END IF
277: *
278: * Rotate the matrix if it is lower bidiagonal.
279: *
280: IF( UPLO.EQ.'L' ) THEN
281: DO 10 I = 1, N - 1
282: CALL DLARTG( D( I ), E( I ), CS, SN, R )
283: D( I ) = R
284: E( I ) = SN*D( I+1 )
285: D( I+1 ) = CS*D( I+1 )
286: IF( NRHS.EQ.1 ) THEN
287: CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
288: ELSE
289: RWORK( I*2-1 ) = CS
290: RWORK( I*2 ) = SN
291: END IF
292: 10 CONTINUE
293: IF( NRHS.GT.1 ) THEN
294: DO 30 I = 1, NRHS
295: DO 20 J = 1, N - 1
296: CS = RWORK( J*2-1 )
297: SN = RWORK( J*2 )
298: CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
299: 20 CONTINUE
300: 30 CONTINUE
301: END IF
302: END IF
303: *
304: * Scale.
305: *
306: NM1 = N - 1
307: ORGNRM = DLANST( 'M', N, D, E )
308: IF( ORGNRM.EQ.ZERO ) THEN
309: CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
310: RETURN
311: END IF
312: *
313: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
314: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
315: *
316: * If N is smaller than the minimum divide size SMLSIZ, then solve
317: * the problem with another solver.
318: *
319: IF( N.LE.SMLSIZ ) THEN
320: IRWU = 1
321: IRWVT = IRWU + N*N
322: IRWWRK = IRWVT + N*N
323: IRWRB = IRWWRK
324: IRWIB = IRWRB + N*NRHS
325: IRWB = IRWIB + N*NRHS
326: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
327: CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
328: CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
329: $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
330: $ RWORK( IRWWRK ), INFO )
331: IF( INFO.NE.0 ) THEN
332: RETURN
333: END IF
334: *
335: * In the real version, B is passed to DLASDQ and multiplied
1.9 bertrand 336: * internally by Q**H. Here B is complex and that product is
1.1 bertrand 337: * computed below in two steps (real and imaginary parts).
338: *
339: J = IRWB - 1
340: DO 50 JCOL = 1, NRHS
341: DO 40 JROW = 1, N
342: J = J + 1
343: RWORK( J ) = DBLE( B( JROW, JCOL ) )
344: 40 CONTINUE
345: 50 CONTINUE
346: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
347: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
348: J = IRWB - 1
349: DO 70 JCOL = 1, NRHS
350: DO 60 JROW = 1, N
351: J = J + 1
352: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
353: 60 CONTINUE
354: 70 CONTINUE
355: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
356: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
357: JREAL = IRWRB - 1
358: JIMAG = IRWIB - 1
359: DO 90 JCOL = 1, NRHS
360: DO 80 JROW = 1, N
361: JREAL = JREAL + 1
362: JIMAG = JIMAG + 1
363: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
364: $ RWORK( JIMAG ) )
365: 80 CONTINUE
366: 90 CONTINUE
367: *
368: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
369: DO 100 I = 1, N
370: IF( D( I ).LE.TOL ) THEN
371: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
372: ELSE
373: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
374: $ LDB, INFO )
375: RANK = RANK + 1
376: END IF
377: 100 CONTINUE
378: *
379: * Since B is complex, the following call to DGEMM is performed
380: * in two steps (real and imaginary parts). That is for V * B
1.9 bertrand 381: * (in the real version of the code V**H is stored in WORK).
1.1 bertrand 382: *
383: * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
384: * $ WORK( NWORK ), N )
385: *
386: J = IRWB - 1
387: DO 120 JCOL = 1, NRHS
388: DO 110 JROW = 1, N
389: J = J + 1
390: RWORK( J ) = DBLE( B( JROW, JCOL ) )
391: 110 CONTINUE
392: 120 CONTINUE
393: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
394: $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
395: J = IRWB - 1
396: DO 140 JCOL = 1, NRHS
397: DO 130 JROW = 1, N
398: J = J + 1
399: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
400: 130 CONTINUE
401: 140 CONTINUE
402: CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
403: $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
404: JREAL = IRWRB - 1
405: JIMAG = IRWIB - 1
406: DO 160 JCOL = 1, NRHS
407: DO 150 JROW = 1, N
408: JREAL = JREAL + 1
409: JIMAG = JIMAG + 1
410: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
411: $ RWORK( JIMAG ) )
412: 150 CONTINUE
413: 160 CONTINUE
414: *
415: * Unscale.
416: *
417: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
418: CALL DLASRT( 'D', N, D, INFO )
419: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
420: *
421: RETURN
422: END IF
423: *
424: * Book-keeping and setting up some constants.
425: *
426: NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
427: *
428: SMLSZP = SMLSIZ + 1
429: *
430: U = 1
431: VT = 1 + SMLSIZ*N
432: DIFL = VT + SMLSZP*N
433: DIFR = DIFL + NLVL*N
434: Z = DIFR + NLVL*N*2
435: C = Z + NLVL*N
436: S = C + N
437: POLES = S + N
438: GIVNUM = POLES + 2*NLVL*N
439: NRWORK = GIVNUM + 2*NLVL*N
440: BX = 1
441: *
442: IRWRB = NRWORK
443: IRWIB = IRWRB + SMLSIZ*NRHS
444: IRWB = IRWIB + SMLSIZ*NRHS
445: *
446: SIZEI = 1 + N
447: K = SIZEI + N
448: GIVPTR = K + N
449: PERM = GIVPTR + N
450: GIVCOL = PERM + NLVL*N
451: IWK = GIVCOL + NLVL*N*2
452: *
453: ST = 1
454: SQRE = 0
455: ICMPQ1 = 1
456: ICMPQ2 = 0
457: NSUB = 0
458: *
459: DO 170 I = 1, N
460: IF( ABS( D( I ) ).LT.EPS ) THEN
461: D( I ) = SIGN( EPS, D( I ) )
462: END IF
463: 170 CONTINUE
464: *
465: DO 240 I = 1, NM1
466: IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
467: NSUB = NSUB + 1
468: IWORK( NSUB ) = ST
469: *
470: * Subproblem found. First determine its size and then
471: * apply divide and conquer on it.
472: *
473: IF( I.LT.NM1 ) THEN
474: *
475: * A subproblem with E(I) small for I < NM1.
476: *
477: NSIZE = I - ST + 1
478: IWORK( SIZEI+NSUB-1 ) = NSIZE
479: ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
480: *
481: * A subproblem with E(NM1) not too small but I = NM1.
482: *
483: NSIZE = N - ST + 1
484: IWORK( SIZEI+NSUB-1 ) = NSIZE
485: ELSE
486: *
487: * A subproblem with E(NM1) small. This implies an
488: * 1-by-1 subproblem at D(N), which is not solved
489: * explicitly.
490: *
491: NSIZE = I - ST + 1
492: IWORK( SIZEI+NSUB-1 ) = NSIZE
493: NSUB = NSUB + 1
494: IWORK( NSUB ) = N
495: IWORK( SIZEI+NSUB-1 ) = 1
496: CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
497: END IF
498: ST1 = ST - 1
499: IF( NSIZE.EQ.1 ) THEN
500: *
501: * This is a 1-by-1 subproblem and is not solved
502: * explicitly.
503: *
504: CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
505: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
506: *
507: * This is a small subproblem and is solved by DLASDQ.
508: *
509: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
510: $ RWORK( VT+ST1 ), N )
511: CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
512: $ RWORK( U+ST1 ), N )
513: CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
514: $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
515: $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
516: $ INFO )
517: IF( INFO.NE.0 ) THEN
518: RETURN
519: END IF
520: *
521: * In the real version, B is passed to DLASDQ and multiplied
1.9 bertrand 522: * internally by Q**H. Here B is complex and that product is
1.1 bertrand 523: * computed below in two steps (real and imaginary parts).
524: *
525: J = IRWB - 1
526: DO 190 JCOL = 1, NRHS
527: DO 180 JROW = ST, ST + NSIZE - 1
528: J = J + 1
529: RWORK( J ) = DBLE( B( JROW, JCOL ) )
530: 180 CONTINUE
531: 190 CONTINUE
532: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
533: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
534: $ ZERO, RWORK( IRWRB ), NSIZE )
535: J = IRWB - 1
536: DO 210 JCOL = 1, NRHS
537: DO 200 JROW = ST, ST + NSIZE - 1
538: J = J + 1
539: RWORK( J ) = DIMAG( B( JROW, JCOL ) )
540: 200 CONTINUE
541: 210 CONTINUE
542: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
543: $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
544: $ ZERO, RWORK( IRWIB ), NSIZE )
545: JREAL = IRWRB - 1
546: JIMAG = IRWIB - 1
547: DO 230 JCOL = 1, NRHS
548: DO 220 JROW = ST, ST + NSIZE - 1
549: JREAL = JREAL + 1
550: JIMAG = JIMAG + 1
551: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
552: $ RWORK( JIMAG ) )
553: 220 CONTINUE
554: 230 CONTINUE
555: *
556: CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
557: $ WORK( BX+ST1 ), N )
558: ELSE
559: *
560: * A large problem. Solve it using divide and conquer.
561: *
562: CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
563: $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
564: $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
565: $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
566: $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
567: $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
568: $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
569: $ RWORK( S+ST1 ), RWORK( NRWORK ),
570: $ IWORK( IWK ), INFO )
571: IF( INFO.NE.0 ) THEN
572: RETURN
573: END IF
574: BXST = BX + ST1
575: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
576: $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
577: $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
578: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
579: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
580: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
581: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
582: $ RWORK( C+ST1 ), RWORK( S+ST1 ),
583: $ RWORK( NRWORK ), IWORK( IWK ), INFO )
584: IF( INFO.NE.0 ) THEN
585: RETURN
586: END IF
587: END IF
588: ST = I + 1
589: END IF
590: 240 CONTINUE
591: *
592: * Apply the singular values and treat the tiny ones as zero.
593: *
594: TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
595: *
596: DO 250 I = 1, N
597: *
598: * Some of the elements in D can be negative because 1-by-1
599: * subproblems were not solved explicitly.
600: *
601: IF( ABS( D( I ) ).LE.TOL ) THEN
602: CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
603: ELSE
604: RANK = RANK + 1
605: CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
606: $ WORK( BX+I-1 ), N, INFO )
607: END IF
608: D( I ) = ABS( D( I ) )
609: 250 CONTINUE
610: *
611: * Now apply back the right singular vectors.
612: *
613: ICMPQ2 = 1
614: DO 320 I = 1, NSUB
615: ST = IWORK( I )
616: ST1 = ST - 1
617: NSIZE = IWORK( SIZEI+I-1 )
618: BXST = BX + ST1
619: IF( NSIZE.EQ.1 ) THEN
620: CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
621: ELSE IF( NSIZE.LE.SMLSIZ ) THEN
622: *
623: * Since B and BX are complex, the following call to DGEMM
624: * is performed in two steps (real and imaginary parts).
625: *
626: * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
627: * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
628: * $ B( ST, 1 ), LDB )
629: *
630: J = BXST - N - 1
631: JREAL = IRWB - 1
632: DO 270 JCOL = 1, NRHS
633: J = J + N
634: DO 260 JROW = 1, NSIZE
635: JREAL = JREAL + 1
636: RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
637: 260 CONTINUE
638: 270 CONTINUE
639: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
640: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
641: $ RWORK( IRWRB ), NSIZE )
642: J = BXST - N - 1
643: JIMAG = IRWB - 1
644: DO 290 JCOL = 1, NRHS
645: J = J + N
646: DO 280 JROW = 1, NSIZE
647: JIMAG = JIMAG + 1
648: RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
649: 280 CONTINUE
650: 290 CONTINUE
651: CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
652: $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
653: $ RWORK( IRWIB ), NSIZE )
654: JREAL = IRWRB - 1
655: JIMAG = IRWIB - 1
656: DO 310 JCOL = 1, NRHS
657: DO 300 JROW = ST, ST + NSIZE - 1
658: JREAL = JREAL + 1
659: JIMAG = JIMAG + 1
660: B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
661: $ RWORK( JIMAG ) )
662: 300 CONTINUE
663: 310 CONTINUE
664: ELSE
665: CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
666: $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
667: $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
668: $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
669: $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
670: $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
671: $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
672: $ RWORK( C+ST1 ), RWORK( S+ST1 ),
673: $ RWORK( NRWORK ), IWORK( IWK ), INFO )
674: IF( INFO.NE.0 ) THEN
675: RETURN
676: END IF
677: END IF
678: 320 CONTINUE
679: *
680: * Unscale and sort the singular values.
681: *
682: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
683: CALL DLASRT( 'D', N, D, INFO )
684: CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
685: *
686: RETURN
687: *
688: * End of ZLALSD
689: *
690: 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