Annotation of rpl/lapack/lapack/dlasd0.f, revision 1.7
1.1 bertrand 1: SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
2: $ WORK, INFO )
3: *
1.5 bertrand 4: * -- LAPACK auxiliary routine (version 3.2.2) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.5 bertrand 7: * June 2010
1.1 bertrand 8: *
9: * .. Scalar Arguments ..
10: INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE
11: * ..
12: * .. Array Arguments ..
13: INTEGER IWORK( * )
14: DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
15: $ WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * Using a divide and conquer approach, DLASD0 computes the singular
22: * value decomposition (SVD) of a real upper bidiagonal N-by-M
23: * matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
24: * The algorithm computes orthogonal matrices U and VT such that
25: * B = U * S * VT. The singular values S are overwritten on D.
26: *
27: * A related subroutine, DLASDA, computes only the singular values,
28: * and optionally, the singular vectors in compact form.
29: *
30: * Arguments
31: * =========
32: *
33: * N (input) INTEGER
34: * On entry, the row dimension of the upper bidiagonal matrix.
35: * This is also the dimension of the main diagonal array D.
36: *
37: * SQRE (input) INTEGER
38: * Specifies the column dimension of the bidiagonal matrix.
39: * = 0: The bidiagonal matrix has column dimension M = N;
40: * = 1: The bidiagonal matrix has column dimension M = N+1;
41: *
42: * D (input/output) DOUBLE PRECISION array, dimension (N)
43: * On entry D contains the main diagonal of the bidiagonal
44: * matrix.
45: * On exit D, if INFO = 0, contains its singular values.
46: *
47: * E (input) DOUBLE PRECISION array, dimension (M-1)
48: * Contains the subdiagonal entries of the bidiagonal matrix.
49: * On exit, E has been destroyed.
50: *
51: * U (output) DOUBLE PRECISION array, dimension at least (LDQ, N)
52: * On exit, U contains the left singular vectors.
53: *
54: * LDU (input) INTEGER
55: * On entry, leading dimension of U.
56: *
57: * VT (output) DOUBLE PRECISION array, dimension at least (LDVT, M)
58: * On exit, VT' contains the right singular vectors.
59: *
60: * LDVT (input) INTEGER
61: * On entry, leading dimension of VT.
62: *
63: * SMLSIZ (input) INTEGER
64: * On entry, maximum size of the subproblems at the
65: * bottom of the computation tree.
66: *
67: * IWORK (workspace) INTEGER work array.
68: * Dimension must be at least (8 * N)
69: *
70: * WORK (workspace) DOUBLE PRECISION work array.
71: * Dimension must be at least (3 * M**2 + 2 * M)
72: *
73: * INFO (output) INTEGER
74: * = 0: successful exit.
75: * < 0: if INFO = -i, the i-th argument had an illegal value.
1.5 bertrand 76: * > 0: if INFO = 1, a singular value did not converge
1.1 bertrand 77: *
78: * Further Details
79: * ===============
80: *
81: * Based on contributions by
82: * Ming Gu and Huan Ren, Computer Science Division, University of
83: * California at Berkeley, USA
84: *
85: * =====================================================================
86: *
87: * .. Local Scalars ..
88: INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
89: $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
90: $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
91: DOUBLE PRECISION ALPHA, BETA
92: * ..
93: * .. External Subroutines ..
94: EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA
95: * ..
96: * .. Executable Statements ..
97: *
98: * Test the input parameters.
99: *
100: INFO = 0
101: *
102: IF( N.LT.0 ) THEN
103: INFO = -1
104: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
105: INFO = -2
106: END IF
107: *
108: M = N + SQRE
109: *
110: IF( LDU.LT.N ) THEN
111: INFO = -6
112: ELSE IF( LDVT.LT.M ) THEN
113: INFO = -8
114: ELSE IF( SMLSIZ.LT.3 ) THEN
115: INFO = -9
116: END IF
117: IF( INFO.NE.0 ) THEN
118: CALL XERBLA( 'DLASD0', -INFO )
119: RETURN
120: END IF
121: *
122: * If the input matrix is too small, call DLASDQ to find the SVD.
123: *
124: IF( N.LE.SMLSIZ ) THEN
125: CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
126: $ LDU, WORK, INFO )
127: RETURN
128: END IF
129: *
130: * Set up the computation tree.
131: *
132: INODE = 1
133: NDIML = INODE + N
134: NDIMR = NDIML + N
135: IDXQ = NDIMR + N
136: IWK = IDXQ + N
137: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
138: $ IWORK( NDIMR ), SMLSIZ )
139: *
140: * For the nodes on bottom level of the tree, solve
141: * their subproblems by DLASDQ.
142: *
143: NDB1 = ( ND+1 ) / 2
144: NCC = 0
145: DO 30 I = NDB1, ND
146: *
147: * IC : center row of each node
148: * NL : number of rows of left subproblem
149: * NR : number of rows of right subproblem
150: * NLF: starting row of the left subproblem
151: * NRF: starting row of the right subproblem
152: *
153: I1 = I - 1
154: IC = IWORK( INODE+I1 )
155: NL = IWORK( NDIML+I1 )
156: NLP1 = NL + 1
157: NR = IWORK( NDIMR+I1 )
158: NRP1 = NR + 1
159: NLF = IC - NL
160: NRF = IC + 1
161: SQREI = 1
162: CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
163: $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
164: $ U( NLF, NLF ), LDU, WORK, INFO )
165: IF( INFO.NE.0 ) THEN
166: RETURN
167: END IF
168: ITEMP = IDXQ + NLF - 2
169: DO 10 J = 1, NL
170: IWORK( ITEMP+J ) = J
171: 10 CONTINUE
172: IF( I.EQ.ND ) THEN
173: SQREI = SQRE
174: ELSE
175: SQREI = 1
176: END IF
177: NRP1 = NR + SQREI
178: CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
179: $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
180: $ U( NRF, NRF ), LDU, WORK, INFO )
181: IF( INFO.NE.0 ) THEN
182: RETURN
183: END IF
184: ITEMP = IDXQ + IC
185: DO 20 J = 1, NR
186: IWORK( ITEMP+J-1 ) = J
187: 20 CONTINUE
188: 30 CONTINUE
189: *
190: * Now conquer each subproblem bottom-up.
191: *
192: DO 50 LVL = NLVL, 1, -1
193: *
194: * Find the first node LF and last node LL on the
195: * current level LVL.
196: *
197: IF( LVL.EQ.1 ) THEN
198: LF = 1
199: LL = 1
200: ELSE
201: LF = 2**( LVL-1 )
202: LL = 2*LF - 1
203: END IF
204: DO 40 I = LF, LL
205: IM1 = I - 1
206: IC = IWORK( INODE+IM1 )
207: NL = IWORK( NDIML+IM1 )
208: NR = IWORK( NDIMR+IM1 )
209: NLF = IC - NL
210: IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
211: SQREI = SQRE
212: ELSE
213: SQREI = 1
214: END IF
215: IDXQC = IDXQ + NLF - 1
216: ALPHA = D( IC )
217: BETA = E( IC )
218: CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
219: $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
220: $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
221: IF( INFO.NE.0 ) THEN
222: RETURN
223: END IF
224: 40 CONTINUE
225: 50 CONTINUE
226: *
227: RETURN
228: *
229: * End of DLASD0
230: *
231: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack