1: *> \brief \b DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. 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 DLASD1 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd1.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd1.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd1.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
22: * IDXQ, IWORK, WORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INFO, LDU, LDVT, NL, NR, SQRE
26: * DOUBLE PRECISION ALPHA, BETA
27: * ..
28: * .. Array Arguments ..
29: * INTEGER IDXQ( * ), IWORK( * )
30: * DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
40: *> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
41: *>
42: *> A related subroutine DLASD7 handles the case in which the singular
43: *> values (and the singular vectors in factored form) are desired.
44: *>
45: *> DLASD1 computes the SVD as follows:
46: *>
47: *> ( D1(in) 0 0 0 )
48: *> B = U(in) * ( Z1**T a Z2**T b ) * VT(in)
49: *> ( 0 0 D2(in) 0 )
50: *>
51: *> = U(out) * ( D(out) 0) * VT(out)
52: *>
53: *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
54: *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
55: *> elsewhere; and the entry b is empty if SQRE = 0.
56: *>
57: *> The left singular vectors of the original matrix are stored in U, and
58: *> the transpose of the right singular vectors are stored in VT, and the
59: *> singular values are in D. The algorithm consists of three stages:
60: *>
61: *> The first stage consists of deflating the size of the problem
62: *> when there are multiple singular values or when there are zeros in
63: *> the Z vector. For each such occurrence the dimension of the
64: *> secular equation problem is reduced by one. This stage is
65: *> performed by the routine DLASD2.
66: *>
67: *> The second stage consists of calculating the updated
68: *> singular values. This is done by finding the square roots of the
69: *> roots of the secular equation via the routine DLASD4 (as called
70: *> by DLASD3). This routine also calculates the singular vectors of
71: *> the current problem.
72: *>
73: *> The final stage consists of computing the updated singular vectors
74: *> directly using the updated singular values. The singular vectors
75: *> for the current problem are multiplied with the singular vectors
76: *> from the overall problem.
77: *> \endverbatim
78: *
79: * Arguments:
80: * ==========
81: *
82: *> \param[in] NL
83: *> \verbatim
84: *> NL is INTEGER
85: *> The row dimension of the upper block. NL >= 1.
86: *> \endverbatim
87: *>
88: *> \param[in] NR
89: *> \verbatim
90: *> NR is INTEGER
91: *> The row dimension of the lower block. NR >= 1.
92: *> \endverbatim
93: *>
94: *> \param[in] SQRE
95: *> \verbatim
96: *> SQRE is INTEGER
97: *> = 0: the lower block is an NR-by-NR square matrix.
98: *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
99: *>
100: *> The bidiagonal matrix has row dimension N = NL + NR + 1,
101: *> and column dimension M = N + SQRE.
102: *> \endverbatim
103: *>
104: *> \param[in,out] D
105: *> \verbatim
106: *> D is DOUBLE PRECISION array,
107: *> dimension (N = NL+NR+1).
108: *> On entry D(1:NL,1:NL) contains the singular values of the
109: *> upper block; and D(NL+2:N) contains the singular values of
110: *> the lower block. On exit D(1:N) contains the singular values
111: *> of the modified matrix.
112: *> \endverbatim
113: *>
114: *> \param[in,out] ALPHA
115: *> \verbatim
116: *> ALPHA is DOUBLE PRECISION
117: *> Contains the diagonal element associated with the added row.
118: *> \endverbatim
119: *>
120: *> \param[in,out] BETA
121: *> \verbatim
122: *> BETA is DOUBLE PRECISION
123: *> Contains the off-diagonal element associated with the added
124: *> row.
125: *> \endverbatim
126: *>
127: *> \param[in,out] U
128: *> \verbatim
129: *> U is DOUBLE PRECISION array, dimension(LDU,N)
130: *> On entry U(1:NL, 1:NL) contains the left singular vectors of
131: *> the upper block; U(NL+2:N, NL+2:N) contains the left singular
132: *> vectors of the lower block. On exit U contains the left
133: *> singular vectors of the bidiagonal matrix.
134: *> \endverbatim
135: *>
136: *> \param[in] LDU
137: *> \verbatim
138: *> LDU is INTEGER
139: *> The leading dimension of the array U. LDU >= max( 1, N ).
140: *> \endverbatim
141: *>
142: *> \param[in,out] VT
143: *> \verbatim
144: *> VT is DOUBLE PRECISION array, dimension(LDVT,M)
145: *> where M = N + SQRE.
146: *> On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
147: *> vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
148: *> the right singular vectors of the lower block. On exit
149: *> VT**T contains the right singular vectors of the
150: *> bidiagonal matrix.
151: *> \endverbatim
152: *>
153: *> \param[in] LDVT
154: *> \verbatim
155: *> LDVT is INTEGER
156: *> The leading dimension of the array VT. LDVT >= max( 1, M ).
157: *> \endverbatim
158: *>
159: *> \param[in,out] IDXQ
160: *> \verbatim
161: *> IDXQ is INTEGER array, dimension(N)
162: *> This contains the permutation which will reintegrate the
163: *> subproblem just solved back into sorted order, i.e.
164: *> D( IDXQ( I = 1, N ) ) will be in ascending order.
165: *> \endverbatim
166: *>
167: *> \param[out] IWORK
168: *> \verbatim
169: *> IWORK is INTEGER array, dimension( 4 * N )
170: *> \endverbatim
171: *>
172: *> \param[out] WORK
173: *> \verbatim
174: *> WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
175: *> \endverbatim
176: *>
177: *> \param[out] INFO
178: *> \verbatim
179: *> INFO is INTEGER
180: *> = 0: successful exit.
181: *> < 0: if INFO = -i, the i-th argument had an illegal value.
182: *> > 0: if INFO = 1, a singular value did not converge
183: *> \endverbatim
184: *
185: * Authors:
186: * ========
187: *
188: *> \author Univ. of Tennessee
189: *> \author Univ. of California Berkeley
190: *> \author Univ. of Colorado Denver
191: *> \author NAG Ltd.
192: *
193: *> \ingroup OTHERauxiliary
194: *
195: *> \par Contributors:
196: * ==================
197: *>
198: *> Ming Gu and Huan Ren, Computer Science Division, University of
199: *> California at Berkeley, USA
200: *>
201: * =====================================================================
202: SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
203: $ IDXQ, IWORK, WORK, INFO )
204: *
205: * -- LAPACK auxiliary routine --
206: * -- LAPACK is a software package provided by Univ. of Tennessee, --
207: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208: *
209: * .. Scalar Arguments ..
210: INTEGER INFO, LDU, LDVT, NL, NR, SQRE
211: DOUBLE PRECISION ALPHA, BETA
212: * ..
213: * .. Array Arguments ..
214: INTEGER IDXQ( * ), IWORK( * )
215: DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
216: * ..
217: *
218: * =====================================================================
219: *
220: * .. Parameters ..
221: *
222: DOUBLE PRECISION ONE, ZERO
223: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
224: * ..
225: * .. Local Scalars ..
226: INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
227: $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
228: DOUBLE PRECISION ORGNRM
229: * ..
230: * .. External Subroutines ..
231: EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
232: * ..
233: * .. Intrinsic Functions ..
234: INTRINSIC ABS, MAX
235: * ..
236: * .. Executable Statements ..
237: *
238: * Test the input parameters.
239: *
240: INFO = 0
241: *
242: IF( NL.LT.1 ) THEN
243: INFO = -1
244: ELSE IF( NR.LT.1 ) THEN
245: INFO = -2
246: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
247: INFO = -3
248: END IF
249: IF( INFO.NE.0 ) THEN
250: CALL XERBLA( 'DLASD1', -INFO )
251: RETURN
252: END IF
253: *
254: N = NL + NR + 1
255: M = N + SQRE
256: *
257: * The following values are for bookkeeping purposes only. They are
258: * integer pointers which indicate the portion of the workspace
259: * used by a particular array in DLASD2 and DLASD3.
260: *
261: LDU2 = N
262: LDVT2 = M
263: *
264: IZ = 1
265: ISIGMA = IZ + M
266: IU2 = ISIGMA + N
267: IVT2 = IU2 + LDU2*N
268: IQ = IVT2 + LDVT2*M
269: *
270: IDX = 1
271: IDXC = IDX + N
272: COLTYP = IDXC + N
273: IDXP = COLTYP + N
274: *
275: * Scale.
276: *
277: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
278: D( NL+1 ) = ZERO
279: DO 10 I = 1, N
280: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
281: ORGNRM = ABS( D( I ) )
282: END IF
283: 10 CONTINUE
284: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
285: ALPHA = ALPHA / ORGNRM
286: BETA = BETA / ORGNRM
287: *
288: * Deflate singular values.
289: *
290: CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
291: $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
292: $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
293: $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
294: *
295: * Solve Secular Equation and update singular vectors.
296: *
297: LDQ = K
298: CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
299: $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
300: $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
301: $ INFO )
302: *
303: * Report the convergence failure.
304: *
305: IF( INFO.NE.0 ) THEN
306: RETURN
307: END IF
308: *
309: * Unscale.
310: *
311: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
312: *
313: * Prepare the IDXQ sorting permutation.
314: *
315: N1 = K
316: N2 = N - K
317: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
318: *
319: RETURN
320: *
321: * End of DLASD1
322: *
323: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack