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.
243: *> Dimension must be at least (7 * N).
244: *> \endverbatim
245: *>
246: *> \param[out] INFO
247: *> \verbatim
248: *> INFO is INTEGER
249: *> = 0: successful exit.
250: *> < 0: if INFO = -i, the i-th argument had an illegal value.
251: *> > 0: if INFO = 1, a singular value did not converge
252: *> \endverbatim
253: *
254: * Authors:
255: * ========
256: *
257: *> \author Univ. of Tennessee
258: *> \author Univ. of California Berkeley
259: *> \author Univ. of Colorado Denver
260: *> \author NAG Ltd.
261: *
262: *> \date September 2012
263: *
264: *> \ingroup auxOTHERauxiliary
265: *
266: *> \par Contributors:
267: * ==================
268: *>
269: *> Ming Gu and Huan Ren, Computer Science Division, University of
270: *> California at Berkeley, USA
271: *>
272: * =====================================================================
273: SUBROUTINE DLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
274: $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
275: $ PERM, GIVNUM, C, S, WORK, IWORK, INFO )
276: *
277: * -- LAPACK auxiliary routine (version 3.4.2) --
278: * -- LAPACK is a software package provided by Univ. of Tennessee, --
279: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280: * September 2012
281: *
282: * .. Scalar Arguments ..
283: INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
284: * ..
285: * .. Array Arguments ..
286: INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
287: $ K( * ), PERM( LDGCOL, * )
288: DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
289: $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
290: $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ),
291: $ Z( LDU, * )
292: * ..
293: *
294: * =====================================================================
295: *
296: * .. Parameters ..
297: DOUBLE PRECISION ZERO, ONE
298: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
299: * ..
300: * .. Local Scalars ..
301: INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
302: $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
303: $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU,
304: $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI
305: DOUBLE PRECISION ALPHA, BETA
306: * ..
307: * .. External Subroutines ..
308: EXTERNAL DCOPY, DLASD6, DLASDQ, DLASDT, DLASET, XERBLA
309: * ..
310: * .. Executable Statements ..
311: *
312: * Test the input parameters.
313: *
314: INFO = 0
315: *
316: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
317: INFO = -1
318: ELSE IF( SMLSIZ.LT.3 ) THEN
319: INFO = -2
320: ELSE IF( N.LT.0 ) THEN
321: INFO = -3
322: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
323: INFO = -4
324: ELSE IF( LDU.LT.( N+SQRE ) ) THEN
325: INFO = -8
326: ELSE IF( LDGCOL.LT.N ) THEN
327: INFO = -17
328: END IF
329: IF( INFO.NE.0 ) THEN
330: CALL XERBLA( 'DLASDA', -INFO )
331: RETURN
332: END IF
333: *
334: M = N + SQRE
335: *
336: * If the input matrix is too small, call DLASDQ to find the SVD.
337: *
338: IF( N.LE.SMLSIZ ) THEN
339: IF( ICOMPQ.EQ.0 ) THEN
340: CALL DLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU,
341: $ U, LDU, WORK, INFO )
342: ELSE
343: CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU,
344: $ U, LDU, WORK, INFO )
345: END IF
346: RETURN
347: END IF
348: *
349: * Book-keeping and set up the computation tree.
350: *
351: INODE = 1
352: NDIML = INODE + N
353: NDIMR = NDIML + N
354: IDXQ = NDIMR + N
355: IWK = IDXQ + N
356: *
357: NCC = 0
358: NRU = 0
359: *
360: SMLSZP = SMLSIZ + 1
361: VF = 1
362: VL = VF + M
363: NWORK1 = VL + M
364: NWORK2 = NWORK1 + SMLSZP*SMLSZP
365: *
366: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
367: $ IWORK( NDIMR ), SMLSIZ )
368: *
369: * for the nodes on bottom level of the tree, solve
370: * their subproblems by DLASDQ.
371: *
372: NDB1 = ( ND+1 ) / 2
373: DO 30 I = NDB1, ND
374: *
375: * IC : center row of each node
376: * NL : number of rows of left subproblem
377: * NR : number of rows of right subproblem
378: * NLF: starting row of the left subproblem
379: * NRF: starting row of the right subproblem
380: *
381: I1 = I - 1
382: IC = IWORK( INODE+I1 )
383: NL = IWORK( NDIML+I1 )
384: NLP1 = NL + 1
385: NR = IWORK( NDIMR+I1 )
386: NLF = IC - NL
387: NRF = IC + 1
388: IDXQI = IDXQ + NLF - 2
389: VFI = VF + NLF - 1
390: VLI = VL + NLF - 1
391: SQREI = 1
392: IF( ICOMPQ.EQ.0 ) THEN
393: CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ),
394: $ SMLSZP )
395: CALL DLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ),
396: $ E( NLF ), WORK( NWORK1 ), SMLSZP,
397: $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL,
398: $ WORK( NWORK2 ), INFO )
399: ITEMP = NWORK1 + NL*SMLSZP
400: CALL DCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
401: CALL DCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
402: ELSE
403: CALL DLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU )
404: CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU )
405: CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ),
406: $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU,
407: $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO )
408: CALL DCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 )
409: CALL DCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 )
410: END IF
411: IF( INFO.NE.0 ) THEN
412: RETURN
413: END IF
414: DO 10 J = 1, NL
415: IWORK( IDXQI+J ) = J
416: 10 CONTINUE
417: IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN
418: SQREI = 0
419: ELSE
420: SQREI = 1
421: END IF
422: IDXQI = IDXQI + NLP1
423: VFI = VFI + NLP1
424: VLI = VLI + NLP1
425: NRP1 = NR + SQREI
426: IF( ICOMPQ.EQ.0 ) THEN
427: CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ),
428: $ SMLSZP )
429: CALL DLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ),
430: $ E( NRF ), WORK( NWORK1 ), SMLSZP,
431: $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR,
432: $ WORK( NWORK2 ), INFO )
433: ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP
434: CALL DCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 )
435: CALL DCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 )
436: ELSE
437: CALL DLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU )
438: CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU )
439: CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ),
440: $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU,
441: $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO )
442: CALL DCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 )
443: CALL DCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 )
444: END IF
445: IF( INFO.NE.0 ) THEN
446: RETURN
447: END IF
448: DO 20 J = 1, NR
449: IWORK( IDXQI+J ) = J
450: 20 CONTINUE
451: 30 CONTINUE
452: *
453: * Now conquer each subproblem bottom-up.
454: *
455: J = 2**NLVL
456: DO 50 LVL = NLVL, 1, -1
457: LVL2 = LVL*2 - 1
458: *
459: * Find the first node LF and last node LL on
460: * the current level LVL.
461: *
462: IF( LVL.EQ.1 ) THEN
463: LF = 1
464: LL = 1
465: ELSE
466: LF = 2**( LVL-1 )
467: LL = 2*LF - 1
468: END IF
469: DO 40 I = LF, LL
470: IM1 = I - 1
471: IC = IWORK( INODE+IM1 )
472: NL = IWORK( NDIML+IM1 )
473: NR = IWORK( NDIMR+IM1 )
474: NLF = IC - NL
475: NRF = IC + 1
476: IF( I.EQ.LL ) THEN
477: SQREI = SQRE
478: ELSE
479: SQREI = 1
480: END IF
481: VFI = VF + NLF - 1
482: VLI = VL + NLF - 1
483: IDXQI = IDXQ + NLF - 1
484: ALPHA = D( IC )
485: BETA = E( IC )
486: IF( ICOMPQ.EQ.0 ) THEN
487: CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
488: $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
489: $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL,
490: $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z,
491: $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ),
492: $ IWORK( IWK ), INFO )
493: ELSE
494: J = J - 1
495: CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ),
496: $ WORK( VFI ), WORK( VLI ), ALPHA, BETA,
497: $ IWORK( IDXQI ), PERM( NLF, LVL ),
498: $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
499: $ GIVNUM( NLF, LVL2 ), LDU,
500: $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ),
501: $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ),
502: $ C( J ), S( J ), WORK( NWORK1 ),
503: $ IWORK( IWK ), INFO )
504: END IF
505: IF( INFO.NE.0 ) THEN
506: RETURN
507: END IF
508: 40 CONTINUE
509: 50 CONTINUE
510: *
511: RETURN
512: *
513: * End of DLASDA
514: *
515: 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