1: *> \brief \b DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagonal d and off-diagonal e. Used by sbdsdc.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASDA + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasda.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasda.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasda.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
22: * DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
23: * PERM, GIVNUM, C, S, WORK, IWORK, INFO )
24: *
25: * .. Scalar Arguments ..
26: * INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
27: * ..
28: * .. Array Arguments ..
29: * INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
30: * $ K( * ), PERM( LDGCOL, * )
31: * DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
32: * $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
33: * $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
34: * $ Z( LDU, * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> Using a divide and conquer approach, DLASDA computes the singular
44: *> value decomposition (SVD) of a real upper bidiagonal N-by-M matrix
45: *> B with diagonal D and offdiagonal E, where M = N + SQRE. The
46: *> algorithm computes the singular values in the SVD B = U * S * VT.
47: *> The orthogonal matrices U and VT are optionally computed in
48: *> compact form.
49: *>
50: *> A related subroutine, DLASD0, computes the singular values and
51: *> the singular vectors in explicit form.
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] ICOMPQ
58: *> \verbatim
59: *> ICOMPQ is INTEGER
60: *> Specifies whether singular vectors are to be computed
61: *> in compact form, as follows
62: *> = 0: Compute singular values only.
63: *> = 1: Compute singular vectors of upper bidiagonal
64: *> matrix in compact form.
65: *> \endverbatim
66: *>
67: *> \param[in] SMLSIZ
68: *> \verbatim
69: *> SMLSIZ is INTEGER
70: *> The maximum size of the subproblems at the bottom of the
71: *> computation tree.
72: *> \endverbatim
73: *>
74: *> \param[in] N
75: *> \verbatim
76: *> N is INTEGER
77: *> The row dimension of the upper bidiagonal matrix. This is
78: *> also the dimension of the main diagonal array D.
79: *> \endverbatim
80: *>
81: *> \param[in] SQRE
82: *> \verbatim
83: *> SQRE is INTEGER
84: *> Specifies the column dimension of the bidiagonal matrix.
85: *> = 0: The bidiagonal matrix has column dimension M = N;
86: *> = 1: The bidiagonal matrix has column dimension M = N + 1.
87: *> \endverbatim
88: *>
89: *> \param[in,out] D
90: *> \verbatim
91: *> D is DOUBLE PRECISION array, dimension ( N )
92: *> On entry D contains the main diagonal of the bidiagonal
93: *> matrix. On exit D, if INFO = 0, contains its singular values.
94: *> \endverbatim
95: *>
96: *> \param[in] E
97: *> \verbatim
98: *> E is DOUBLE PRECISION array, dimension ( M-1 )
99: *> Contains the subdiagonal entries of the bidiagonal matrix.
100: *> On exit, E has been destroyed.
101: *> \endverbatim
102: *>
103: *> \param[out] U
104: *> \verbatim
105: *> U is DOUBLE PRECISION array,
106: *> dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced
107: *> if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left
108: *> singular vector matrices of all subproblems at the bottom
109: *> level.
110: *> \endverbatim
111: *>
112: *> \param[in] LDU
113: *> \verbatim
114: *> LDU is INTEGER, LDU = > N.
115: *> The leading dimension of arrays U, VT, DIFL, DIFR, POLES,
116: *> GIVNUM, and Z.
117: *> \endverbatim
118: *>
119: *> \param[out] VT
120: *> \verbatim
121: *> VT is DOUBLE PRECISION array,
122: *> dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced
123: *> if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right
124: *> singular vector matrices of all subproblems at the bottom
125: *> level.
126: *> \endverbatim
127: *>
128: *> \param[out] K
129: *> \verbatim
130: *> K is INTEGER array,
131: *> dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0.
132: *> If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th
133: *> secular equation on the computation tree.
134: *> \endverbatim
135: *>
136: *> \param[out] DIFL
137: *> \verbatim
138: *> DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ),
139: *> where NLVL = floor(log_2 (N/SMLSIZ))).
140: *> \endverbatim
141: *>
142: *> \param[out] DIFR
143: *> \verbatim
144: *> DIFR is DOUBLE PRECISION array,
145: *> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and
146: *> dimension ( N ) if ICOMPQ = 0.
147: *> If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1)
148: *> record distances between singular values on the I-th
149: *> level and singular values on the (I -1)-th level, and
150: *> DIFR(1:N, 2 * I ) contains the normalizing factors for
151: *> the right singular vector matrix. See DLASD8 for details.
152: *> \endverbatim
153: *>
154: *> \param[out] Z
155: *> \verbatim
156: *> Z is DOUBLE PRECISION array,
157: *> dimension ( LDU, NLVL ) if ICOMPQ = 1 and
158: *> dimension ( N ) if ICOMPQ = 0.
159: *> The first K elements of Z(1, I) contain the components of
160: *> the deflation-adjusted updating row vector for subproblems
161: *> on the I-th level.
162: *> \endverbatim
163: *>
164: *> \param[out] POLES
165: *> \verbatim
166: *> POLES is DOUBLE PRECISION array,
167: *> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced
168: *> if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and
169: *> POLES(1, 2*I) contain the new and old singular values
170: *> involved in the secular equations on the I-th level.
171: *> \endverbatim
172: *>
173: *> \param[out] GIVPTR
174: *> \verbatim
175: *> GIVPTR is INTEGER array,
176: *> dimension ( N ) if ICOMPQ = 1, and not referenced if
177: *> ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records
178: *> the number of Givens rotations performed on the I-th
179: *> problem on the computation tree.
180: *> \endverbatim
181: *>
182: *> \param[out] GIVCOL
183: *> \verbatim
184: *> GIVCOL is INTEGER array,
185: *> dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not
186: *> referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
187: *> GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations
188: *> of Givens rotations performed on the I-th level on the
189: *> computation tree.
190: *> \endverbatim
191: *>
192: *> \param[in] LDGCOL
193: *> \verbatim
194: *> LDGCOL is INTEGER, LDGCOL = > N.
195: *> The leading dimension of arrays GIVCOL and PERM.
196: *> \endverbatim
197: *>
198: *> \param[out] PERM
199: *> \verbatim
200: *> PERM is INTEGER array,
201: *> dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced
202: *> if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records
203: *> permutations done on the I-th level of the computation tree.
204: *> \endverbatim
205: *>
206: *> \param[out] GIVNUM
207: *> \verbatim
208: *> GIVNUM is DOUBLE PRECISION array,
209: *> dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not
210: *> referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I,
211: *> GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S-
212: *> values of Givens rotations performed on the I-th level on
213: *> the computation tree.
214: *> \endverbatim
215: *>
216: *> \param[out] C
217: *> \verbatim
218: *> C is DOUBLE PRECISION array,
219: *> dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0.
220: *> If ICOMPQ = 1 and the I-th subproblem is not square, on exit,
221: *> C( I ) contains the C-value of a Givens rotation related to
222: *> the right null space of the I-th subproblem.
223: *> \endverbatim
224: *>
225: *> \param[out] S
226: *> \verbatim
227: *> S is DOUBLE PRECISION array, dimension ( N ) if
228: *> ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1
229: *> and the I-th subproblem is not square, on exit, S( I )
230: *> contains the S-value of a Givens rotation related to
231: *> the right null space of the I-th subproblem.
232: *> \endverbatim
233: *>
234: *> \param[out] WORK
235: *> \verbatim
236: *> WORK is DOUBLE PRECISION array, dimension
237: *> (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)).
238: *> \endverbatim
239: *>
240: *> \param[out] IWORK
241: *> \verbatim
242: *> IWORK is INTEGER array, dimension (7*N)
243: *> \endverbatim
244: *>
245: *> \param[out] INFO
246: *> \verbatim
247: *> INFO is INTEGER
248: *> = 0: successful exit.
249: *> < 0: if INFO = -i, the i-th argument had an illegal value.
250: *> > 0: if INFO = 1, a singular value did not converge
251: *> \endverbatim
252: *
253: * Authors:
254: * ========
255: *
256: *> \author Univ. of Tennessee
257: *> \author Univ. of California Berkeley
258: *> \author Univ. of Colorado Denver
259: *> \author NAG Ltd.
260: *
261: *> \ingroup OTHERauxiliary
262: *
263: *> \par Contributors:
264: * ==================
265: *>
266: *> Ming Gu and Huan Ren, Computer Science Division, University of
267: *> California at Berkeley, USA
268: *>
269: * =====================================================================
270: SUBROUTINE DLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
271: $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
272: $ PERM, GIVNUM, C, S, WORK, IWORK, INFO )
273: *
274: * -- LAPACK auxiliary routine --
275: * -- LAPACK is a software package provided by Univ. of Tennessee, --
276: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
277: *
278: * .. Scalar Arguments ..
279: INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
280: * ..
281: * .. Array Arguments ..
282: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
283: $ K( * ), PERM( LDGCOL, * )
284: DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
285: $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
286: $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
287: $ Z( LDU, * )
288: * ..
289: *
290: * =====================================================================
291: *
292: * .. Parameters ..
293: DOUBLE PRECISION ZERO, ONE
294: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
295: * ..
296: * .. Local Scalars ..
297: INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
298: $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
299: $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU,
300: $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI
301: DOUBLE PRECISION ALPHA, BETA
302: * ..
303: * .. External Subroutines ..
304: EXTERNAL DCOPY, DLASD6, DLASDQ, DLASDT, DLASET, XERBLA
305: * ..
306: * .. Executable Statements ..
307: *
308: * Test the input parameters.
309: *
310: INFO = 0
311: *
312: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
313: INFO = -1
314: ELSE IF( SMLSIZ.LT.3 ) THEN
315: INFO = -2
316: ELSE IF( N.LT.0 ) THEN
317: INFO = -3
318: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
319: INFO = -4
320: ELSE IF( LDU.LT.( N+SQRE ) ) THEN
321: INFO = -8
322: ELSE IF( LDGCOL.LT.N ) THEN
323: INFO = -17
324: END IF
325: IF( INFO.NE.0 ) THEN
326: CALL XERBLA( 'DLASDA', -INFO )
327: RETURN
328: END IF
329: *
330: M = N + SQRE
331: *
332: * If the input matrix is too small, call DLASDQ to find the SVD.
333: *
334: IF( N.LE.SMLSIZ ) THEN
335: IF( ICOMPQ.EQ.0 ) THEN
336: CALL DLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU,
337: $ U, LDU, WORK, INFO )
338: ELSE
339: CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU,
340: $ U, LDU, WORK, INFO )
341: END IF
342: RETURN
343: END IF
344: *
345: * Book-keeping and set up the computation tree.
346: *
347: INODE = 1
348: NDIML = INODE + N
349: NDIMR = NDIML + N
350: IDXQ = NDIMR + N
351: IWK = IDXQ + N
352: *
353: NCC = 0
354: NRU = 0
355: *
356: SMLSZP = SMLSIZ + 1
357: VF = 1
358: VL = VF + M
359: NWORK1 = VL + M
360: NWORK2 = NWORK1 + SMLSZP*SMLSZP
361: *
362: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
363: $ IWORK( NDIMR ), SMLSIZ )
364: *
365: * for the nodes on bottom level of the tree, solve
366: * their subproblems by DLASDQ.
367: *
368: NDB1 = ( ND+1 ) / 2
369: DO 30 I = NDB1, ND
370: *
371: * IC : center row of each node
372: * NL : number of rows of left subproblem
373: * NR : number of rows of right subproblem
374: * NLF: starting row of the left subproblem
375: * NRF: starting row of the right subproblem
376: *
377: I1 = I - 1
378: IC = IWORK( INODE+I1 )
379: NL = IWORK( NDIML+I1 )
380: NLP1 = NL + 1
381: NR = IWORK( NDIMR+I1 )
382: NLF = IC - NL
383: NRF = IC + 1
384: IDXQI = IDXQ + NLF - 2
385: VFI = VF + NLF - 1
386: VLI = VL + NLF - 1
387: SQREI = 1
388: IF( ICOMPQ.EQ.0 ) THEN
389: CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ),
390: $ SMLSZP )
391: CALL DLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ),
392: $ E( NLF ), WORK( NWORK1 ), SMLSZP,
393: $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL,
394: $ WORK( NWORK2 ), INFO )
395: ITEMP = NWORK1 + NL*SMLSZP
396: CALL DCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
397: CALL DCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
398: ELSE
399: CALL DLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU )
400: CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU )
401: CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ),
402: $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU,
403: $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO )
404: CALL DCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 )
405: CALL DCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 )
406: END IF
407: IF( INFO.NE.0 ) THEN
408: RETURN
409: END IF
410: DO 10 J = 1, NL
411: IWORK( IDXQI+J ) = J
412: 10 CONTINUE
413: IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN
414: SQREI = 0
415: ELSE
416: SQREI = 1
417: END IF
418: IDXQI = IDXQI + NLP1
419: VFI = VFI + NLP1
420: VLI = VLI + NLP1
421: NRP1 = NR + SQREI
422: IF( ICOMPQ.EQ.0 ) THEN
423: CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ),
424: $ SMLSZP )
425: CALL DLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ),
426: $ E( NRF ), WORK( NWORK1 ), SMLSZP,
427: $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR,
428: $ WORK( NWORK2 ), INFO )
429: ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP
430: CALL DCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
431: CALL DCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
432: ELSE
433: CALL DLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU )
434: CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU )
435: CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ),
436: $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU,
437: $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO )
438: CALL DCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 )
439: CALL DCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 )
440: END IF
441: IF( INFO.NE.0 ) THEN
442: RETURN
443: END IF
444: DO 20 J = 1, NR
445: IWORK( IDXQI+J ) = J
446: 20 CONTINUE
447: 30 CONTINUE
448: *
449: * Now conquer each subproblem bottom-up.
450: *
451: J = 2**NLVL
452: DO 50 LVL = NLVL, 1, -1
453: LVL2 = LVL*2 - 1
454: *
455: * Find the first node LF and last node LL on
456: * the current level LVL.
457: *
458: IF( LVL.EQ.1 ) THEN
459: LF = 1
460: LL = 1
461: ELSE
462: LF = 2**( LVL-1 )
463: LL = 2*LF - 1
464: END IF
465: DO 40 I = LF, LL
466: IM1 = I - 1
467: IC = IWORK( INODE+IM1 )
468: NL = IWORK( NDIML+IM1 )
469: NR = IWORK( NDIMR+IM1 )
470: NLF = IC - NL
471: NRF = IC + 1
472: IF( I.EQ.LL ) THEN
473: SQREI = SQRE
474: ELSE
475: SQREI = 1
476: END IF
477: VFI = VF + NLF - 1
478: VLI = VL + NLF - 1
479: IDXQI = IDXQ + NLF - 1
480: ALPHA = D( IC )
481: BETA = E( IC )
482: IF( ICOMPQ.EQ.0 ) THEN
483: CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
484: $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
485: $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL,
486: $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z,
487: $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ),
488: $ IWORK( IWK ), INFO )
489: ELSE
490: J = J - 1
491: CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
492: $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
493: $ IWORK( IDXQI ), PERM( NLF, LVL ),
494: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
495: $ GIVNUM( NLF, LVL2 ), LDU,
496: $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ),
497: $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ),
498: $ C( J ), S( J ), WORK( NWORK1 ),
499: $ IWORK( IWK ), INFO )
500: END IF
501: IF( INFO.NE.0 ) THEN
502: RETURN
503: END IF
504: 40 CONTINUE
505: 50 CONTINUE
506: *
507: RETURN
508: *
509: * End of DLASDA
510: *
511: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasda.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack