vendor/gonum.org/v1/gonum/lapack/internal/testdata/dsterftest/dsterf.f

   1 *> \brief \b DSTERF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSTERF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSTERF( N, D, E, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INFO, N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       DOUBLE PRECISION   D( * ), E( * )
  28 *       ..
  29 *
  30 *
  31 *> \par Purpose:
  32 *  =============
  33 *>
  34 *> \verbatim
  35 *>
  36 *> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
  37 *> using the Pal-Walker-Kahan variant of the QL or QR algorithm.
  38 *> \endverbatim
  39 *
  40 *  Arguments:
  41 *  ==========
  42 *
  43 *> \param[in] N
  44 *> \verbatim
  45 *>          N is INTEGER
  46 *>          The order of the matrix.  N >= 0.
  47 *> \endverbatim
  48 *>
  49 *> \param[in,out] D
  50 *> \verbatim
  51 *>          D is DOUBLE PRECISION array, dimension (N)
  52 *>          On entry, the n diagonal elements of the tridiagonal matrix.
  53 *>          On exit, if INFO = 0, the eigenvalues in ascending order.
  54 *> \endverbatim
  55 *>
  56 *> \param[in,out] E
  57 *> \verbatim
  58 *>          E is DOUBLE PRECISION array, dimension (N-1)
  59 *>          On entry, the (n-1) subdiagonal elements of the tridiagonal
  60 *>          matrix.
  61 *>          On exit, E has been destroyed.
  62 *> \endverbatim
  63 *>
  64 *> \param[out] INFO
  65 *> \verbatim
  66 *>          INFO is INTEGER
  67 *>          = 0:  successful exit
  68 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  69 *>          > 0:  the algorithm failed to find all of the eigenvalues in
  70 *>                a total of 30*N iterations; if INFO = i, then i
  71 *>                elements of E have not converged to zero.
  72 *> \endverbatim
  73 *
  74 *  Authors:
  75 *  ========
  76 *
  77 *> \author Univ. of Tennessee
  78 *> \author Univ. of California Berkeley
  79 *> \author Univ. of Colorado Denver
  80 *> \author NAG Ltd.
  81 *
  82 *> \date November 2011
  83 *
  84 *> \ingroup auxOTHERcomputational
  85 *
  86 *  =====================================================================
  87       SUBROUTINE DSTERF( N, D, E, INFO )
  88 *
  89 *  -- LAPACK computational routine (version 3.4.0) --
  90 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  91 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  92 *     November 2011
  93 *
  94 *     .. Scalar Arguments ..
  95       INTEGER            INFO, N
  96 *     ..
  97 *     .. Array Arguments ..
  98       DOUBLE PRECISION   D( * ), E( * )
  99 *     ..
 100 *
 101 *  =====================================================================
 102 *
 103 *     .. Parameters ..
 104       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
 105       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
 106      $                   THREE = 3.0D0 )
 107       INTEGER            MAXIT
 108       PARAMETER          ( MAXIT = 30 )
 109 *     ..
 110 *     .. Local Scalars ..
 111       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,
 112      $                   NMAXIT
 113       DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
 114      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
 115      $                   SIGMA, SSFMAX, SSFMIN, RMAX
 116 *     ..
 117 *     .. External Functions ..
 118       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
 119       EXTERNAL           DLAMCH, DLANST, DLAPY2
 120 *     ..
 121 *     .. External Subroutines ..
 122       EXTERNAL           DLAE2, DLASCL, DLASRT, XERBLA
 123 *     ..
 124 *     .. Intrinsic Functions ..
 125       INTRINSIC          ABS, SIGN, SQRT
 126 *     ..
 127 *     .. Executable Statements ..
 128 *
 129 *     Test the input parameters.
 130 *
 131       INFO = 0
 132 *
 133 *     Quick return if possible
 134 *
 135       IF( N.LT.0 ) THEN
 136          INFO = -1
 137          CALL XERBLA( 'DSTERF', -INFO )
 138          RETURN
 139       END IF
 140       IF( N.LE.1 )
 141      $   RETURN
 142 *
 143 *     Determine the unit roundoff for this environment.
 144 *
 145       EPS = DLAMCH( 'E' )
 146       EPS2 = EPS**2
 147       SAFMIN = DLAMCH( 'S' )
 148       SAFMAX = ONE / SAFMIN
 149       SSFMAX = SQRT( SAFMAX ) / THREE
 150       SSFMIN = SQRT( SAFMIN ) / EPS2
 151       RMAX = DLAMCH( 'O' )
 152 *
 153 *     Compute the eigenvalues of the tridiagonal matrix.
 154 *
 155       NMAXIT = N*MAXIT
 156       SIGMA = ZERO
 157       JTOT = 0
 158 *
 159 *     Determine where the matrix splits and choose QL or QR iteration
 160 *     for each block, according to whether top or bottom diagonal
 161 *     element is smaller.
 162 *
 163       L1 = 1
 164 *
 165    10 CONTINUE
 166       print *, "l1 = ", l1
 167       IF( L1.GT.N ) THEN
 168         print *, "going to 170"
 169         GO TO 170
 170       end if
 171       IF( L1.GT.1 )
 172      $   E( L1-1 ) = ZERO
 173       DO 20 M = L1, N - 1
 174          IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
 175      $       1 ) ) ) )*EPS ) THEN
 176             E( M ) = ZERO
 177             GO TO 30
 178          END IF
 179    20 CONTINUE
 180       M = N
 181 *
 182    30 CONTINUE
 183       print *, "30, d"
 184       print *, d(1:n)
 185       L = L1
 186       LSV = L
 187       LEND = M
 188       LENDSV = LEND
 189       L1 = M + 1
 190       IF( LEND.EQ.L )
 191      $   GO TO 10
 192 *
 193 *     Scale submatrix in rows and columns L to LEND
 194 *
 195       ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) )
 196       ISCALE = 0
 197       IF( ANORM.EQ.ZERO )
 198      $   GO TO 10
 199       IF( (ANORM.GT.SSFMAX) ) THEN
 200          ISCALE = 1
 201          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
 202      $                INFO )
 203          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
 204      $                INFO )
 205       ELSE IF( ANORM.LT.SSFMIN ) THEN
 206          ISCALE = 2
 207          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
 208      $                INFO )
 209          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
 210      $                INFO )
 211       END IF
 212 *
 213       DO 40 I = L, LEND - 1
 214          E( I ) = E( I )**2
 215    40 CONTINUE
 216 *
 217 *     Choose between QL and QR iteration
 218 *
 219       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
 220          LEND = LSV
 221          L = LENDSV
 222       END IF
 223 *
 224       IF( LEND.GE.L ) THEN
 225         print *, "ql, d"
 226         print *, d(1:n)
 227 *
 228 *        QL Iteration
 229 *
 230 *        Look for small subdiagonal element.
 231 *
 232    50    CONTINUE
 233          IF( L.NE.LEND ) THEN
 234             DO 60 M = L, LEND - 1
 235                IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )
 236      $            GO TO 70
 237    60       CONTINUE
 238          END IF
 239          M = LEND
 240 *
 241    70    CONTINUE
 242          IF( M.LT.LEND )
 243      $      E( M ) = ZERO
 244          P = D( L )
 245          IF( M.EQ.L )
 246      $      GO TO 90
 247 *
 248 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
 249 *        eigenvalues.
 250 *
 251          IF( M.EQ.L+1 ) THEN
 252             RTE = SQRT( E( L ) )
 253             CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
 254             D( L ) = RT1
 255             D( L+1 ) = RT2
 256             E( L ) = ZERO
 257             L = L + 2
 258             IF( L.LE.LEND )
 259      $         GO TO 50
 260             GO TO 150
 261          END IF
 262 *
 263          IF( JTOT.EQ.NMAXIT )
 264      $      GO TO 150
 265          JTOT = JTOT + 1
 266 *
 267 *        Form shift.
 268 *
 269          RTE = SQRT( E( L ) )
 270          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
 271          R = DLAPY2( SIGMA, ONE )
 272          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 273 *
 274          C = ONE
 275          S = ZERO
 276          GAMMA = D( M ) - SIGMA
 277          P = GAMMA*GAMMA
 278 *
 279 *        Inner loop
 280 *
 281           print *, "inner loop d before"
 282           print *, d(1:n)
 283
 284          DO 80 I = M - 1, L, -1
 285             print *, "inner loop"
 286             print *, "ei", e(i)
 287             BB = E( I )
 288             R = P + BB
 289             print *, "bb,p,r"
 290             print *, bb,p,r
 291             IF( I.NE.M-1 ) THEN
 292               print *, s,r
 293               E( I+1 ) = S*R
 294             end if
 295             OLDC = C
 296             C = P / R
 297             S = BB / R
 298             OLDGAM = GAMMA
 299             print *, "di", d(i)
 300             ALPHA = D( I )
 301             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 302             print *,"og, a, ga", OLDGAM, ALPHA, GAMMA
 303             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
 304             IF( C.NE.ZERO ) THEN
 305                P = ( GAMMA*GAMMA ) / C
 306             ELSE
 307                P = OLDC*BB
 308             END IF
 309             print *, "p, gamma = ", p,GAMMA
 310    80    CONTINUE
 311 *
 312          E( L ) = S*P
 313          D( L ) = SIGMA + GAMMA
 314
 315         print *, "inner loop d after"
 316         print *, d(1:n)
 317          GO TO 50
 318 *
 319 *        Eigenvalue found.
 320 *
 321    90    CONTINUE
 322          D( L ) = P
 323 *
 324          L = L + 1
 325          IF( L.LE.LEND )
 326      $      GO TO 50
 327          GO TO 150
 328 *
 329       ELSE
 330 *
 331 *        QR Iteration
 332 *
 333 *        Look for small superdiagonal element.
 334 *
 335   100    CONTINUE
 336          DO 110 M = L, LEND + 1, -1
 337             IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )
 338      $         GO TO 120
 339   110    CONTINUE
 340          M = LEND
 341 *
 342   120    CONTINUE
 343          IF( M.GT.LEND )
 344      $      E( M-1 ) = ZERO
 345          P = D( L )
 346          IF( M.EQ.L )
 347      $      GO TO 140
 348 *
 349 *        If remaining matrix is 2 by 2, use DLAE2 to compute its
 350 *        eigenvalues.
 351 *
 352          IF( M.EQ.L-1 ) THEN
 353             RTE = SQRT( E( L-1 ) )
 354             CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
 355             D( L ) = RT1
 356             D( L-1 ) = RT2
 357             E( L-1 ) = ZERO
 358             L = L - 2
 359             IF( L.GE.LEND )
 360      $         GO TO 100
 361             GO TO 150
 362          END IF
 363 *
 364          IF( JTOT.EQ.NMAXIT )
 365      $      GO TO 150
 366          JTOT = JTOT + 1
 367 *
 368 *        Form shift.
 369 *
 370          RTE = SQRT( E( L-1 ) )
 371          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
 372          R = DLAPY2( SIGMA, ONE )
 373          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
 374 *
 375          C = ONE
 376          S = ZERO
 377          GAMMA = D( M ) - SIGMA
 378          P = GAMMA*GAMMA
 379 *
 380 *        Inner loop
 381 *
 382          DO 130 I = M, L - 1
 383             BB = E( I )
 384             R = P + BB
 385             IF( I.NE.M )
 386      $         E( I-1 ) = S*R
 387             OLDC = C
 388             C = P / R
 389             S = BB / R
 390             OLDGAM = GAMMA
 391             ALPHA = D( I+1 )
 392             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
 393             D( I ) = OLDGAM + ( ALPHA-GAMMA )
 394             IF( C.NE.ZERO ) THEN
 395                P = ( GAMMA*GAMMA ) / C
 396             ELSE
 397                P = OLDC*BB
 398             END IF
 399   130    CONTINUE
 400 *
 401          E( L-1 ) = S*P
 402          D( L ) = SIGMA + GAMMA
 403          GO TO 100
 404 *
 405 *        Eigenvalue found.
 406 *
 407   140    CONTINUE
 408          D( L ) = P
 409 *
 410          L = L - 1
 411          IF( L.GE.LEND )
 412      $      GO TO 100
 413          GO TO 150
 414 *
 415       END IF
 416 *
 417 *     Undo scaling if necessary
 418 *
 419   150 CONTINUE
 420       IF( ISCALE.EQ.1 )
 421      $   CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
 422      $                D( LSV ), N, INFO )
 423       IF( ISCALE.EQ.2 )
 424      $   CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
 425      $                D( LSV ), N, INFO )
 426 *
 427 *     Check for no convergence to an eigenvalue after a total
 428 *     of N*MAXIT iterations.
 429 *
 430       IF( JTOT.LT.NMAXIT )
 431      $   GO TO 10
 432       DO 160 I = 1, N - 1
 433          IF( E( I ).NE.ZERO )
 434      $      INFO = INFO + 1
 435   160 CONTINUE
 436       GO TO 180
 437 *
 438 *     Sort eigenvalues in increasing order.
 439 *
 440   170 CONTINUE
 441       CALL DLASRT( 'I', N, D, INFO )
 442 *
 443   180 CONTINUE
 444       RETURN
 445 *
 446 *     End of DSTERF
 447 *
 448       END