LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slasq1.f
Go to the documentation of this file.
00001 *> \brief \b SLASQ1
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLASQ1 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasq1.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasq1.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasq1.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLASQ1( N, D, E, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       REAL               D( * ), E( * ), WORK( * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> SLASQ1 computes the singular values of a real N-by-N bidiagonal
00037 *> matrix with diagonal D and off-diagonal E. The singular values
00038 *> are computed to high relative accuracy, in the absence of
00039 *> denormalization, underflow and overflow. The algorithm was first
00040 *> presented in
00041 *>
00042 *> "Accurate singular values and differential qd algorithms" by K. V.
00043 *> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
00044 *> 1994,
00045 *>
00046 *> and the present implementation is described in "An implementation of
00047 *> the dqds Algorithm (Positive Case)", LAPACK Working Note.
00048 *> \endverbatim
00049 *
00050 *  Arguments:
00051 *  ==========
00052 *
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>        The number of rows and columns in the matrix. N >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in,out] D
00060 *> \verbatim
00061 *>          D is REAL array, dimension (N)
00062 *>        On entry, D contains the diagonal elements of the
00063 *>        bidiagonal matrix whose SVD is desired. On normal exit,
00064 *>        D contains the singular values in decreasing order.
00065 *> \endverbatim
00066 *>
00067 *> \param[in,out] E
00068 *> \verbatim
00069 *>          E is REAL array, dimension (N)
00070 *>        On entry, elements E(1:N-1) contain the off-diagonal elements
00071 *>        of the bidiagonal matrix whose SVD is desired.
00072 *>        On exit, E is overwritten.
00073 *> \endverbatim
00074 *>
00075 *> \param[out] WORK
00076 *> \verbatim
00077 *>          WORK is REAL array, dimension (4*N)
00078 *> \endverbatim
00079 *>
00080 *> \param[out] INFO
00081 *> \verbatim
00082 *>          INFO is INTEGER
00083 *>        = 0: successful exit
00084 *>        < 0: if INFO = -i, the i-th argument had an illegal value
00085 *>        > 0: the algorithm failed
00086 *>             = 1, a split was marked by a positive value in E
00087 *>             = 2, current block of Z not diagonalized after 100*N
00088 *>                  iterations (in inner while loop)  On exit D and E
00089 *>                  represent a matrix with the same singular values
00090 *>                  which the calling subroutine could use to finish the
00091 *>                  computation, or even feed back into SLASQ1
00092 *>             = 3, termination criterion of outer while loop not met 
00093 *>                  (program created more than N unreduced blocks)
00094 *> \endverbatim
00095 *
00096 *  Authors:
00097 *  ========
00098 *
00099 *> \author Univ. of Tennessee 
00100 *> \author Univ. of California Berkeley 
00101 *> \author Univ. of Colorado Denver 
00102 *> \author NAG Ltd. 
00103 *
00104 *> \date November 2011
00105 *
00106 *> \ingroup auxOTHERcomputational
00107 *
00108 *  =====================================================================
00109       SUBROUTINE SLASQ1( N, D, E, WORK, INFO )
00110 *
00111 *  -- LAPACK computational routine (version 3.4.0) --
00112 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00113 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00114 *     November 2011
00115 *
00116 *     .. Scalar Arguments ..
00117       INTEGER            INFO, N
00118 *     ..
00119 *     .. Array Arguments ..
00120       REAL               D( * ), E( * ), WORK( * )
00121 *     ..
00122 *
00123 *  =====================================================================
00124 *
00125 *     .. Parameters ..
00126       REAL               ZERO
00127       PARAMETER          ( ZERO = 0.0E0 )
00128 *     ..
00129 *     .. Local Scalars ..
00130       INTEGER            I, IINFO
00131       REAL               EPS, SCALE, SAFMIN, SIGMN, SIGMX
00132 *     ..
00133 *     .. External Subroutines ..
00134       EXTERNAL           SCOPY, SLAS2, SLASCL, SLASQ2, SLASRT, XERBLA
00135 *     ..
00136 *     .. External Functions ..
00137       REAL               SLAMCH
00138       EXTERNAL           SLAMCH
00139 *     ..
00140 *     .. Intrinsic Functions ..
00141       INTRINSIC          ABS, MAX, SQRT
00142 *     ..
00143 *     .. Executable Statements ..
00144 *
00145       INFO = 0
00146       IF( N.LT.0 ) THEN
00147          INFO = -2
00148          CALL XERBLA( 'SLASQ1', -INFO )
00149          RETURN
00150       ELSE IF( N.EQ.0 ) THEN
00151          RETURN
00152       ELSE IF( N.EQ.1 ) THEN
00153          D( 1 ) = ABS( D( 1 ) )
00154          RETURN
00155       ELSE IF( N.EQ.2 ) THEN
00156          CALL SLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
00157          D( 1 ) = SIGMX
00158          D( 2 ) = SIGMN
00159          RETURN
00160       END IF
00161 *
00162 *     Estimate the largest singular value.
00163 *
00164       SIGMX = ZERO
00165       DO 10 I = 1, N - 1
00166          D( I ) = ABS( D( I ) )
00167          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
00168    10 CONTINUE
00169       D( N ) = ABS( D( N ) )
00170 *
00171 *     Early return if SIGMX is zero (matrix is already diagonal).
00172 *
00173       IF( SIGMX.EQ.ZERO ) THEN
00174          CALL SLASRT( 'D', N, D, IINFO )
00175          RETURN
00176       END IF
00177 *
00178       DO 20 I = 1, N
00179          SIGMX = MAX( SIGMX, D( I ) )
00180    20 CONTINUE
00181 *
00182 *     Copy D and E into WORK (in the Z format) and scale (squaring the
00183 *     input data makes scaling by a power of the radix pointless).
00184 *
00185       EPS = SLAMCH( 'Precision' )
00186       SAFMIN = SLAMCH( 'Safe minimum' )
00187       SCALE = SQRT( EPS / SAFMIN )
00188       CALL SCOPY( N, D, 1, WORK( 1 ), 2 )
00189       CALL SCOPY( N-1, E, 1, WORK( 2 ), 2 )
00190       CALL SLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
00191      $             IINFO )
00192 *         
00193 *     Compute the q's and e's.
00194 *
00195       DO 30 I = 1, 2*N - 1
00196          WORK( I ) = WORK( I )**2
00197    30 CONTINUE
00198       WORK( 2*N ) = ZERO
00199 *
00200       CALL SLASQ2( N, WORK, INFO )
00201 *
00202       IF( INFO.EQ.0 ) THEN
00203          DO 40 I = 1, N
00204             D( I ) = SQRT( WORK( I ) )
00205    40    CONTINUE
00206          CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
00207       ELSE IF( INFO.EQ.2 ) THEN
00208 *
00209 *     Maximum number of iterations exceeded.  Move data from WORK
00210 *     into D and E so the calling subroutine can try to finish
00211 *
00212          DO I = 1, N
00213             D( I ) = SQRT( WORK( 2*I-1 ) )
00214             E( I ) = SQRT( WORK( 2*I ) )
00215          END DO
00216          CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
00217          CALL SLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
00218       END IF
00219 *
00220       RETURN
00221 *
00222 *     End of SLASQ1
00223 *
00224       END
 All Files Functions