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