Annotation of rpl/lapack/lapack/dlasd6.f, revision 1.17
1.14 bertrand 1: *> \brief \b DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc.
1.11 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLASD6 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd6.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd6.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd6.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
22: * IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
23: * LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
24: * IWORK, INFO )
25: *
26: * .. Scalar Arguments ..
27: * INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
28: * $ NR, SQRE
29: * DOUBLE PRECISION ALPHA, BETA, C, S
30: * ..
31: * .. Array Arguments ..
32: * INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
33: * $ PERM( * )
34: * DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ),
35: * $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
36: * $ VF( * ), VL( * ), WORK( * ), Z( * )
37: * ..
38: *
39: *
40: *> \par Purpose:
41: * =============
42: *>
43: *> \verbatim
44: *>
45: *> DLASD6 computes the SVD of an updated upper bidiagonal matrix B
46: *> obtained by merging two smaller ones by appending a row. This
47: *> routine is used only for the problem which requires all singular
48: *> values and optionally singular vector matrices in factored form.
49: *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
50: *> A related subroutine, DLASD1, handles the case in which all singular
51: *> values and singular vectors of the bidiagonal matrix are desired.
52: *>
53: *> DLASD6 computes the SVD as follows:
54: *>
55: *> ( D1(in) 0 0 0 )
56: *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
57: *> ( 0 0 D2(in) 0 )
58: *>
59: *> = U(out) * ( D(out) 0) * VT(out)
60: *>
61: *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
62: *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
63: *> elsewhere; and the entry b is empty if SQRE = 0.
64: *>
65: *> The singular values of B can be computed using D1, D2, the first
66: *> components of all the right singular vectors of the lower block, and
67: *> the last components of all the right singular vectors of the upper
68: *> block. These components are stored and updated in VF and VL,
69: *> respectively, in DLASD6. Hence U and VT are not explicitly
70: *> referenced.
71: *>
72: *> The singular values are stored in D. The algorithm consists of two
73: *> stages:
74: *>
75: *> The first stage consists of deflating the size of the problem
76: *> when there are multiple singular values or if there is a zero
77: *> in the Z vector. For each such occurence the dimension of the
78: *> secular equation problem is reduced by one. This stage is
79: *> performed by the routine DLASD7.
80: *>
81: *> The second stage consists of calculating the updated
82: *> singular values. This is done by finding the roots of the
83: *> secular equation via the routine DLASD4 (as called by DLASD8).
84: *> This routine also updates VF and VL and computes the distances
85: *> between the updated singular values and the old singular
86: *> values.
87: *>
88: *> DLASD6 is called from DLASDA.
89: *> \endverbatim
90: *
91: * Arguments:
92: * ==========
93: *
94: *> \param[in] ICOMPQ
95: *> \verbatim
96: *> ICOMPQ is INTEGER
97: *> Specifies whether singular vectors are to be computed in
98: *> factored form:
99: *> = 0: Compute singular values only.
100: *> = 1: Compute singular vectors in factored form as well.
101: *> \endverbatim
102: *>
103: *> \param[in] NL
104: *> \verbatim
105: *> NL is INTEGER
106: *> The row dimension of the upper block. NL >= 1.
107: *> \endverbatim
108: *>
109: *> \param[in] NR
110: *> \verbatim
111: *> NR is INTEGER
112: *> The row dimension of the lower block. NR >= 1.
113: *> \endverbatim
114: *>
115: *> \param[in] SQRE
116: *> \verbatim
117: *> SQRE is INTEGER
118: *> = 0: the lower block is an NR-by-NR square matrix.
119: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
120: *>
121: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
122: *> and column dimension M = N + SQRE.
123: *> \endverbatim
124: *>
125: *> \param[in,out] D
126: *> \verbatim
127: *> D is DOUBLE PRECISION array, dimension ( NL+NR+1 ).
128: *> On entry D(1:NL,1:NL) contains the singular values of the
129: *> upper block, and D(NL+2:N) contains the singular values
130: *> of the lower block. On exit D(1:N) contains the singular
131: *> values of the modified matrix.
132: *> \endverbatim
133: *>
134: *> \param[in,out] VF
135: *> \verbatim
136: *> VF is DOUBLE PRECISION array, dimension ( M )
137: *> On entry, VF(1:NL+1) contains the first components of all
138: *> right singular vectors of the upper block; and VF(NL+2:M)
139: *> contains the first components of all right singular vectors
140: *> of the lower block. On exit, VF contains the first components
141: *> of all right singular vectors of the bidiagonal matrix.
142: *> \endverbatim
143: *>
144: *> \param[in,out] VL
145: *> \verbatim
146: *> VL is DOUBLE PRECISION array, dimension ( M )
147: *> On entry, VL(1:NL+1) contains the last components of all
148: *> right singular vectors of the upper block; and VL(NL+2:M)
149: *> contains the last components of all right singular vectors of
150: *> the lower block. On exit, VL contains the last components of
151: *> all right singular vectors of the bidiagonal matrix.
152: *> \endverbatim
153: *>
154: *> \param[in,out] ALPHA
155: *> \verbatim
156: *> ALPHA is DOUBLE PRECISION
157: *> Contains the diagonal element associated with the added row.
158: *> \endverbatim
159: *>
160: *> \param[in,out] BETA
161: *> \verbatim
162: *> BETA is DOUBLE PRECISION
163: *> Contains the off-diagonal element associated with the added
164: *> row.
165: *> \endverbatim
166: *>
1.17 ! bertrand 167: *> \param[in,out] IDXQ
1.11 bertrand 168: *> \verbatim
169: *> IDXQ is INTEGER array, dimension ( N )
170: *> This contains the permutation which will reintegrate the
171: *> subproblem just solved back into sorted order, i.e.
172: *> D( IDXQ( I = 1, N ) ) will be in ascending order.
173: *> \endverbatim
174: *>
175: *> \param[out] PERM
176: *> \verbatim
177: *> PERM is INTEGER array, dimension ( N )
178: *> The permutations (from deflation and sorting) to be applied
179: *> to each block. Not referenced if ICOMPQ = 0.
180: *> \endverbatim
181: *>
182: *> \param[out] GIVPTR
183: *> \verbatim
184: *> GIVPTR is INTEGER
185: *> The number of Givens rotations which took place in this
186: *> subproblem. Not referenced if ICOMPQ = 0.
187: *> \endverbatim
188: *>
189: *> \param[out] GIVCOL
190: *> \verbatim
191: *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
192: *> Each pair of numbers indicates a pair of columns to take place
193: *> in a Givens rotation. Not referenced if ICOMPQ = 0.
194: *> \endverbatim
195: *>
196: *> \param[in] LDGCOL
197: *> \verbatim
198: *> LDGCOL is INTEGER
199: *> leading dimension of GIVCOL, must be at least N.
200: *> \endverbatim
201: *>
202: *> \param[out] GIVNUM
203: *> \verbatim
204: *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
205: *> Each number indicates the C or S value to be used in the
206: *> corresponding Givens rotation. Not referenced if ICOMPQ = 0.
207: *> \endverbatim
208: *>
209: *> \param[in] LDGNUM
210: *> \verbatim
211: *> LDGNUM is INTEGER
212: *> The leading dimension of GIVNUM and POLES, must be at least N.
213: *> \endverbatim
214: *>
215: *> \param[out] POLES
216: *> \verbatim
217: *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
218: *> On exit, POLES(1,*) is an array containing the new singular
219: *> values obtained from solving the secular equation, and
220: *> POLES(2,*) is an array containing the poles in the secular
221: *> equation. Not referenced if ICOMPQ = 0.
222: *> \endverbatim
223: *>
224: *> \param[out] DIFL
225: *> \verbatim
226: *> DIFL is DOUBLE PRECISION array, dimension ( N )
227: *> On exit, DIFL(I) is the distance between I-th updated
228: *> (undeflated) singular value and the I-th (undeflated) old
229: *> singular value.
230: *> \endverbatim
231: *>
232: *> \param[out] DIFR
233: *> \verbatim
234: *> DIFR is DOUBLE PRECISION array,
235: *> dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and
236: *> dimension ( N ) if ICOMPQ = 0.
237: *> On exit, DIFR(I, 1) is the distance between I-th updated
238: *> (undeflated) singular value and the I+1-th (undeflated) old
239: *> singular value.
240: *>
241: *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
242: *> normalizing factors for the right singular vector matrix.
243: *>
244: *> See DLASD8 for details on DIFL and DIFR.
245: *> \endverbatim
246: *>
247: *> \param[out] Z
248: *> \verbatim
249: *> Z is DOUBLE PRECISION array, dimension ( M )
250: *> The first elements of this array contain the components
251: *> of the deflation-adjusted updating row vector.
252: *> \endverbatim
253: *>
254: *> \param[out] K
255: *> \verbatim
256: *> K is INTEGER
257: *> Contains the dimension of the non-deflated matrix,
258: *> This is the order of the related secular equation. 1 <= K <=N.
259: *> \endverbatim
260: *>
261: *> \param[out] C
262: *> \verbatim
263: *> C is DOUBLE PRECISION
264: *> C contains garbage if SQRE =0 and the C-value of a Givens
265: *> rotation related to the right null space if SQRE = 1.
266: *> \endverbatim
267: *>
268: *> \param[out] S
269: *> \verbatim
270: *> S is DOUBLE PRECISION
271: *> S contains garbage if SQRE =0 and the S-value of a Givens
272: *> rotation related to the right null space if SQRE = 1.
273: *> \endverbatim
274: *>
275: *> \param[out] WORK
276: *> \verbatim
277: *> WORK is DOUBLE PRECISION array, dimension ( 4 * M )
278: *> \endverbatim
279: *>
280: *> \param[out] IWORK
281: *> \verbatim
282: *> IWORK is INTEGER array, dimension ( 3 * N )
283: *> \endverbatim
284: *>
285: *> \param[out] INFO
286: *> \verbatim
287: *> INFO is INTEGER
288: *> = 0: successful exit.
289: *> < 0: if INFO = -i, the i-th argument had an illegal value.
290: *> > 0: if INFO = 1, a singular value did not converge
291: *> \endverbatim
292: *
293: * Authors:
294: * ========
295: *
296: *> \author Univ. of Tennessee
297: *> \author Univ. of California Berkeley
298: *> \author Univ. of Colorado Denver
299: *> \author NAG Ltd.
300: *
1.17 ! bertrand 301: *> \date November 2015
1.11 bertrand 302: *
303: *> \ingroup auxOTHERauxiliary
304: *
305: *> \par Contributors:
306: * ==================
307: *>
308: *> Ming Gu and Huan Ren, Computer Science Division, University of
309: *> California at Berkeley, USA
310: *>
311: * =====================================================================
1.1 bertrand 312: SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
313: $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
314: $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
315: $ IWORK, INFO )
316: *
1.17 ! bertrand 317: * -- LAPACK auxiliary routine (version 3.6.0) --
1.1 bertrand 318: * -- LAPACK is a software package provided by Univ. of Tennessee, --
319: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.17 ! bertrand 320: * November 2015
1.1 bertrand 321: *
322: * .. Scalar Arguments ..
323: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
324: $ NR, SQRE
325: DOUBLE PRECISION ALPHA, BETA, C, S
326: * ..
327: * .. Array Arguments ..
328: INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
329: $ PERM( * )
330: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ),
331: $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
332: $ VF( * ), VL( * ), WORK( * ), Z( * )
333: * ..
334: *
335: * =====================================================================
336: *
337: * .. Parameters ..
338: DOUBLE PRECISION ONE, ZERO
339: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
340: * ..
341: * .. Local Scalars ..
342: INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
343: $ N, N1, N2
344: DOUBLE PRECISION ORGNRM
345: * ..
346: * .. External Subroutines ..
347: EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA
348: * ..
349: * .. Intrinsic Functions ..
350: INTRINSIC ABS, MAX
351: * ..
352: * .. Executable Statements ..
353: *
354: * Test the input parameters.
355: *
356: INFO = 0
357: N = NL + NR + 1
358: M = N + SQRE
359: *
360: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
361: INFO = -1
362: ELSE IF( NL.LT.1 ) THEN
363: INFO = -2
364: ELSE IF( NR.LT.1 ) THEN
365: INFO = -3
366: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
367: INFO = -4
368: ELSE IF( LDGCOL.LT.N ) THEN
369: INFO = -14
370: ELSE IF( LDGNUM.LT.N ) THEN
371: INFO = -16
372: END IF
373: IF( INFO.NE.0 ) THEN
374: CALL XERBLA( 'DLASD6', -INFO )
375: RETURN
376: END IF
377: *
378: * The following values are for bookkeeping purposes only. They are
379: * integer pointers which indicate the portion of the workspace
380: * used by a particular array in DLASD7 and DLASD8.
381: *
382: ISIGMA = 1
383: IW = ISIGMA + N
384: IVFW = IW + M
385: IVLW = IVFW + M
386: *
387: IDX = 1
388: IDXC = IDX + N
389: IDXP = IDXC + N
390: *
391: * Scale.
392: *
393: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
394: D( NL+1 ) = ZERO
395: DO 10 I = 1, N
396: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
397: ORGNRM = ABS( D( I ) )
398: END IF
399: 10 CONTINUE
400: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
401: ALPHA = ALPHA / ORGNRM
402: BETA = BETA / ORGNRM
403: *
404: * Sort and Deflate singular values.
405: *
406: CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
407: $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
408: $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
409: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
410: $ INFO )
411: *
412: * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
413: *
414: CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
415: $ WORK( ISIGMA ), WORK( IW ), INFO )
416: *
1.17 ! bertrand 417: * Report the possible convergence failure.
1.8 bertrand 418: *
419: IF( INFO.NE.0 ) THEN
420: RETURN
421: END IF
422: *
1.1 bertrand 423: * Save the poles if ICOMPQ = 1.
424: *
425: IF( ICOMPQ.EQ.1 ) THEN
426: CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 )
427: CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
428: END IF
429: *
430: * Unscale.
431: *
432: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
433: *
434: * Prepare the IDXQ sorting permutation.
435: *
436: N1 = K
437: N2 = N - K
438: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
439: *
440: RETURN
441: *
442: * End of DLASD6
443: *
444: END
CVSweb interface <joel.bertrand@systella.fr>