Annotation of rpl/lapack/lapack/dlasq1.f, revision 1.4
1.1 bertrand 1: SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: *
5: * -- Contributed by Osni Marques of the Lawrence Berkeley National --
6: * -- Laboratory and Beresford Parlett of the Univ. of California at --
7: * -- Berkeley --
8: * -- November 2008 --
9: *
10: * -- LAPACK is a software package provided by Univ. of Tennessee, --
11: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
12: *
13: * .. Scalar Arguments ..
14: INTEGER INFO, N
15: * ..
16: * .. Array Arguments ..
17: DOUBLE PRECISION D( * ), E( * ), WORK( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * DLASQ1 computes the singular values of a real N-by-N bidiagonal
24: * matrix with diagonal D and off-diagonal E. The singular values
25: * are computed to high relative accuracy, in the absence of
26: * denormalization, underflow and overflow. The algorithm was first
27: * presented in
28: *
29: * "Accurate singular values and differential qd algorithms" by K. V.
30: * Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
31: * 1994,
32: *
33: * and the present implementation is described in "An implementation of
34: * the dqds Algorithm (Positive Case)", LAPACK Working Note.
35: *
36: * Arguments
37: * =========
38: *
39: * N (input) INTEGER
40: * The number of rows and columns in the matrix. N >= 0.
41: *
42: * D (input/output) DOUBLE PRECISION array, dimension (N)
43: * On entry, D contains the diagonal elements of the
44: * bidiagonal matrix whose SVD is desired. On normal exit,
45: * D contains the singular values in decreasing order.
46: *
47: * E (input/output) DOUBLE PRECISION array, dimension (N)
48: * On entry, elements E(1:N-1) contain the off-diagonal elements
49: * of the bidiagonal matrix whose SVD is desired.
50: * On exit, E is overwritten.
51: *
52: * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
53: *
54: * INFO (output) INTEGER
55: * = 0: successful exit
56: * < 0: if INFO = -i, the i-th argument had an illegal value
57: * > 0: the algorithm failed
58: * = 1, a split was marked by a positive value in E
59: * = 2, current block of Z not diagonalized after 30*N
60: * iterations (in inner while loop)
61: * = 3, termination criterion of outer while loop not met
62: * (program created more than N unreduced blocks)
63: *
64: * =====================================================================
65: *
66: * .. Parameters ..
67: DOUBLE PRECISION ZERO
68: PARAMETER ( ZERO = 0.0D0 )
69: * ..
70: * .. Local Scalars ..
71: INTEGER I, IINFO
72: DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX
73: * ..
74: * .. External Subroutines ..
75: EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
76: * ..
77: * .. External Functions ..
78: DOUBLE PRECISION DLAMCH
79: EXTERNAL DLAMCH
80: * ..
81: * .. Intrinsic Functions ..
82: INTRINSIC ABS, MAX, SQRT
83: * ..
84: * .. Executable Statements ..
85: *
86: INFO = 0
87: IF( N.LT.0 ) THEN
88: INFO = -2
89: CALL XERBLA( 'DLASQ1', -INFO )
90: RETURN
91: ELSE IF( N.EQ.0 ) THEN
92: RETURN
93: ELSE IF( N.EQ.1 ) THEN
94: D( 1 ) = ABS( D( 1 ) )
95: RETURN
96: ELSE IF( N.EQ.2 ) THEN
97: CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
98: D( 1 ) = SIGMX
99: D( 2 ) = SIGMN
100: RETURN
101: END IF
102: *
103: * Estimate the largest singular value.
104: *
105: SIGMX = ZERO
106: DO 10 I = 1, N - 1
107: D( I ) = ABS( D( I ) )
108: SIGMX = MAX( SIGMX, ABS( E( I ) ) )
109: 10 CONTINUE
110: D( N ) = ABS( D( N ) )
111: *
112: * Early return if SIGMX is zero (matrix is already diagonal).
113: *
114: IF( SIGMX.EQ.ZERO ) THEN
115: CALL DLASRT( 'D', N, D, IINFO )
116: RETURN
117: END IF
118: *
119: DO 20 I = 1, N
120: SIGMX = MAX( SIGMX, D( I ) )
121: 20 CONTINUE
122: *
123: * Copy D and E into WORK (in the Z format) and scale (squaring the
124: * input data makes scaling by a power of the radix pointless).
125: *
126: EPS = DLAMCH( 'Precision' )
127: SAFMIN = DLAMCH( 'Safe minimum' )
128: SCALE = SQRT( EPS / SAFMIN )
129: CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
130: CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
131: CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
132: $ IINFO )
133: *
134: * Compute the q's and e's.
135: *
136: DO 30 I = 1, 2*N - 1
137: WORK( I ) = WORK( I )**2
138: 30 CONTINUE
139: WORK( 2*N ) = ZERO
140: *
141: CALL DLASQ2( N, WORK, INFO )
142: *
143: IF( INFO.EQ.0 ) THEN
144: DO 40 I = 1, N
145: D( I ) = SQRT( WORK( I ) )
146: 40 CONTINUE
147: CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
148: END IF
149: *
150: RETURN
151: *
152: * End of DLASQ1
153: *
154: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasq1.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack