diff libcruft/arpack/src/zneigh.f @ 12274:9f5d2ef078e8 release-3-4-x

import ARPACK sources to libcruft from Debian package libarpack2 2.1+parpack96.dfsg-3+b1
author John W. Eaton <jwe@octave.org>
date Fri, 28 Jan 2011 14:04:33 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/arpack/src/zneigh.f	Fri Jan 28 14:04:33 2011 -0500
@@ -0,0 +1,257 @@
+c\BeginDoc
+c
+c\Name: zneigh
+c
+c\Description:
+c  Compute the eigenvalues of the current upper Hessenberg matrix
+c  and the corresponding Ritz estimates given the current residual norm.
+c
+c\Usage:
+c  call zneigh
+c     ( RNORM, N, H, LDH, RITZ, BOUNDS, Q, LDQ, WORKL, RWORK, IERR )
+c
+c\Arguments
+c  RNORM   Double precision scalar.  (INPUT)
+c          Residual norm corresponding to the current upper Hessenberg 
+c          matrix H.
+c
+c  N       Integer.  (INPUT)
+c          Size of the matrix H.
+c
+c  H       Complex*16 N by N array.  (INPUT)
+c          H contains the current upper Hessenberg matrix.
+c
+c  LDH     Integer.  (INPUT)
+c          Leading dimension of H exactly as declared in the calling
+c          program.
+c
+c  RITZ    Complex*16 array of length N.  (OUTPUT)
+c          On output, RITZ(1:N) contains the eigenvalues of H.
+c
+c  BOUNDS  Complex*16 array of length N.  (OUTPUT)
+c          On output, BOUNDS contains the Ritz estimates associated with
+c          the eigenvalues held in RITZ.  This is equal to RNORM 
+c          times the last components of the eigenvectors corresponding 
+c          to the eigenvalues in RITZ.
+c
+c  Q       Complex*16 N by N array.  (WORKSPACE)
+c          Workspace needed to store the eigenvectors of H.
+c
+c  LDQ     Integer.  (INPUT)
+c          Leading dimension of Q exactly as declared in the calling
+c          program.
+c
+c  WORKL   Complex*16 work array of length N**2 + 3*N.  (WORKSPACE)
+c          Private (replicated) array on each PE or array allocated on
+c          the front end.  This is needed to keep the full Schur form
+c          of H and also in the calculation of the eigenvectors of H.
+c
+c  RWORK   Double precision  work array of length N (WORKSPACE)
+c          Private (replicated) array on each PE or array allocated on
+c          the front end. 
+c
+c  IERR    Integer.  (OUTPUT)
+c          Error exit flag from zlahqr or ztrevc.
+c
+c\EndDoc
+c
+c-----------------------------------------------------------------------
+c
+c\BeginLib
+c
+c\Local variables:
+c     xxxxxx  Complex*16
+c
+c\Routines called:
+c     ivout   ARPACK utility routine that prints integers.
+c     arscnd  ARPACK utility routine for timing.
+c     zmout   ARPACK utility routine that prints matrices
+c     zvout   ARPACK utility routine that prints vectors.
+c     dvout   ARPACK utility routine that prints vectors.
+c     zlacpy  LAPACK matrix copy routine.
+c     zlahqr  LAPACK routine to compute the Schur form of an
+c             upper Hessenberg matrix.
+c     zlaset  LAPACK matrix initialization routine.
+c     ztrevc  LAPACK routine to compute the eigenvectors of a matrix
+c             in upper triangular form
+c     zcopy   Level 1 BLAS that copies one vector to another. 
+c     zdscal  Level 1 BLAS that scales a complex vector by a real number.
+c     dznrm2  Level 1 BLAS that computes the norm of a vector.
+c     
+c
+c\Author
+c     Danny Sorensen               Phuong Vu
+c     Richard Lehoucq              CRPC / Rice University
+c     Dept. of Computational &     Houston, Texas
+c     Applied Mathematics 
+c     Rice University           
+c     Houston, Texas 
+c
+c\SCCS Information: @(#)
+c FILE: neigh.F   SID: 2.2   DATE OF SID: 4/20/96   RELEASE: 2
+c
+c\Remarks
+c     None
+c
+c\EndLib
+c
+c-----------------------------------------------------------------------
+c
+      subroutine zneigh (rnorm, n, h, ldh, ritz, bounds, 
+     &                   q, ldq, workl, rwork, ierr)
+c
+c     %----------------------------------------------------%
+c     | Include files for debugging and timing information |
+c     %----------------------------------------------------%
+c
+      include   'debug.h'
+      include   'stat.h'
+c
+c     %------------------%
+c     | Scalar Arguments |
+c     %------------------%
+c
+      integer    ierr, n, ldh, ldq
+      Double precision     
+     &           rnorm
+c
+c     %-----------------%
+c     | Array Arguments |
+c     %-----------------%
+c
+      Complex*16     
+     &           bounds(n), h(ldh,n), q(ldq,n), ritz(n),
+     &           workl(n*(n+3)) 
+      Double precision 
+     &           rwork(n)
+c 
+c     %------------%
+c     | Parameters |
+c     %------------%
+c
+      Complex*16     
+     &           one, zero
+      Double precision
+     &           rone
+      parameter  (one = (1.0D+0, 0.0D+0), zero = (0.0D+0, 0.0D+0),
+     &           rone = 1.0D+0)
+c 
+c     %------------------------%
+c     | Local Scalars & Arrays |
+c     %------------------------%
+c
+      logical    select(1)
+      integer    j,  msglvl
+      Complex*16     
+     &           vl(1)
+      Double precision
+     &           temp
+c
+c     %----------------------%
+c     | External Subroutines |
+c     %----------------------%
+c
+      external   zlacpy, zlahqr, ztrevc, zcopy, 
+     &           zdscal, zmout, zvout, arscnd
+c
+c     %--------------------%
+c     | External Functions |
+c     %--------------------%
+c
+      Double precision 
+     &           dznrm2
+      external   dznrm2
+c
+c     %-----------------------%
+c     | Executable Statements |
+c     %-----------------------%
+c
+c     %-------------------------------%
+c     | Initialize timing statistics  |
+c     | & message level for debugging |
+c     %-------------------------------%
+c
+      call arscnd (t0)
+      msglvl = mceigh
+c 
+      if (msglvl .gt. 2) then
+          call zmout (logfil, n, n, h, ldh, ndigit, 
+     &         '_neigh: Entering upper Hessenberg matrix H ')
+      end if
+c 
+c     %----------------------------------------------------------%
+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     |    zlahqr returns the full Schur form of H               | 
+c     |    in WORKL(1:N**2), and the Schur vectors in q.         |
+c     %----------------------------------------------------------%
+c
+      call zlacpy ('All', n, n, h, ldh, workl, n)
+      call zlaset ('All', n, n, zero, one, q, ldq)
+      call zlahqr (.true., .true., n, 1, n, workl, ldh, ritz,
+     &             1, n, q, ldq, ierr)
+      if (ierr .ne. 0) go to 9000
+c
+      call zcopy (n, q(n-1,1), ldq, bounds, 1)
+      if (msglvl .gt. 1) then
+         call zvout (logfil, n, bounds, ndigit,
+     &              '_neigh: last row of the Schur matrix for H')
+      end if
+c
+c     %----------------------------------------------------------%
+c     | 2. Compute the eigenvectors of the full Schur form T and |
+c     |    apply the Schur vectors to get the corresponding      |
+c     |    eigenvectors.                                         |
+c     %----------------------------------------------------------%
+c
+      call ztrevc ('Right', 'Back', select, n, workl, n, vl, n, q, 
+     &             ldq, n, n, workl(n*n+1), rwork, ierr)
+c
+      if (ierr .ne. 0) go to 9000
+c
+c     %------------------------------------------------%
+c     | Scale the returning eigenvectors so that their |
+c     | Euclidean norms are all one. LAPACK subroutine |
+c     | ztrevc 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
+      do 10 j=1, n
+            temp = dznrm2( n, q(1,j), 1 )
+            call zdscal ( n, rone / temp, q(1,j), 1 )
+   10 continue
+c
+      if (msglvl .gt. 1) then
+         call zcopy(n, q(n,1), ldq, workl, 1)
+         call zvout (logfil, n, workl, ndigit,
+     &              '_neigh: Last row of the eigenvector matrix for H')
+      end if
+c
+c     %----------------------------%
+c     | Compute the Ritz estimates |
+c     %----------------------------%
+c
+      call zcopy(n, q(n,1), n, bounds, 1)
+      call zdscal(n, rnorm, bounds, 1)
+c
+      if (msglvl .gt. 2) then
+         call zvout (logfil, n, ritz, ndigit,
+     &              '_neigh: The eigenvalues of H')
+         call zvout (logfil, n, bounds, ndigit,
+     &              '_neigh: Ritz estimates for the eigenvalues of H')
+      end if
+c
+      call arscnd(t1)
+      tceigh = tceigh + (t1 - t0)
+c
+ 9000 continue
+      return
+c
+c     %---------------%
+c     | End of zneigh |
+c     %---------------%
+c
+      end