Numérotation des lignes :

dneigh
C DNEIGH    SOURCE    FANDEUR   22/05/02    21:15:11     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dneighcc\Description:c  Compute the eigenvalues of the current upper Hessenberg matrixc  and the corresponding Ritz estimates given the current residual norm.cc\Usage:c  call dneigh(RNORM, N, H, LDH, RITZR, RITZI, BOUNDS,Q ,LDQ,WORKL,IERR)cc\Argumentsc  RNORM   Double precision scalar.  (INPUT)c          Residual norm corresponding to the current upper Hessenbergc          matrix H.cc  N       Integer.  (INPUT)c          Size of the matrix H.cc  H       REAL*8 N by N array.  (INPUT)c          H contains the current upper Hessenberg matrix.cc  LDH     Integer.  (INPUT)c          Leading dimension of H exactly as declared in the callingc          program.cc  RITZR,  Double precision arrays of length N.  (OUTPUT)c  RITZI   On output, RITZR(1:N) (resp. RITZI(1:N)) contains the realc          (respectively imaginary) parts of the eigenvalues of H.cc  BOUNDS  Double precision array of length N.  (OUTPUT)c c          On output, BOUNDS contains the Ritz estimates associated withc          the eigenvalues RITZR and RITZI.  This is equal to RNORMc          times the last components of the eigenvectors correspondingc          to the eigenvalues in RITZR and RITZI.cc  Q       REAL*8 N by N array.  (WORKSPACE)c          Workspace needed to store the eigenvectors of H.cc  LDQ     Integer.  (INPUT)c          Leading dimension of Q exactly as declared in the callingc          program.cc  WORKL   Double precision work array of length N**2 + 3*N.  (WORKSPACE)c          Private (replicated) array on each PE or array allocated onc          the front end.  This is needed to keep the full Schur formc          of H and also in the calculation of the eigenvectors of H.cc  IERR    Integer.  (OUTPUT)c          Error exit flag from dlaqrb or dtrevc.cc\EndDoccc-----------------------------------------------------------------------cc\BeginLibcc\Local variables:c     xxxxxx  realcc\Routines called:c     dlaqrb  ARPACK routine to compute the real Schur form of anc             upper Hessenberg matrix and last row of the Schur vectors.c     arscnd  ARPACK utility routine for timing.c     dmout   ARPACK utility routine that prints matricesc     dvout   ARPACK utility routine that prints vectors.c     dlacpy  LAPACK matrix copy routine.c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.c     dtrevc  LAPACK routine to compute the eigenvectors of a matrixc             in upper quasi-triangular formc     dgemv   Level 2 BLAS routine for matrix vector multiplication.c     dcopy   Level 1 BLAS that copies one vector to another .c     dnrm2   Level 1 BLAS that computes the norm of a vector.c     dscal   Level 1 BLAS that scales a vector.ccc\Authorc     Danny Sorensen               Phuong Vuc     Richard Lehoucq              CRPC / Rice Universityc     Dept. of Computational &     Houston, Texasc     Applied Mathematicsc     Rice Universityc     Houston, Texascc\Revision history:c     xx/xx/92: Version ' 2.1'cc\SCCS Information: @(#)c FILE: neigh.F   SID: 2.3   DATE OF SID: 4/20/96   RELEASE: 2cc\Remarksc     Nonecc\EndLibcc-----------------------------------------------------------------------c      subroutine dneigh (rnorm, n, h, ldh, ritzr, ritzi, bounds,     &                   q, ldq, workl, ierr)cc     %----------------------------------------------------%c     | Include files for debugging and timing information |c -INC TARTRAKc     %----------------------------------------------------%ccc     %------------------%c     | Scalar Arguments |c     %------------------%c      integer    ierr, n, ldh, ldq      REAL*8     &           rnormcc     %-----------------%c     | Array Arguments |c     %-----------------%c      REAL*8     &           bounds(n), h(ldh,n), q(ldq,n), ritzi(n), ritzr(n),     &           workl(n*(n+3))cc     %------------%c     | Parameters |c     %------------%c      REAL*8     &           one, zero      parameter (one = 1.0D+0, zero = 0.0D+0)cc     %------------------------%c     | Local Scalars & Arrays |c     %------------------------%c      logical    select(1)      integer    i, iconj, msglvl      REAL*8     &           temp, vl(1)      parameter (msglvl=0)cc     %----------------------%c     | External Subroutines |c     %----------------------%c      external   dcopy, dlacpy, dlaqrb, dtrevc,     &           dvout, arscndcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           dlapy2, dnrm2      external   dlapy2, dnrm2cc     %---------------------%**c     | Intrinsic Functions |**c     %---------------------%**c**      intrinsic  abs**c**c     %-----------------------%**c     | Executable Statements |c     %-----------------------%ccc     %-------------------------------%c     | Initialize timing statistics  |c     | & message level for debugging |c     %-------------------------------%c*      call arscnd (t0)  c       msglvl = mneighcc       if (msglvl .gt. 2) thenc           call dmout (logfil, n, n, h, ldh, ndigit,c      &         '_neigh: Entering upper Hessenberg matrix H ')c       end ifcc     %-----------------------------------------------------------%c     | 1. Compute the eigenvalues, the last components of the    |c     |    corresponding Schur vectors and the full Schur form T  |c     |    of the current upper Hessenberg matrix H.              |c     | dlaqrb returns the full Schur form of H in WORKL(1:N**2)  |c     | and the last components of the Schur vectors in BOUNDS.   |c     %-----------------------------------------------------------%c      call dlacpy ('All', n, n, h, ldh, workl, n)       do 5 j = 1, n-1        bounds(j) = zero5     continue      bounds(n) = one       call dlahqr(.true., .true., n, 1, n, workl, n, ritzr, ritzi, 1,     &                                        1, bounds, 1, ierr) c       call dlaqrb (.true., n, 1, n, workl, n, ritzr, ritzi, bounds,c      &             ierr)        if (ierr .ne. 0) go to 9000cc       if (msglvl .gt. 1) thenc          call dvout ( n, bounds, ndigit,c      &              '_neigh: last row of the Schur matrix for H')c       end ifcc     %-----------------------------------------------------------%c     | 2. Compute the eigenvectors of the full Schur form T and  |c     |    apply the last components of the Schur vectors to get  |c     |    the last components of the corresponding eigenvectors. |c     | Remember that if the i-th and (i+1)-st eigenvalues are    |c     | complex conjugate pairs, then the real & imaginary part   |c     | of the eigenvector components are split across adjacent   |c     | columns of Q.                                             |c     %-----------------------------------------------------------%c      call dtrevc ('R', 'A', select, n, workl, n, vl, n, q, ldq,     &             n, n, workl(n*n+1), ierr)c      if (ierr .ne. 0) go to 9000cc     %------------------------------------------------%c     | Scale the returning eigenvectors so that their |c     | euclidean norms are all one. LAPACK subroutine |c     | dtrevc returns each eigenvector normalized so  |c     | that the element of largest magnitude has      |c     | magnitude 1; here the magnitude of a complex   |c     | number (x,y) is taken to be |x| + |y|.         |c     %------------------------------------------------%c      iconj = 0      do 10 i=1, n         if ( abs( ritzi(i) ) .le. zero ) thencc           %----------------------%c           | Real eigenvalue case |c           %----------------------%c            temp = dnrm2( n, q(1,i), 1 )            call dscal ( n, one / temp, q(1,i), 1 )         elsecc           %-------------------------------------------%c           | Complex conjugate pair case. Note that    |c           | since the real and imaginary part of      |c           | the eigenvector are stored in consecutive |c           | columns, we further normalize by the      |c           | square root of two.                       |c           %-------------------------------------------%c            if (iconj .eq. 0) then               temp = dlapy2( dnrm2( n, q(1,i), 1 ),     &                        dnrm2( n, q(1,i+1), 1 ) )               call dscal ( n, one / temp, q(1,i), 1 )               call dscal ( n, one / temp, q(1,i+1), 1 )               iconj = 1            else               iconj = 0            end if         end if   10 continuec      call dgemv ('T', n, n, one, q, ldq, bounds, 1, zero, workl, 1)cc       if (msglvl .gt. 1) thenc          call dvout ( n, workl, ndigit,c      &              '_neigh: Last row of the eigenvector matrix for H')c       end ifcc     %----------------------------%c     | Compute the Ritz estimates |c     %----------------------------%c      iconj = 0      do 20 i = 1, n         if ( abs( ritzi(i) ) .le. zero ) thencc           %----------------------%c           | Real eigenvalue case |c           %----------------------%c            bounds(i) = rnorm * abs( workl(i) )         elsecc           %-------------------------------------------%c           | Complex conjugate pair case. Note that    |c           | since the real and imaginary part of      |c           | the eigenvector are stored in consecutive |c           | columns, we need to take the magnitude    |c           | of the last components of the two vectors |c           %-------------------------------------------%c            if (iconj .eq. 0) then               bounds(i) = rnorm * dlapy2( workl(i), workl(i+1) )               bounds(i+1) = bounds(i)               iconj = 1            else               iconj = 0            end if         end if   20 continuecc       if (msglvl .gt. 2) thenc          call dvout ( n, ritzr, ndigit,c      &              '_neigh: Real part of the eigenvalues of H')c          call dvout ( n, ritzi, ndigit,c      &              '_neigh: Imaginary part of the eigenvalues of H')c          call dvout ( n, bounds, ndigit,c      &              '_neigh: Ritz estimates for the eigenvalues of H')c       end ifc*      call arscnd (t1)*      tneigh = tneigh + (t1 - t0)c 9000 continue      returncc     %---------------%c     | End of dneigh |c     %---------------%c      end    

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