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