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.
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 occurrence 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: *>
167: *> \param[in,out] IDXQ
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 ( LDDIFR, 2 ) if ICOMPQ = 1 and
236: *> dimension ( K ) if ICOMPQ = 0.
237: *> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
238: *> defined and will not be referenced.
239: *>
240: *> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
241: *> normalizing factors for the right singular vector matrix.
242: *>
243: *> See DLASD8 for details on DIFL and DIFR.
244: *> \endverbatim
245: *>
246: *> \param[out] Z
247: *> \verbatim
248: *> Z is DOUBLE PRECISION array, dimension ( M )
249: *> The first elements of this array contain the components
250: *> of the deflation-adjusted updating row vector.
251: *> \endverbatim
252: *>
253: *> \param[out] K
254: *> \verbatim
255: *> K is INTEGER
256: *> Contains the dimension of the non-deflated matrix,
257: *> This is the order of the related secular equation. 1 <= K <=N.
258: *> \endverbatim
259: *>
260: *> \param[out] C
261: *> \verbatim
262: *> C is DOUBLE PRECISION
263: *> C contains garbage if SQRE =0 and the C-value of a Givens
264: *> rotation related to the right null space if SQRE = 1.
265: *> \endverbatim
266: *>
267: *> \param[out] S
268: *> \verbatim
269: *> S is DOUBLE PRECISION
270: *> S contains garbage if SQRE =0 and the S-value of a Givens
271: *> rotation related to the right null space if SQRE = 1.
272: *> \endverbatim
273: *>
274: *> \param[out] WORK
275: *> \verbatim
276: *> WORK is DOUBLE PRECISION array, dimension ( 4 * M )
277: *> \endverbatim
278: *>
279: *> \param[out] IWORK
280: *> \verbatim
281: *> IWORK is INTEGER array, dimension ( 3 * N )
282: *> \endverbatim
283: *>
284: *> \param[out] INFO
285: *> \verbatim
286: *> INFO is INTEGER
287: *> = 0: successful exit.
288: *> < 0: if INFO = -i, the i-th argument had an illegal value.
289: *> > 0: if INFO = 1, a singular value did not converge
290: *> \endverbatim
291: *
292: * Authors:
293: * ========
294: *
295: *> \author Univ. of Tennessee
296: *> \author Univ. of California Berkeley
297: *> \author Univ. of Colorado Denver
298: *> \author NAG Ltd.
299: *
300: *> \date June 2016
301: *
302: *> \ingroup auxOTHERauxiliary
303: *
304: *> \par Contributors:
305: * ==================
306: *>
307: *> Ming Gu and Huan Ren, Computer Science Division, University of
308: *> California at Berkeley, USA
309: *>
310: * =====================================================================
311: SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
312: $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
313: $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
314: $ IWORK, INFO )
315: *
316: * -- LAPACK auxiliary routine (version 3.6.1) --
317: * -- LAPACK is a software package provided by Univ. of Tennessee, --
318: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
319: * June 2016
320: *
321: * .. Scalar Arguments ..
322: INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
323: $ NR, SQRE
324: DOUBLE PRECISION ALPHA, BETA, C, S
325: * ..
326: * .. Array Arguments ..
327: INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
328: $ PERM( * )
329: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ),
330: $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
331: $ VF( * ), VL( * ), WORK( * ), Z( * )
332: * ..
333: *
334: * =====================================================================
335: *
336: * .. Parameters ..
337: DOUBLE PRECISION ONE, ZERO
338: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
339: * ..
340: * .. Local Scalars ..
341: INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
342: $ N, N1, N2
343: DOUBLE PRECISION ORGNRM
344: * ..
345: * .. External Subroutines ..
346: EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA
347: * ..
348: * .. Intrinsic Functions ..
349: INTRINSIC ABS, MAX
350: * ..
351: * .. Executable Statements ..
352: *
353: * Test the input parameters.
354: *
355: INFO = 0
356: N = NL + NR + 1
357: M = N + SQRE
358: *
359: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
360: INFO = -1
361: ELSE IF( NL.LT.1 ) THEN
362: INFO = -2
363: ELSE IF( NR.LT.1 ) THEN
364: INFO = -3
365: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
366: INFO = -4
367: ELSE IF( LDGCOL.LT.N ) THEN
368: INFO = -14
369: ELSE IF( LDGNUM.LT.N ) THEN
370: INFO = -16
371: END IF
372: IF( INFO.NE.0 ) THEN
373: CALL XERBLA( 'DLASD6', -INFO )
374: RETURN
375: END IF
376: *
377: * The following values are for bookkeeping purposes only. They are
378: * integer pointers which indicate the portion of the workspace
379: * used by a particular array in DLASD7 and DLASD8.
380: *
381: ISIGMA = 1
382: IW = ISIGMA + N
383: IVFW = IW + M
384: IVLW = IVFW + M
385: *
386: IDX = 1
387: IDXC = IDX + N
388: IDXP = IDXC + N
389: *
390: * Scale.
391: *
392: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
393: D( NL+1 ) = ZERO
394: DO 10 I = 1, N
395: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
396: ORGNRM = ABS( D( I ) )
397: END IF
398: 10 CONTINUE
399: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
400: ALPHA = ALPHA / ORGNRM
401: BETA = BETA / ORGNRM
402: *
403: * Sort and Deflate singular values.
404: *
405: CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
406: $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
407: $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
408: $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
409: $ INFO )
410: *
411: * Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
412: *
413: CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
414: $ WORK( ISIGMA ), WORK( IW ), INFO )
415: *
416: * Report the possible convergence failure.
417: *
418: IF( INFO.NE.0 ) THEN
419: RETURN
420: END IF
421: *
422: * Save the poles if ICOMPQ = 1.
423: *
424: IF( ICOMPQ.EQ.1 ) THEN
425: CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 )
426: CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
427: END IF
428: *
429: * Unscale.
430: *
431: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
432: *
433: * Prepare the IDXQ sorting permutation.
434: *
435: N1 = K
436: N2 = N - K
437: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
438: *
439: RETURN
440: *
441: * End of DLASD6
442: *
443: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack