vendor/gonum.org/v1/gonum/lapack/internal/testdata/dlasqtest/dlasq1.f

   1 *> \brief \b DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASQ1 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq1.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq1.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq1.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DLASQ1 computes the singular values of a real N-by-N bidiagonal
  37 *> matrix with diagonal D and off-diagonal E. The singular values
  38 *> are computed to high relative accuracy, in the absence of
  39 *> denormalization, underflow and overflow. The algorithm was first
  40 *> presented in
  41 *>
  42 *> "Accurate singular values and differential qd algorithms" by K. V.
  43 *> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
  44 *> 1994,
  45 *>
  46 *> and the present implementation is described in "An implementation of
  47 *> the dqds Algorithm (Positive Case)", LAPACK Working Note.
  48 *> \endverbatim
  49 *
  50 *  Arguments:
  51 *  ==========
  52 *
  53 *> \param[in] N
  54 *> \verbatim
  55 *>          N is INTEGER
  56 *>        The number of rows and columns in the matrix. N >= 0.
  57 *> \endverbatim
  58 *>
  59 *> \param[in,out] D
  60 *> \verbatim
  61 *>          D is DOUBLE PRECISION array, dimension (N)
  62 *>        On entry, D contains the diagonal elements of the
  63 *>        bidiagonal matrix whose SVD is desired. On normal exit,
  64 *>        D contains the singular values in decreasing order.
  65 *> \endverbatim
  66 *>
  67 *> \param[in,out] E
  68 *> \verbatim
  69 *>          E is DOUBLE PRECISION array, dimension (N)
  70 *>        On entry, elements E(1:N-1) contain the off-diagonal elements
  71 *>        of the bidiagonal matrix whose SVD is desired.
  72 *>        On exit, E is overwritten.
  73 *> \endverbatim
  74 *>
  75 *> \param[out] WORK
  76 *> \verbatim
  77 *>          WORK is DOUBLE PRECISION array, dimension (4*N)
  78 *> \endverbatim
  79 *>
  80 *> \param[out] INFO
  81 *> \verbatim
  82 *>          INFO is INTEGER
  83 *>        = 0: successful exit
  84 *>        < 0: if INFO = -i, the i-th argument had an illegal value
  85 *>        > 0: the algorithm failed
  86 *>             = 1, a split was marked by a positive value in E
  87 *>             = 2, current block of Z not diagonalized after 100*N
  88 *>                  iterations (in inner while loop)  On exit D and E
  89 *>                  represent a matrix with the same singular values
  90 *>                  which the calling subroutine could use to finish the
  91 *>                  computation, or even feed back into DLASQ1
  92 *>             = 3, termination criterion of outer while loop not met
  93 *>                  (program created more than N unreduced blocks)
  94 *> \endverbatim
  95 *
  96 *  Authors:
  97 *  ========
  98 *
  99 *> \author Univ. of Tennessee
 100 *> \author Univ. of California Berkeley
 101 *> \author Univ. of Colorado Denver
 102 *> \author NAG Ltd.
 103 *
 104 *> \date September 2012
 105 *
 106 *> \ingroup auxOTHERcomputational
 107 *
 108 *  =====================================================================
 109       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
 110 *
 111 *  -- LAPACK computational routine (version 3.4.2) --
 112 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 113 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 114 *     September 2012
 115 *
 116 *     .. Scalar Arguments ..
 117       INTEGER            INFO, N
 118 *     ..
 119 *     .. Array Arguments ..
 120       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
 121 *     ..
 122 *
 123 *  =====================================================================
 124 *
 125 *     .. Parameters ..
 126       DOUBLE PRECISION   ZERO
 127       PARAMETER          ( ZERO = 0.0D0 )
 128 *     ..
 129 *     .. Local Scalars ..
 130       INTEGER            I, IINFO
 131       DOUBLE PRECISION   EPS, SCALE, SAFMIN, SIGMN, SIGMX
 132 *     ..
 133 *     .. External Subroutines ..
 134       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
 135 *     ..
 136 *     .. External Functions ..
 137       DOUBLE PRECISION   DLAMCH
 138       EXTERNAL           DLAMCH
 139 *     ..
 140 *     .. Intrinsic Functions ..
 141       INTRINSIC          ABS, MAX, SQRT
 142 *     ..
 143 *     .. Executable Statements ..
 144 *
 145       INFO = 0
 146       IF( N.LT.0 ) THEN
 147          INFO = -2
 148          CALL XERBLA( 'DLASQ1', -INFO )
 149          RETURN
 150       ELSE IF( N.EQ.0 ) THEN
 151          RETURN
 152       ELSE IF( N.EQ.1 ) THEN
 153          D( 1 ) = ABS( D( 1 ) )
 154          RETURN
 155       ELSE IF( N.EQ.2 ) THEN
 156          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
 157          D( 1 ) = SIGMX
 158          D( 2 ) = SIGMN
 159          RETURN
 160       END IF
 161 *
 162 *     Estimate the largest singular value.
 163 *
 164       SIGMX = ZERO
 165       DO 10 I = 1, N - 1
 166          D( I ) = ABS( D( I ) )
 167          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
 168    10 CONTINUE
 169       D( N ) = ABS( D( N ) )
 170 *
 171 *     Early return if SIGMX is zero (matrix is already diagonal).
 172 *
 173       IF( SIGMX.EQ.ZERO ) THEN
 174          CALL DLASRT( 'D', N, D, IINFO )
 175          RETURN
 176       END IF
 177 *
 178       DO 20 I = 1, N
 179          SIGMX = MAX( SIGMX, D( I ) )
 180    20 CONTINUE
 181 *
 182 *     Copy D and E into WORK (in the Z format) and scale (squaring the
 183 *     input data makes scaling by a power of the radix pointless).
 184 *
 185       EPS = DLAMCH( 'Precision' )
 186       SAFMIN = DLAMCH( 'Safe minimum' )
 187       SCALE = SQRT( EPS / SAFMIN )
 188
 189       CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
 190       CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
 191       CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
 192      $             IINFO )
 193 *
 194 *     Compute the q's and e's.
 195 *
 196       DO 30 I = 1, 2*N - 1
 197          WORK( I ) = WORK( I )**2
 198    30 CONTINUE
 199       WORK( 2*N ) = ZERO
 200 *
 201
 202       CALL DLASQ2( N, WORK, INFO )
 203 *
 204       IF( INFO.EQ.0 ) THEN
 205          DO 40 I = 1, N
 206             D( I ) = SQRT( WORK( I ) )
 207    40    CONTINUE
 208          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
 209       ELSE IF( INFO.EQ.2 ) THEN
 210 *
 211 *     Maximum number of iterations exceeded.  Move data from WORK
 212 *     into D and E so the calling subroutine can try to finish
 213 *
 214          DO I = 1, N
 215             D( I ) = SQRT( WORK( 2*I-1 ) )
 216             E( I ) = SQRT( WORK( 2*I ) )
 217          END DO
 218          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
 219          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
 220       END IF
 221 *
 222       RETURN
 223 *
 224 *     End of DLASQ1
 225 *
 226       END