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