Numérotation des lignes :

dsterf
C DSTERF    SOURCE    BP208322  22/09/16    21:15:06     11454          *> \brief \b DSTERF**  =========== DOCUMENTATION ===========** Online html documentation available at*            http://www.netlib.org/lapack/explore-html/**> \htmlonly*> Download DSTERF + dependencies*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsterf.f">*> [TGZ]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsterf.f">*> [ZIP]&lt;/a>*> &lt;a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsterf.f">*> [TXT]&lt;/a>*> \endhtmlonly**  Definition:*  ===========**       SUBROUTINE DSTERF( N, D, E, INFO )**       .. Scalar Arguments ..*       INTEGER            INFO, N*       ..*       .. Array Arguments ..*       REAL*8   D( * ), E( * )*       ..***> \par Purpose:*  =============*>*> \verbatim*>*> DSTERF computes all eigenvalues of a symmetric tridiagonal matrix*> using the Pal-Walker-Kahan variant of the QL or QR algorithm.*> \endverbatim**  Arguments:*  ==========**> \param[in] N*> \verbatim*>          N is INTEGER*>          The order of the matrix.  N >= 0.*> \endverbatim*>*> \param[in,out] D*> \verbatim*>          D is REAL*8 array, dimension (N)*>          On entry, the n diagonal elements of the tridiagonal matrix.*>          On exit, if INFO = 0, the eigenvalues in ascending order.*> \endverbatim*>*> \param[in,out] E*> \verbatim*>          E is REAL*8 array, dimension (N-1)*>          On entry, the (n-1) subdiagonal elements of the tridiagonal*>          matrix.*>          On exit, E has been destroyed.*> \endverbatim*>*> \param[out] INFO*> \verbatim*>          INFO is INTEGER*>          = 0:  successful exit*>          &lt; 0:  if INFO = -i, the i-th argument had an illegal value*>          > 0:  the algorithm failed to find all of the eigenvalues in*>                a total of 30*N iterations; if INFO = i, then i*>                elements of E have not converged to zero.*> \endverbatim**  Authors:*  ========**> \author Univ. of Tennessee*> \author Univ. of California Berkeley*> \author Univ. of Colorado Denver*> \author NAG Ltd.**> \ingroup auxOTHERcomputational**  =====================================================================      SUBROUTINE DSTERF( N, D, E, INFO )**  -- LAPACK computational routine --*  -- LAPACK is a software package provided by Univ. of Tennessee,    --*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--**     .. Scalar Arguments ..      INTEGER            INFO, N*     ..*     .. Array Arguments ..      REAL*8   D( * ), E( * )*     ..**  =====================================================================**     .. Parameters ..      REAL*8   ZERO, ONE, TWO, THREE      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,     $THREE = 3.0D0 ) INTEGER MAXIT PARAMETER ( MAXIT = 30 )* ..* .. Local Scalars .. INTEGER I, ISCALE, JTOT, L, L1, LEND, LENDSV, LSV, M,$                   NMAXIT      REAL*8   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,     $OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,$                   SIGMA, SSFMAX, SSFMIN, RMAX*     ..*     .. External Functions ..      REAL*8   DLAMCH, DLANST, DLAPY2      EXTERNAL           DLAMCH, DLANST, DLAPY2*     ..*     .. External Subroutines ..      EXTERNAL           DLAE2, DLASCL, DLASRT, XERBLA*     ..*     .. Intrinsic Functions ..*      INTRINSIC          ABS, SIGN, SQRT*     ..*     .. Executable Statements ..**     Test the input parameters.*      INFO = 0**     Quick return if possible*      IF( N.LT.0 ) THEN         INFO = -1         CALL XERBLA( 'DSTERF', -INFO )         RETURN      END IF      IF( N.LE.1 )     $RETURN** Determine the unit roundoff for this environment.* EPS = DLAMCH( 'E' ) EPS2 = EPS**2 SAFMIN = DLAMCH( 'S' ) SAFMAX = ONE / SAFMIN SSFMAX = SQRT( SAFMAX ) / THREE SSFMIN = SQRT( SAFMIN ) / EPS2 RMAX = DLAMCH( 'O' )** Compute the eigenvalues of the tridiagonal matrix.* NMAXIT = N*MAXIT SIGMA = ZERO JTOT = 0** Determine where the matrix splits and choose QL or QR iteration* for each block, according to whether top or bottom diagonal* element is smaller.* L1 = 1* 10 CONTINUE IF( L1.GT.N )$   GO TO 170      IF( L1.GT.1 )     $E( L1-1 ) = ZERO DO 20 M = L1, N - 1 IF( ABS( E( M ) ).LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+$       1 ) ) ) )*EPS ) THEN            E( M ) = ZERO            GO TO 30         END IF   20 CONTINUE      M = N*   30 CONTINUE      L = L1      LSV = L      LEND = M      LENDSV = LEND      L1 = M + 1      IF( LEND.EQ.L )     $GO TO 10** Scale submatrix in rows and columns L to LEND* ANORM = DLANST( 'M', LEND-L+1, D( L ), E( L ) ) ISCALE = 0 IF( ANORM.EQ.ZERO )$   GO TO 10      IF( (ANORM.GT.SSFMAX) ) THEN         ISCALE = 1         CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,     $INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,$                INFO )      ELSE IF( ANORM.LT.SSFMIN ) THEN         ISCALE = 2         CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,     $INFO ) CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,$                INFO )      END IF*      DO 40 I = L, LEND - 1         E( I ) = E( I )**2   40 CONTINUE**     Choose between QL and QR iteration*      IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN         LEND = LSV         L = LENDSV      END IF*      IF( LEND.GE.L ) THEN**        QL Iteration**        Look for small subdiagonal element.*   50    CONTINUE         IF( L.NE.LEND ) THEN            DO 60 M = L, LEND - 1               IF( ABS( E( M ) ).LE.EPS2*ABS( D( M )*D( M+1 ) ) )     $GO TO 70 60 CONTINUE END IF M = LEND* 70 CONTINUE IF( M.LT.LEND )$      E( M ) = ZERO         P = D( L )         IF( M.EQ.L )     $GO TO 90** If remaining matrix is 2 by 2, use DLAE2 to compute its* eigenvalues.* IF( M.EQ.L+1 ) THEN RTE = SQRT( E( L ) ) CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 ) D( L ) = RT1 D( L+1 ) = RT2 E( L ) = ZERO L = L + 2 IF( L.LE.LEND )$         GO TO 50            GO TO 150         END IF*         IF( JTOT.EQ.NMAXIT )     $GO TO 150 JTOT = JTOT + 1** Form shift.* RTE = SQRT( E( L ) ) SIGMA = ( D( L+1 )-P ) / ( TWO*RTE ) R = DLAPY2( SIGMA, ONE ) SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )* C = ONE S = ZERO GAMMA = D( M ) - SIGMA P = GAMMA*GAMMA** Inner loop* DO 80 I = M - 1, L, -1 BB = E( I ) R = P + BB IF( I.NE.M-1 )$         E( I+1 ) = S*R            OLDC = C            C = P / R            S = BB / R            OLDGAM = GAMMA            ALPHA = D( I )            GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM            D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )            IF( C.NE.ZERO ) THEN               P = ( GAMMA*GAMMA ) / C            ELSE               P = OLDC*BB            END IF   80    CONTINUE*         E( L ) = S*P         D( L ) = SIGMA + GAMMA         GO TO 50**        Eigenvalue found.*   90    CONTINUE         D( L ) = P*         L = L + 1         IF( L.LE.LEND )     $GO TO 50 GO TO 150* ELSE** QR Iteration** Look for small superdiagonal element.* 100 CONTINUE DO 110 M = L, LEND + 1, -1 IF( ABS( E( M-1 ) ).LE.EPS2*ABS( D( M )*D( M-1 ) ) )$         GO TO 120  110    CONTINUE         M = LEND*  120    CONTINUE         IF( M.GT.LEND )     $E( M-1 ) = ZERO P = D( L ) IF( M.EQ.L )$      GO TO 140**        If remaining matrix is 2 by 2, use DLAE2 to compute its*        eigenvalues.*         IF( M.EQ.L-1 ) THEN            RTE = SQRT( E( L-1 ) )            CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )            D( L ) = RT1            D( L-1 ) = RT2            E( L-1 ) = ZERO            L = L - 2            IF( L.GE.LEND )     $GO TO 100 GO TO 150 END IF* IF( JTOT.EQ.NMAXIT )$      GO TO 150         JTOT = JTOT + 1**        Form shift.*         RTE = SQRT( E( L-1 ) )         SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )         R = DLAPY2( SIGMA, ONE )         SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )*         C = ONE         S = ZERO         GAMMA = D( M ) - SIGMA         P = GAMMA*GAMMA**        Inner loop*         DO 130 I = M, L - 1            BB = E( I )            R = P + BB            IF( I.NE.M )     $E( I-1 ) = S*R OLDC = C C = P / R S = BB / R OLDGAM = GAMMA ALPHA = D( I+1 ) GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM D( I ) = OLDGAM + ( ALPHA-GAMMA ) IF( C.NE.ZERO ) THEN P = ( GAMMA*GAMMA ) / C ELSE P = OLDC*BB END IF 130 CONTINUE* E( L-1 ) = S*P D( L ) = SIGMA + GAMMA GO TO 100** Eigenvalue found.* 140 CONTINUE D( L ) = P* L = L - 1 IF( L.GE.LEND )$      GO TO 100         GO TO 150*      END IF**     Undo scaling if necessary*  150 CONTINUE      IF( ISCALE.EQ.1 )     $CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,$                D( LSV ), N, INFO )      IF( ISCALE.EQ.2 )     $CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,$                D( LSV ), N, INFO )**     Check for no convergence to an eigenvalue after a total*     of N*MAXIT iterations.*      IF( JTOT.LT.NMAXIT )     $GO TO 10 DO 160 I = 1, N - 1 IF( E( I ).NE.ZERO )$      INFO = INFO + 1  160 CONTINUE      GO TO 180**     Sort eigenvalues in increasing order.*  170 CONTINUE      CALL DLASRT( 'I', N, D, INFO )*  180 CONTINUE      RETURN**     End of DSTERF*      END  

© Cast3M 2003 - Tous droits réservés.
Mentions légales