Numérotation des lignes :

dstqrb
C DSTQRB    SOURCE    FANDEUR   22/05/02    21:15:16     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dstqrbcc\Description:c  Computes all eigenvalues and the last component of the eigenvectorsc  of a symmetric tridiagonal matrix using the implicit QL or QR method.cc  This is mostly a modification of the LAPACK routine dsteqr.c  See Remarks.cc\Usage:c  call dstqrbc     ( N, D, E, Z, WORK, INFO )cc\Argumentsc  N       Integer.  (INPUT)c          The number of rows and columns in the matrix.  N >= 0.cc  D       REAL*8 array, dimension (N).  (INPUT/OUTPUT)c          On entry, D contains the diagonal elements of thec          tridiagonal matrix.c          On exit, D contains the eigenvalues, in ascending order.c          If an error exit is made, the eigenvalues are correctc          for indices 1,2,...,INFO-1, but they are unordered andc          may not be the smallest eigenvalues of the matrix.cc  E       REAL*8 array, dimension (N-1).  (INPUT/OUTPUT)c          On entry, E contains the subdiagonal elements of thec          tridiagonal matrix in positions 1 through N-1.c          On exit, E has been destroyed.cc  Z       REAL*8 array, dimension (N).  (OUTPUT)c          On exit, Z contains the last row of the orthonormalc          eigenvector matrix of the symmetric tridiagonal matrix.c          If an error exit is made, Z contains the last row of thec          eigenvector matrix associated with the stored eigenvalues.cc  WORK    Double precision array, dimension (max(1,2*N-2)).  (WORKSPACE)c          Workspace used in accumulating the transformation forc          computing the last components of the eigenvectors.cc  INFO    Integer.  (OUTPUT)c          = 0:  normal return.c          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.c          > 0:  if INFO = +i, the i-th eigenvalue has not convergedc                              after a total of  30*N  iterations.cc\Remarksc  1. None.cc-----------------------------------------------------------------------cc\BeginLibcc\Local variables:c     xxxxxx  realcc\Routines called:c     daxpy   Level 1 BLAS that computes a vector triad.c     dcopy   Level 1 BLAS that copies one vector to another.c     dswap   Level 1 BLAS that swaps the contents of two vectors.c     lsame   LAPACK character comparison routine.c     dlae2   LAPACK routine that computes the eigenvalues of a 2-by-2c             symmetric matrix.c     dlaev2  LAPACK routine that eigendecomposition of a 2-by-2 symmetricc             matrix.c     dlamch  LAPACK routine that determines machine constants.c     dlanst  LAPACK routine that computes the norm of a matrix.c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.c     dlartg  LAPACK Givens rotation construction routine.c     dlascl  LAPACK routine for careful scaling of a matrix.c     dlaset  LAPACK matrix initialization routine.c     dlasr   LAPACK routine that applies an orthogonal transformation toc             a matrix.c     dlasrt  LAPACK sorting routine.c     dsteqr  LAPACK routine that computes eigenvalues and eigenvectorsc             of a symmetric tridiagonal matrix.c     xerbla  LAPACK error handler routine.cc\Authorsc     Danny Sorensen               Phuong Vuc     Richard Lehoucq              CRPC / Rice Universityc     Dept. of Computational &     Houston, Texasc     Applied Mathematicsc     Rice Universityc     Houston, Texascc\SCCS Information: @(#)c FILE: stqrb.F   SID: 2.5   DATE OF SID: 8/27/96   RELEASE: 2cc\Remarksc     1. Starting with version 2.5, this routine is a modified versionc        of LAPACK version 2.0 subroutine SSTEQR. No lines are deleted,c        only commeted out and new lines inserted.c        All lines commented out have "c$$" at the beginning.c Note that the LAPACK version 1.0 subroutine SSTEQR containedc bugs.cc\EndLibcc-----------------------------------------------------------------------c subroutine dstqrb ( n, d, e, z, work, info )cc %------------------%c | Scalar Arguments |c %------------------%c integer info, ncc %-----------------%c | Array Arguments |c %-----------------%c REAL*8 & d( n ), e( n-1 ), z( n ), work( 2*n-2 )cc .. parameters .. REAL*8 & zero, one, two, three parameter ( zero = 0.0D+0, one = 1.0D+0, & two = 2.0D+0, three = 3.0D+0 ) integer maxit parameter ( maxit = 30 )c ..c .. local scalars .. integer i, icompz, ii, iscale, j, jtot, k, l, l1, lend, & lendm1, lendp1, lendsv, lm1, lsv, m, mm, mm1, & nm1, nmaxit REAL*8 & anorm, b, c, eps, eps2, f, g, p, r, rt1, rt2, & s, safmax, safmin, ssfmax, ssfmin, tstc ..c .. external functions .. logical lsame REAL*8 & dlamch, dlanst, dlapy2 external lsame, dlamch, dlanst, dlapy2c ..c .. external subroutines .. external dlae2, dlaev2, dlartg, & dlascl, dlaset, dlasr, & dlasrt, dswap, xerblac ..*c .. intrinsic functions ..* intrinsic abs, max, sign, sqrt*c ..*c .. executable statements ..cc test the input parameters.c info = 0cc$$$IF( LSAME( COMPZ, 'N' ) ) THENc$$ICOMPZ = 0c$$$      ELSE IF( LSAME( COMPZ, 'V' ) ) THENc$$ICOMPZ = 1c$$$ELSE IF( LSAME( COMPZ, 'I' ) ) THENc$$ICOMPZ = 2c$$$      ELSEc$$ICOMPZ = -1c$$$END IFc$$IF( ICOMPZ.LT.0 ) THENc$$$         INFO = -1c$$ELSE IF( N.LT.0 ) THENc$$$INFO = -2c$$ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,c$$$     $N ) ) ) THENc$$INFO = -6c$$$      END IFc$$IF( INFO.NE.0 ) THENc$$$CALL XERBLA( 'SSTEQR', -INFO )c$$RETURNc$$$      END IFcc    *** New starting with version 2.5 ***c      icompz = 2c    *************************************cc     quick return if possiblec      if( n.eq.0 )     $returnc if( n.eq.1 ) then if( icompz.eq.2 ) z( 1 ) = one return end ifcc determine the unit roundoff and over/underflow thresholds.c eps = dlamch( 'e' ) eps2 = eps**2 safmin = dlamch( 's' ) safmax = one / safmin ssfmax = sqrt( safmax ) / three ssfmin = sqrt( safmin ) / eps2cc compute the eigenvalues and eigenvectors of the tridiagonalc matrix.cc$$if( icompz.eq.2 )c$$$     $call dlaset( 'full', n, n, zero, one, z, ldz )cc *** New starting with version 2.5 ***c if ( icompz .eq. 2 ) then do 5 j = 1, n-1 z(j) = zero 5 continue z( n ) = one end ifc *************************************c nmaxit = n*maxit jtot = 0cc determine where the matrix splits and choose ql or qr iterationc for each block, according to whether top or bottom diagonalc element is smaller.c l1 = 1 nm1 = n - 1c 10 continue if( l1.gt.n )$   go to 160      if( l1.gt.1 )     $e( l1-1 ) = zero if( l1.le.nm1 ) then do 20 m = l1, nm1 tst = abs( e( m ) ) if( tst.eq.zero )$         go to 30            if( tst.le.( sqrt( abs( d( m ) ) )*sqrt( abs( d( m+     $1 ) ) ) )*eps ) then e( m ) = zero go to 30 end if 20 continue end if m = nc 30 continue l = l1 lsv = l lend = m lendsv = lend l1 = m + 1 if( lend.eq.l )$   go to 10cc     scale submatrix in rows and columns l to lendc      anorm = dlanst( 'i', 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 ifcc choose between ql and qr iterationc if( abs( d( lend ) ).lt.abs( d( l ) ) ) then lend = lsv l = lendsv end ifc if( lend.gt.l ) thencc ql iterationcc look for small subdiagonal element.c 40 continue if( l.ne.lend ) then lendm1 = lend - 1 do 50 m = l, lendm1 tst = abs( e( m ) )**2 if( tst.le.( eps2*abs( d( m ) ) )*abs( d( m+1 ) )+$             safmin )go to 60   50       continue         end ifc         m = lendc   60    continue         if( m.lt.lend )     $e( m ) = zero p = d( l ) if( m.eq.l )$      go to 80cc        if remaining matrix is 2-by-2, use dlae2 or dlaev2c        to compute its eigensystem.c         if( m.eq.l+1 ) then            if( icompz.gt.0 ) then               call dlaev2( d( l ), e( l ), d( l+1 ), rt1, rt2, c, s )               work( l ) = c               work( n-1+l ) = sc$$call dlasr( 'r', 'v', 'b', n, 2, work( l ),c$$                     work( n-1+l ), z( 1, l ), ldz )cc              *** New starting with version 2.5 ***c               tst      = z(l+1)               z(l+1) = c*tst - s*z(l)               z(l)   = s*tst + c*z(l)c              *************************************            else               call dlae2( d( l ), e( l ), d( l+1 ), rt1, rt2 )            end if            d( l ) = rt1            d( l+1 ) = rt2            e( l ) = zero            l = l + 2            if( l.le.lend )     $go to 40 go to 140 end ifc if( jtot.eq.nmaxit )$      go to 140         jtot = jtot + 1cc        form shift.c         g = ( d( l+1 )-p ) / ( two*e( l ) )         r = dlapy2( g, one )         g = d( m ) - p + ( e( l ) / ( g+sign( r, g ) ) )c         s = one         c = one         p = zerocc        inner loopc         mm1 = m - 1         do 70 i = mm1, l, -1            f = s*e( i )            b = c*e( i )            call dlartg( g, f, c, s, r )            if( i.ne.m-1 )     $e( i+1 ) = r g = d( i+1 ) - p r = ( d( i )-g )*s + two*c*b p = s*r d( i+1 ) = g + p g = c*r - bcc if eigenvectors are desired, then save rotations.c if( icompz.gt.0 ) then work( i ) = c work( n-1+i ) = -s end ifc 70 continuecc if eigenvectors are desired, then apply saved rotations.c if( icompz.gt.0 ) then mm = m - l + 1c$$call dlasr( 'r', 'v', 'b', n, mm, work( l ), work( n-1+l ),c$$$     $z( 1, l ), ldz )cc *** New starting with version 2.5 ***c call dlasr( 'r', 'v', 'b', 1, mm, work( l ), & work( n-1+l ), z( l ), 1 )c ************************************* end ifc d( l ) = d( l ) - p e( l ) = g go to 40cc eigenvalue found.c 80 continue d( l ) = pc l = l + 1 if( l.le.lend )$      go to 40         go to 140c      elsecc        qr iterationcc        look for small superdiagonal element.c   90    continue         if( l.ne.lend ) then            lendp1 = lend + 1            do 100 m = l, lendp1, -1               tst = abs( e( m-1 ) )**2               if( tst.le.( eps2*abs( d( m ) ) )*abs( d( m-1 ) )+     $safmin )go to 110 100 continue end ifc m = lendc 110 continue if( m.gt.lend )$      e( m-1 ) = zero         p = d( l )         if( m.eq.l )     $go to 130cc if remaining matrix is 2-by-2, use dlae2 or dlaev2c to compute its eigensystem.c if( m.eq.l-1 ) then if( icompz.gt.0 ) then call dlaev2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2, c, s )c$$work( m ) = cc$$$               work( n-1+m ) = sc$$call dlasr( 'r', 'v', 'f', n, 2, work( m ),c$$                     work( n-1+m ), z( 1, l-1 ), ldz )cc               *** New starting with version 2.5 ***c                tst      = z(l)                z(l)   = c*tst - s*z(l-1)                z(l-1) = s*tst + c*z(l-1)c               *************************************            else               call dlae2( d( l-1 ), e( l-1 ), d( l ), rt1, rt2 )            end if            d( l-1 ) = rt1            d( l ) = rt2            e( l-1 ) = zero            l = l - 2            if( l.ge.lend )     $go to 90 go to 140 end ifc if( jtot.eq.nmaxit )$      go to 140         jtot = jtot + 1cc        form shift.c         g = ( d( l-1 )-p ) / ( two*e( l-1 ) )         r = dlapy2( g, one )         g = d( m ) - p + ( e( l-1 ) / ( g+sign( r, g ) ) )c         s = one         c = one         p = zerocc        inner loopc         lm1 = l - 1         do 120 i = m, lm1            f = s*e( i )            b = c*e( i )            call dlartg( g, f, c, s, r )            if( i.ne.m )     $e( i-1 ) = r g = d( i ) - p r = ( d( i+1 )-g )*s + two*c*b p = s*r d( i ) = g + p g = c*r - bcc if eigenvectors are desired, then save rotations.c if( icompz.gt.0 ) then work( i ) = c work( n-1+i ) = s end ifc 120 continuecc if eigenvectors are desired, then apply saved rotations.c if( icompz.gt.0 ) then mm = l - m + 1c$$call dlasr( 'r', 'v', 'f', n, mm, work( m ), work( n-1+m ),c$$$     $z( 1, m ), ldz )cc *** New starting with version 2.5 ***c call dlasr( 'r', 'v', 'f', 1, mm, work( m ), work( n-1+m ), & z( m ), 1 )c ************************************* end ifc d( l ) = d( l ) - p e( lm1 ) = g go to 90cc eigenvalue found.c 130 continue d( l ) = pc l = l - 1 if( l.ge.lend )$      go to 90         go to 140c      end ifcc     undo scaling if necessaryc  140 continue      if( iscale.eq.1 ) then         call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv+1, 1,     $d( lsv ), n, info ) call dlascl( 'g', 0, 0, ssfmax, anorm, lendsv-lsv, 1, e( lsv ),$                n, info )      else if( iscale.eq.2 ) then         call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv+1, 1,     $d( lsv ), n, info ) call dlascl( 'g', 0, 0, ssfmin, anorm, lendsv-lsv, 1, e( lsv ),$                n, info )      end ifcc     check for no convergence to an eigenvalue after a totalc     of n*maxit iterations.c      if( jtot.lt.nmaxit )     $go to 10 do 150 i = 1, n - 1 if( e( i ).ne.zero )$      info = info + 1  150 continue      go to 190cc     order eigenvalues and eigenvectors.c  160 continue      if( icompz.eq.0 ) thencc        use quick sortc         call dlasrt( 'i', n, d, info )c      elsecc        use selection sort to minimize swaps of eigenvectorsc         do 180 ii = 2, n            i = ii - 1            k = i            p = d( i )            do 170 j = ii, n               if( d( j ).lt.p ) then                  k = j                  p = d( j )               end if  170       continue            if( k.ne.i ) then               d( k ) = d( i )               d( i ) = pc\$               call dswap( n, z( 1, i ), 1, z( 1, k ), 1 )c           *** New starting with version 2.5 ***c               p    = z(k)               z(k) = z(i)               z(i) = pc           *************************************            end if  180    continue      end ifc  190 continue      returncc     %---------------%c     | End of dstqrb |c     %---------------%c      end   

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