Annotation of rpl/lapack/lapack/dlasd1.f, revision 1.2
1.1 bertrand 1: SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
2: $ IDXQ, IWORK, WORK, INFO )
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, LDU, LDVT, NL, NR, SQRE
11: DOUBLE PRECISION ALPHA, BETA
12: * ..
13: * .. Array Arguments ..
14: INTEGER IDXQ( * ), IWORK( * )
15: DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
22: * where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
23: *
24: * A related subroutine DLASD7 handles the case in which the singular
25: * values (and the singular vectors in factored form) are desired.
26: *
27: * DLASD1 computes the SVD as follows:
28: *
29: * ( D1(in) 0 0 0 )
30: * B = U(in) * ( Z1' a Z2' b ) * VT(in)
31: * ( 0 0 D2(in) 0 )
32: *
33: * = U(out) * ( D(out) 0) * VT(out)
34: *
35: * where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
36: * with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
37: * elsewhere; and the entry b is empty if SQRE = 0.
38: *
39: * The left singular vectors of the original matrix are stored in U, and
40: * the transpose of the right singular vectors are stored in VT, and the
41: * singular values are in D. The algorithm consists of three stages:
42: *
43: * The first stage consists of deflating the size of the problem
44: * when there are multiple singular values or when there are zeros in
45: * the Z vector. For each such occurence the dimension of the
46: * secular equation problem is reduced by one. This stage is
47: * performed by the routine DLASD2.
48: *
49: * The second stage consists of calculating the updated
50: * singular values. This is done by finding the square roots of the
51: * roots of the secular equation via the routine DLASD4 (as called
52: * by DLASD3). This routine also calculates the singular vectors of
53: * the current problem.
54: *
55: * The final stage consists of computing the updated singular vectors
56: * directly using the updated singular values. The singular vectors
57: * for the current problem are multiplied with the singular vectors
58: * from the overall problem.
59: *
60: * Arguments
61: * =========
62: *
63: * NL (input) INTEGER
64: * The row dimension of the upper block. NL >= 1.
65: *
66: * NR (input) INTEGER
67: * The row dimension of the lower block. NR >= 1.
68: *
69: * SQRE (input) INTEGER
70: * = 0: the lower block is an NR-by-NR square matrix.
71: * = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
72: *
73: * The bidiagonal matrix has row dimension N = NL + NR + 1,
74: * and column dimension M = N + SQRE.
75: *
76: * D (input/output) DOUBLE PRECISION array,
77: * dimension (N = NL+NR+1).
78: * On entry D(1:NL,1:NL) contains the singular values of the
79: * upper block; and D(NL+2:N) contains the singular values of
80: * the lower block. On exit D(1:N) contains the singular values
81: * of the modified matrix.
82: *
83: * ALPHA (input/output) DOUBLE PRECISION
84: * Contains the diagonal element associated with the added row.
85: *
86: * BETA (input/output) DOUBLE PRECISION
87: * Contains the off-diagonal element associated with the added
88: * row.
89: *
90: * U (input/output) DOUBLE PRECISION array, dimension(LDU,N)
91: * On entry U(1:NL, 1:NL) contains the left singular vectors of
92: * the upper block; U(NL+2:N, NL+2:N) contains the left singular
93: * vectors of the lower block. On exit U contains the left
94: * singular vectors of the bidiagonal matrix.
95: *
96: * LDU (input) INTEGER
97: * The leading dimension of the array U. LDU >= max( 1, N ).
98: *
99: * VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M)
100: * where M = N + SQRE.
101: * On entry VT(1:NL+1, 1:NL+1)' contains the right singular
102: * vectors of the upper block; VT(NL+2:M, NL+2:M)' contains
103: * the right singular vectors of the lower block. On exit
104: * VT' contains the right singular vectors of the
105: * bidiagonal matrix.
106: *
107: * LDVT (input) INTEGER
108: * The leading dimension of the array VT. LDVT >= max( 1, M ).
109: *
110: * IDXQ (output) INTEGER array, dimension(N)
111: * This contains the permutation which will reintegrate the
112: * subproblem just solved back into sorted order, i.e.
113: * D( IDXQ( I = 1, N ) ) will be in ascending order.
114: *
115: * IWORK (workspace) INTEGER array, dimension( 4 * N )
116: *
117: * WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
118: *
119: * INFO (output) INTEGER
120: * = 0: successful exit.
121: * < 0: if INFO = -i, the i-th argument had an illegal value.
122: * > 0: if INFO = 1, an singular value did not converge
123: *
124: * Further Details
125: * ===============
126: *
127: * Based on contributions by
128: * Ming Gu and Huan Ren, Computer Science Division, University of
129: * California at Berkeley, USA
130: *
131: * =====================================================================
132: *
133: * .. Parameters ..
134: *
135: DOUBLE PRECISION ONE, ZERO
136: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
137: * ..
138: * .. Local Scalars ..
139: INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
140: $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
141: DOUBLE PRECISION ORGNRM
142: * ..
143: * .. External Subroutines ..
144: EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
145: * ..
146: * .. Intrinsic Functions ..
147: INTRINSIC ABS, MAX
148: * ..
149: * .. Executable Statements ..
150: *
151: * Test the input parameters.
152: *
153: INFO = 0
154: *
155: IF( NL.LT.1 ) THEN
156: INFO = -1
157: ELSE IF( NR.LT.1 ) THEN
158: INFO = -2
159: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
160: INFO = -3
161: END IF
162: IF( INFO.NE.0 ) THEN
163: CALL XERBLA( 'DLASD1', -INFO )
164: RETURN
165: END IF
166: *
167: N = NL + NR + 1
168: M = N + SQRE
169: *
170: * The following values are for bookkeeping purposes only. They are
171: * integer pointers which indicate the portion of the workspace
172: * used by a particular array in DLASD2 and DLASD3.
173: *
174: LDU2 = N
175: LDVT2 = M
176: *
177: IZ = 1
178: ISIGMA = IZ + M
179: IU2 = ISIGMA + N
180: IVT2 = IU2 + LDU2*N
181: IQ = IVT2 + LDVT2*M
182: *
183: IDX = 1
184: IDXC = IDX + N
185: COLTYP = IDXC + N
186: IDXP = COLTYP + N
187: *
188: * Scale.
189: *
190: ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
191: D( NL+1 ) = ZERO
192: DO 10 I = 1, N
193: IF( ABS( D( I ) ).GT.ORGNRM ) THEN
194: ORGNRM = ABS( D( I ) )
195: END IF
196: 10 CONTINUE
197: CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
198: ALPHA = ALPHA / ORGNRM
199: BETA = BETA / ORGNRM
200: *
201: * Deflate singular values.
202: *
203: CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
204: $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
205: $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
206: $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
207: *
208: * Solve Secular Equation and update singular vectors.
209: *
210: LDQ = K
211: CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
212: $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
213: $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
214: $ INFO )
215: IF( INFO.NE.0 ) THEN
216: RETURN
217: END IF
218: *
219: * Unscale.
220: *
221: CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
222: *
223: * Prepare the IDXQ sorting permutation.
224: *
225: N1 = K
226: N2 = N - K
227: CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
228: *
229: RETURN
230: *
231: * End of DLASD1
232: *
233: 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