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
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