Annotation of rpl/lapack/lapack/dlasd0.f, revision 1.10
1.10 ! bertrand 1: *> \brief \b DLASD0
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLASD0 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd0.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd0.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd0.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
! 22: * WORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * INTEGER IWORK( * )
! 29: * DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
! 30: * $ WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> Using a divide and conquer approach, DLASD0 computes the singular
! 40: *> value decomposition (SVD) of a real upper bidiagonal N-by-M
! 41: *> matrix B with diagonal D and offdiagonal E, where M = N + SQRE.
! 42: *> The algorithm computes orthogonal matrices U and VT such that
! 43: *> B = U * S * VT. The singular values S are overwritten on D.
! 44: *>
! 45: *> A related subroutine, DLASDA, computes only the singular values,
! 46: *> and optionally, the singular vectors in compact form.
! 47: *> \endverbatim
! 48: *
! 49: * Arguments:
! 50: * ==========
! 51: *
! 52: *> \param[in] N
! 53: *> \verbatim
! 54: *> N is INTEGER
! 55: *> On entry, the row dimension of the upper bidiagonal matrix.
! 56: *> This is also the dimension of the main diagonal array D.
! 57: *> \endverbatim
! 58: *>
! 59: *> \param[in] SQRE
! 60: *> \verbatim
! 61: *> SQRE is INTEGER
! 62: *> Specifies the column dimension of the bidiagonal matrix.
! 63: *> = 0: The bidiagonal matrix has column dimension M = N;
! 64: *> = 1: The bidiagonal matrix has column dimension M = N+1;
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in,out] D
! 68: *> \verbatim
! 69: *> D is DOUBLE PRECISION array, dimension (N)
! 70: *> On entry D contains the main diagonal of the bidiagonal
! 71: *> matrix.
! 72: *> On exit D, if INFO = 0, contains its singular values.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] E
! 76: *> \verbatim
! 77: *> E is DOUBLE PRECISION array, dimension (M-1)
! 78: *> Contains the subdiagonal entries of the bidiagonal matrix.
! 79: *> On exit, E has been destroyed.
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[out] U
! 83: *> \verbatim
! 84: *> U is DOUBLE PRECISION array, dimension at least (LDQ, N)
! 85: *> On exit, U contains the left singular vectors.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in] LDU
! 89: *> \verbatim
! 90: *> LDU is INTEGER
! 91: *> On entry, leading dimension of U.
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[out] VT
! 95: *> \verbatim
! 96: *> VT is DOUBLE PRECISION array, dimension at least (LDVT, M)
! 97: *> On exit, VT**T contains the right singular vectors.
! 98: *> \endverbatim
! 99: *>
! 100: *> \param[in] LDVT
! 101: *> \verbatim
! 102: *> LDVT is INTEGER
! 103: *> On entry, leading dimension of VT.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[in] SMLSIZ
! 107: *> \verbatim
! 108: *> SMLSIZ is INTEGER
! 109: *> On entry, maximum size of the subproblems at the
! 110: *> bottom of the computation tree.
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[out] IWORK
! 114: *> \verbatim
! 115: *> IWORK is INTEGER work array.
! 116: *> Dimension must be at least (8 * N)
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[out] WORK
! 120: *> \verbatim
! 121: *> WORK is DOUBLE PRECISION work array.
! 122: *> Dimension must be at least (3 * M**2 + 2 * M)
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[out] INFO
! 126: *> \verbatim
! 127: *> INFO is INTEGER
! 128: *> = 0: successful exit.
! 129: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 130: *> > 0: if INFO = 1, a singular value did not converge
! 131: *> \endverbatim
! 132: *
! 133: * Authors:
! 134: * ========
! 135: *
! 136: *> \author Univ. of Tennessee
! 137: *> \author Univ. of California Berkeley
! 138: *> \author Univ. of Colorado Denver
! 139: *> \author NAG Ltd.
! 140: *
! 141: *> \date November 2011
! 142: *
! 143: *> \ingroup auxOTHERauxiliary
! 144: *
! 145: *> \par Contributors:
! 146: * ==================
! 147: *>
! 148: *> Ming Gu and Huan Ren, Computer Science Division, University of
! 149: *> California at Berkeley, USA
! 150: *>
! 151: * =====================================================================
1.1 bertrand 152: SUBROUTINE DLASD0( N, SQRE, D, E, U, LDU, VT, LDVT, SMLSIZ, IWORK,
153: $ WORK, INFO )
154: *
1.10 ! bertrand 155: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 156: * -- LAPACK is a software package provided by Univ. of Tennessee, --
157: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.10 ! bertrand 158: * November 2011
1.1 bertrand 159: *
160: * .. Scalar Arguments ..
161: INTEGER INFO, LDU, LDVT, N, SMLSIZ, SQRE
162: * ..
163: * .. Array Arguments ..
164: INTEGER IWORK( * )
165: DOUBLE PRECISION D( * ), E( * ), U( LDU, * ), VT( LDVT, * ),
166: $ WORK( * )
167: * ..
168: *
169: * =====================================================================
170: *
171: * .. Local Scalars ..
172: INTEGER I, I1, IC, IDXQ, IDXQC, IM1, INODE, ITEMP, IWK,
173: $ J, LF, LL, LVL, M, NCC, ND, NDB1, NDIML, NDIMR,
174: $ NL, NLF, NLP1, NLVL, NR, NRF, NRP1, SQREI
175: DOUBLE PRECISION ALPHA, BETA
176: * ..
177: * .. External Subroutines ..
178: EXTERNAL DLASD1, DLASDQ, DLASDT, XERBLA
179: * ..
180: * .. Executable Statements ..
181: *
182: * Test the input parameters.
183: *
184: INFO = 0
185: *
186: IF( N.LT.0 ) THEN
187: INFO = -1
188: ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
189: INFO = -2
190: END IF
191: *
192: M = N + SQRE
193: *
194: IF( LDU.LT.N ) THEN
195: INFO = -6
196: ELSE IF( LDVT.LT.M ) THEN
197: INFO = -8
198: ELSE IF( SMLSIZ.LT.3 ) THEN
199: INFO = -9
200: END IF
201: IF( INFO.NE.0 ) THEN
202: CALL XERBLA( 'DLASD0', -INFO )
203: RETURN
204: END IF
205: *
206: * If the input matrix is too small, call DLASDQ to find the SVD.
207: *
208: IF( N.LE.SMLSIZ ) THEN
209: CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDVT, U, LDU, U,
210: $ LDU, WORK, INFO )
211: RETURN
212: END IF
213: *
214: * Set up the computation tree.
215: *
216: INODE = 1
217: NDIML = INODE + N
218: NDIMR = NDIML + N
219: IDXQ = NDIMR + N
220: IWK = IDXQ + N
221: CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
222: $ IWORK( NDIMR ), SMLSIZ )
223: *
224: * For the nodes on bottom level of the tree, solve
225: * their subproblems by DLASDQ.
226: *
227: NDB1 = ( ND+1 ) / 2
228: NCC = 0
229: DO 30 I = NDB1, ND
230: *
231: * IC : center row of each node
232: * NL : number of rows of left subproblem
233: * NR : number of rows of right subproblem
234: * NLF: starting row of the left subproblem
235: * NRF: starting row of the right subproblem
236: *
237: I1 = I - 1
238: IC = IWORK( INODE+I1 )
239: NL = IWORK( NDIML+I1 )
240: NLP1 = NL + 1
241: NR = IWORK( NDIMR+I1 )
242: NRP1 = NR + 1
243: NLF = IC - NL
244: NRF = IC + 1
245: SQREI = 1
246: CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), E( NLF ),
247: $ VT( NLF, NLF ), LDVT, U( NLF, NLF ), LDU,
248: $ U( NLF, NLF ), LDU, WORK, INFO )
249: IF( INFO.NE.0 ) THEN
250: RETURN
251: END IF
252: ITEMP = IDXQ + NLF - 2
253: DO 10 J = 1, NL
254: IWORK( ITEMP+J ) = J
255: 10 CONTINUE
256: IF( I.EQ.ND ) THEN
257: SQREI = SQRE
258: ELSE
259: SQREI = 1
260: END IF
261: NRP1 = NR + SQREI
262: CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), E( NRF ),
263: $ VT( NRF, NRF ), LDVT, U( NRF, NRF ), LDU,
264: $ U( NRF, NRF ), LDU, WORK, INFO )
265: IF( INFO.NE.0 ) THEN
266: RETURN
267: END IF
268: ITEMP = IDXQ + IC
269: DO 20 J = 1, NR
270: IWORK( ITEMP+J-1 ) = J
271: 20 CONTINUE
272: 30 CONTINUE
273: *
274: * Now conquer each subproblem bottom-up.
275: *
276: DO 50 LVL = NLVL, 1, -1
277: *
278: * Find the first node LF and last node LL on the
279: * current level LVL.
280: *
281: IF( LVL.EQ.1 ) THEN
282: LF = 1
283: LL = 1
284: ELSE
285: LF = 2**( LVL-1 )
286: LL = 2*LF - 1
287: END IF
288: DO 40 I = LF, LL
289: IM1 = I - 1
290: IC = IWORK( INODE+IM1 )
291: NL = IWORK( NDIML+IM1 )
292: NR = IWORK( NDIMR+IM1 )
293: NLF = IC - NL
294: IF( ( SQRE.EQ.0 ) .AND. ( I.EQ.LL ) ) THEN
295: SQREI = SQRE
296: ELSE
297: SQREI = 1
298: END IF
299: IDXQC = IDXQ + NLF - 1
300: ALPHA = D( IC )
301: BETA = E( IC )
302: CALL DLASD1( NL, NR, SQREI, D( NLF ), ALPHA, BETA,
303: $ U( NLF, NLF ), LDU, VT( NLF, NLF ), LDVT,
304: $ IWORK( IDXQC ), IWORK( IWK ), WORK, INFO )
305: IF( INFO.NE.0 ) THEN
306: RETURN
307: END IF
308: 40 CONTINUE
309: 50 CONTINUE
310: *
311: RETURN
312: *
313: * End of DLASD0
314: *
315: END
CVSweb interface <joel.bertrand@systella.fr>