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