Annotation of rpl/lapack/lapack/dlasd1.f, revision 1.18
1.13 bertrand 1: *> \brief \b DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc.
1.10 bertrand 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
1.17 bertrand 63: *> the Z vector. For each such occurrence the dimension of the
1.10 bertrand 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: *>
1.17 bertrand 159: *> \param[in,out] IDXQ
1.10 bertrand 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: *
1.17 bertrand 193: *> \date June 2016
1.10 bertrand 194: *
195: *> \ingroup auxOTHERauxiliary
196: *
197: *> \par Contributors:
198: * ==================
199: *>
200: *> Ming Gu and Huan Ren, Computer Science Division, University of
201: *> California at Berkeley, USA
202: *>
203: * =====================================================================
1.1 bertrand 204: SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
205: $ IDXQ, IWORK, WORK, INFO )
206: *
1.17 bertrand 207: * -- LAPACK auxiliary routine (version 3.6.1) --
1.1 bertrand 208: * -- LAPACK is a software package provided by Univ. of Tennessee, --
209: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.17 bertrand 210: * June 2016
1.1 bertrand 211: *
212: * .. Scalar Arguments ..
213: INTEGER INFO, LDU, LDVT, NL, NR, SQRE
214: DOUBLE PRECISION ALPHA, BETA
215: * ..
216: * .. Array Arguments ..
217: INTEGER IDXQ( * ), IWORK( * )
218: DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: *
225: DOUBLE PRECISION ONE, ZERO
226: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
227: * ..
228: * .. Local Scalars ..
229: INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
230: $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
231: DOUBLE PRECISION ORGNRM
232: * ..
233: * .. External Subroutines ..
234: EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
235: * ..
236: * .. Intrinsic Functions ..
237: INTRINSIC ABS, MAX
238: * ..
239: * .. Executable Statements ..
240: *
241: * Test the input parameters.
242: *
243: INFO = 0
244: *
245: IF( NL.LT.1 ) THEN
246: INFO = -1
247: ELSE IF( NR.LT.1 ) THEN
248: INFO = -2
249: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
250: INFO = -3
251: END IF
252: IF( INFO.NE.0 ) THEN
253: CALL XERBLA( 'DLASD1', -INFO )
254: RETURN
255: END IF
256: *
257: N = NL + NR + 1
258: M = N + SQRE
259: *
260: * The following values are for bookkeeping purposes only. They are
261: * integer pointers which indicate the portion of the workspace
262: * used by a particular array in DLASD2 and DLASD3.
263: *
264: LDU2 = N
265: LDVT2 = M
266: *
267: IZ = 1
268: ISIGMA = IZ + M
269: IU2 = ISIGMA + N
270: IVT2 = IU2 + LDU2*N
271: IQ = IVT2 + LDVT2*M
272: *
273: IDX = 1
274: IDXC = IDX + N
275: COLTYP = IDXC + N
276: IDXP = COLTYP + N
277: *
278: * Scale.
279: *
280: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
281: D( NL+1 ) = ZERO
282: DO 10 I = 1, N
283: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
284: ORGNRM = ABS( D( I ) )
285: END IF
286: 10 CONTINUE
287: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
288: ALPHA = ALPHA / ORGNRM
289: BETA = BETA / ORGNRM
290: *
291: * Deflate singular values.
292: *
293: CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
294: $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
295: $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
296: $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
297: *
298: * Solve Secular Equation and update singular vectors.
299: *
300: LDQ = K
301: CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
302: $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
303: $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
304: $ INFO )
1.16 bertrand 305: *
306: * Report the convergence failure.
307: *
1.1 bertrand 308: IF( INFO.NE.0 ) THEN
309: RETURN
310: END IF
311: *
312: * Unscale.
313: *
314: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
315: *
316: * Prepare the IDXQ sorting permutation.
317: *
318: N1 = K
319: N2 = N - K
320: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
321: *
322: RETURN
323: *
324: * End of DLASD1
325: *
326: 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