LAPACK 3.3.0

dlasq1.f

Go to the documentation of this file.
00001       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2)                                    --
00004 *
00005 *  -- Contributed by Osni Marques of the Lawrence Berkeley National   --
00006 *  -- Laboratory and Beresford Parlett of the Univ. of California at  --
00007 *  -- Berkeley                                                        --
00008 *  -- November 2008                                                   --
00009 *
00010 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00011 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00012 *
00013 *     .. Scalar Arguments ..
00014       INTEGER            INFO, N
00015 *     ..
00016 *     .. Array Arguments ..
00017       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  DLASQ1 computes the singular values of a real N-by-N bidiagonal
00024 *  matrix with diagonal D and off-diagonal E. The singular values
00025 *  are computed to high relative accuracy, in the absence of
00026 *  denormalization, underflow and overflow. The algorithm was first
00027 *  presented in
00028 *
00029 *  "Accurate singular values and differential qd algorithms" by K. V.
00030 *  Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
00031 *  1994,
00032 *
00033 *  and the present implementation is described in "An implementation of
00034 *  the dqds Algorithm (Positive Case)", LAPACK Working Note.
00035 *
00036 *  Arguments
00037 *  =========
00038 *
00039 *  N     (input) INTEGER
00040 *        The number of rows and columns in the matrix. N >= 0.
00041 *
00042 *  D     (input/output) DOUBLE PRECISION array, dimension (N)
00043 *        On entry, D contains the diagonal elements of the
00044 *        bidiagonal matrix whose SVD is desired. On normal exit,
00045 *        D contains the singular values in decreasing order.
00046 *
00047 *  E     (input/output) DOUBLE PRECISION array, dimension (N)
00048 *        On entry, elements E(1:N-1) contain the off-diagonal elements
00049 *        of the bidiagonal matrix whose SVD is desired.
00050 *        On exit, E is overwritten.
00051 *
00052 *  WORK  (workspace) DOUBLE PRECISION array, dimension (4*N)
00053 *
00054 *  INFO  (output) INTEGER
00055 *        = 0: successful exit
00056 *        < 0: if INFO = -i, the i-th argument had an illegal value
00057 *        > 0: the algorithm failed
00058 *             = 1, a split was marked by a positive value in E
00059 *             = 2, current block of Z not diagonalized after 30*N
00060 *                  iterations (in inner while loop)
00061 *             = 3, termination criterion of outer while loop not met 
00062 *                  (program created more than N unreduced blocks)
00063 *
00064 *  =====================================================================
00065 *
00066 *     .. Parameters ..
00067       DOUBLE PRECISION   ZERO
00068       PARAMETER          ( ZERO = 0.0D0 )
00069 *     ..
00070 *     .. Local Scalars ..
00071       INTEGER            I, IINFO
00072       DOUBLE PRECISION   EPS, SCALE, SAFMIN, SIGMN, SIGMX
00073 *     ..
00074 *     .. External Subroutines ..
00075       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
00076 *     ..
00077 *     .. External Functions ..
00078       DOUBLE PRECISION   DLAMCH
00079       EXTERNAL           DLAMCH
00080 *     ..
00081 *     .. Intrinsic Functions ..
00082       INTRINSIC          ABS, MAX, SQRT
00083 *     ..
00084 *     .. Executable Statements ..
00085 *
00086       INFO = 0
00087       IF( N.LT.0 ) THEN
00088          INFO = -2
00089          CALL XERBLA( 'DLASQ1', -INFO )
00090          RETURN
00091       ELSE IF( N.EQ.0 ) THEN
00092          RETURN
00093       ELSE IF( N.EQ.1 ) THEN
00094          D( 1 ) = ABS( D( 1 ) )
00095          RETURN
00096       ELSE IF( N.EQ.2 ) THEN
00097          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
00098          D( 1 ) = SIGMX
00099          D( 2 ) = SIGMN
00100          RETURN
00101       END IF
00102 *
00103 *     Estimate the largest singular value.
00104 *
00105       SIGMX = ZERO
00106       DO 10 I = 1, N - 1
00107          D( I ) = ABS( D( I ) )
00108          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
00109    10 CONTINUE
00110       D( N ) = ABS( D( N ) )
00111 *
00112 *     Early return if SIGMX is zero (matrix is already diagonal).
00113 *
00114       IF( SIGMX.EQ.ZERO ) THEN
00115          CALL DLASRT( 'D', N, D, IINFO )
00116          RETURN
00117       END IF
00118 *
00119       DO 20 I = 1, N
00120          SIGMX = MAX( SIGMX, D( I ) )
00121    20 CONTINUE
00122 *
00123 *     Copy D and E into WORK (in the Z format) and scale (squaring the
00124 *     input data makes scaling by a power of the radix pointless).
00125 *
00126       EPS = DLAMCH( 'Precision' )
00127       SAFMIN = DLAMCH( 'Safe minimum' )
00128       SCALE = SQRT( EPS / SAFMIN )
00129       CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
00130       CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
00131       CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
00132      $             IINFO )
00133 *         
00134 *     Compute the q's and e's.
00135 *
00136       DO 30 I = 1, 2*N - 1
00137          WORK( I ) = WORK( I )**2
00138    30 CONTINUE
00139       WORK( 2*N ) = ZERO
00140 *
00141       CALL DLASQ2( N, WORK, INFO )
00142 *
00143       IF( INFO.EQ.0 ) THEN
00144          DO 40 I = 1, N
00145             D( I ) = SQRT( WORK( I ) )
00146    40    CONTINUE
00147          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
00148       END IF
00149 *
00150       RETURN
00151 *
00152 *     End of DLASQ1
00153 *
00154       END
 All Files Functions