Sources Esope | Cast3M

dnapps
C DNAPPS    SOURCE    FANDEUR   22/05/02    21:15:10     11359          
c-----------------------------------------------------------------------
c\BeginDoc
c
c\Name: dnapps
c
c\Description:
c  Given the Arnoldi factorization
c
c     A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,
c
c  apply NP implicit shifts resulting in
c
c     A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q
c
c  where Q is an orthogonal matrix which is the product of rotations
c  and reflections resulting from the NP bulge chage sweeps.
c  The updated Arnoldi factorization becomes:
c
c     A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.
c
c\Usage:
c  call dnapps
c     ( N, KEV, NP, SHIFTR, SHIFTI, V, LDV, H, LDH, RESID, Q, LDQ,
c       WORKL, WORKD )
c
c\Arguments
c  N       Integer.  (INPUT)
c          Problem size, i.e. size of matrix A.
c
c  KEV     Integer.  (INPUT/OUTPUT)
c          KEV+NP is the size of the input matrix H.
c          KEV is the size of the updated matrix HNEW.  KEV is only
c          updated on ouput when fewer than NP shifts are applied in
c          order to keep the conjugate pair together.
c
c  NP      Integer.  (INPUT)
c          Number of implicit shifts to be applied.
c
c  SHIFTR, Double precision array of length NP.  (INPUT)
c  SHIFTI  Real and imaginary part of the shifts to be applied.
c          Upon, entry to dnapps, the shifts must be sorted so that the
c          conjugate pairs are in consecutive locations.
c
c  V       REAL*8 N by (KEV+NP) array.  (INPUT/OUTPUT)
c          On INPUT, V contains the current KEV+NP Arnoldi vectors.
c          On OUTPUT, V contains the updated KEV Arnoldi vectors
c          in the first KEV columns of V.
c
c  LDV     Integer.  (INPUT)
c          Leading dimension of V exactly as declared in the calling
c          program.
c
c  H       REAL*8 (KEV+NP) by (KEV+NP) array.  (INPUT/OUTPUT)
c          On INPUT, H contains the current KEV+NP by KEV+NP upper
c          Hessenber matrix of the Arnoldi factorization.
c          On OUTPUT, H contains the updated KEV by KEV upper Hessenberg
c          matrix in the KEV leading submatrix.
c
c  LDH     Integer.  (INPUT)
c          Leading dimension of H exactly as declared in the calling
c          program.
c
c  RESID   Double precision array of length N.  (INPUT/OUTPUT)
c          On INPUT, RESID contains the the residual vector r_{k+p}.
c          On OUTPUT, RESID is the update residual vector rnew_{k}
c          in the first KEV locations.
c
c  Q       REAL*8 KEV+NP by KEV+NP work array.  (WORKSPACE)
c          Work array used to accumulate the rotations and reflections
c          during the bulge chase sweep.
c
c  LDQ     Integer.  (INPUT)
c          Leading dimension of Q exactly as declared in the calling
c          program.
c
c  WORKL   Double precision work array of length (KEV+NP).  (WORKSPACE)
c          Private (replicated) array on each PE or array allocated on
c          the front end.
c
c  WORKD   Double precision work array of length 2*N.  (WORKSPACE)
c          Distributed array used in the application of the accumulated
c          orthogonal matrix Q.
c
c\EndDoc
c
c-----------------------------------------------------------------------
c
c\BeginLib
c
c\Local variables:
c     xxxxxx  real
c
c\References:
c  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
c     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
c     pp 357-385.
c
c\Routines called:
c     ivout   ARPACK utility routine that prints integers.
c     arscnd  ARPACK utility routine for timing. -> deleted by BP in 2020
c     dmout   ARPACK utility routine that prints matrices.
c     dvout   ARPACK utility routine that prints vectors.
c     dlabad  LAPACK routine that computes machine constants.
c     dlacpy  LAPACK matrix copy routine.
c     dlamch  LAPACK routine that determines machine constants.
c     dlanhs  LAPACK routine that computes various norms of a matrix.
c     dlapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
c     dlarf   LAPACK routine that applies Householder reflection to
c             a matrix.
c     dlarfg  LAPACK Householder reflection construction routine.
c     dlartg  LAPACK Givens rotation construction routine.
c     dlaset  LAPACK matrix initialization routine.
c     dgemv   Level 2 BLAS routine for matrix vector multiplication.
c     daxpy   Level 1 BLAS that computes a vector triad.
c     dcopy   Level 1 BLAS that copies one vector to another .
c     dscal   Level 1 BLAS that scales a vector.
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\Revision history:
c     xx/xx/92: Version ' 2.4'
c
c\SCCS Information: @(#)
c FILE: napps.F   SID: 2.4   DATE OF SID: 3/28/97   RELEASE: 2
c
c\Remarks
c  1. In this version, each shift is applied to all the sublocks of
c     the Hessenberg matrix H and not just to the submatrix that it
c     comes from. Deflation as in LAPACK routine dlahqr (QR algorithm
c     for upper Hessenberg matrices ) is used.
c     The subdiagonals of H are enforced to be non-negative.
c
c\EndLib
c
c-----------------------------------------------------------------------
c
      subroutine dnapps
     &   ( n, kev, np, shiftr, shifti, v, ldv, h, ldh, resid, q, ldq,
     &     workl, workd )
c
c     %----------------------------------------------------%
c     | Include files for debugging and timing information |
c -INC TARTRAK
c     %----------------------------------------------------%
c
c
c     %------------------%
c     | Scalar Arguments |
c     %------------------%
c
      integer    kev, ldh, ldq, ldv, n, np
c
c     %-----------------%
c     | Array Arguments |
c     %-----------------%
c
      REAL*8
     &           h(ldh,kev+np), resid(n), shifti(np), shiftr(np),
     &           v(ldv,kev+np), q(ldq,kev+np), workd(2*n), workl(kev+np)
c
c     %------------%
c     | Parameters |
c     %------------%
c
      REAL*8
     &           one, zero
      parameter (one = 1.0D+0, zero = 0.0D+0)
c
c     %------------------------%
c     | Local Scalars & Arrays |
c     %------------------------%
c
      integer    i, iend, ir, istart, j, jj, kplusp, msglvl, nr
      logical    cconj, first
      REAL*8
     &           c, f, g, h11, h12, h21, h22, h32, ovfl, r, s, sigmai,
     &           sigmar, smlnum, ulp, unfl, u(3), t, tau, tst1
      save       first, ovfl, smlnum, ulp, unfl
      parameter (msglvl=0)
c
c     %----------------------%
c     | External Subroutines |
c     %----------------------%
c
      external   daxpy, dcopy, dscal,
     &           dlacpy, dlarfg, dlarf,
     &           dlaset, dlabad, dlartg
c
c     %--------------------%
c     | External Functions |
c     %--------------------%
c
      REAL*8
     &           dlamch, dlanhs, dlapy2
      external   dlamch, dlanhs, dlapy2
c
c     %----------------------%
**c     | Intrinsics Functions |
**c     %----------------------%
**c
**      intrinsic  abs, max, min
**c
**c     %----------------%
**c     | Data statments |
**c     %----------------%
**c
      data       first / .true. /
**c
**c     %-----------------------%
**c     | Executable Statements |
c     %-----------------------%
c
      if (first) then
c
c        %-----------------------------------------------%
c        | Set machine-dependent constants for the       |
c        | stopping criterion. If norm(H) &lt;= sqrt(OVFL), |
c        | overflow should not occur.                    |
c        | REFERENCE: LAPACK subroutine dlahqr           |
c        %-----------------------------------------------%
c
         unfl = dlamch( 'safe minimum' )
         ovfl = one / unfl
         call dlabad( unfl, ovfl )
         ulp = dlamch( 'precision' )
         smlnum = unfl*( n / ulp )
         first = .false.
      end if
c
c     %-------------------------------%
c     | Initialize timing statistics  |
c     | & message level for debugging |
c     %-------------------------------%
c
*      call arscnd (t0)
c       msglvl = mnapps
      kplusp = kev + np
c
c     %--------------------------------------------%
c     | Initialize Q to the identity to accumulate |
c     | the rotations and reflections              |
c     %--------------------------------------------%
c
      call dlaset ('All', kplusp, kplusp, zero, one, q, ldq)
c
c     %----------------------------------------------%
c     | Quick return if there are no shifts to apply |
c     %----------------------------------------------%
c
      if (np .eq. 0) go to 9000
c
c     %----------------------------------------------%
c     | Chase the bulge with the application of each |
c     | implicit shift. Each shift is applied to the |
c     | whole matrix including each block.           |
c     %----------------------------------------------%
c
      cconj = .false.
      do 110 jj = 1, np
         sigmar = shiftr(jj)
         sigmai = shifti(jj)
c
c          if (msglvl .gt. 2 ) then
c             call ivout ( 1, jj, ndigit,
c      &               '_napps: shift number.')
c             call dvout ( 1, sigmar, ndigit,
c      &               '_napps: The real part of the shift ')
c             call dvout ( 1, sigmai, ndigit,
c      &               '_napps: The imaginary part of the shift ')
c          end if
c
c        %-------------------------------------------------%
c        | The following set of conditionals is necessary  |
c        | in order that complex conjugate pairs of shifts |
c        | are applied together or not at all.             |
c        %-------------------------------------------------%
c
         if ( cconj ) then
c
c           %-----------------------------------------%
c           | cconj = .true. means the previous shift |
c           | had non-zero imaginary part.            |
c           %-----------------------------------------%
c
            cconj = .false.
            go to 110
         else if ( jj .lt. np .and. abs( sigmai ) .gt. zero ) then
c
c           %------------------------------------%
c           | Start of a complex conjugate pair. |
c           %------------------------------------%
c
            cconj = .true.
         else if ( jj .eq. np .and. abs( sigmai ) .gt. zero ) then
c
c           %----------------------------------------------%
c           | The last shift has a nonzero imaginary part. |
c           | Don't apply it; thus the order of the        |
c           | compressed H is order KEV+1 since only np-1  |
c           | were applied.                                |
c           %----------------------------------------------%
c
            kev = kev + 1
            go to 110
         end if
         istart = 1
   20    continue
c
c        %--------------------------------------------------%
c        | if sigmai = 0 then                               |
c        |    Apply the jj-th shift ...                     |
c        | else                                             |
c        |    Apply the jj-th and (jj+1)-th together ...    |
c        |    (Note that jj &lt; np at this point in the code) |
c        | end                                              |
c        | to the current block of H. The next do loop      |
c        | determines the current block ;                   |
c        %--------------------------------------------------%
c
         do 30 i = istart, kplusp-1
c
c           %----------------------------------------%
c           | Check for splitting and deflation. Use |
c           | a standard test as in the QR algorithm |
c           | REFERENCE: LAPACK subroutine dlahqr    |
c           %----------------------------------------%
c
            tst1 = abs( h( i, i ) ) + abs( h( i+1, i+1 ) )
            if( tst1.eq.zero )
     &         tst1 = dlanhs( '1', kplusp-jj+1, h, ldh, workl )
            if( abs( h( i+1,i ) ).le.max( ulp*tst1, smlnum ) ) then
               if (msglvl .gt. 0) then
                  call ivout ( 1, i, ndigit,
     &                 '_napps: matrix splitting at row/column no.')
                  call ivout ( 1, jj, ndigit,
     &                 '_napps: matrix splitting with shift number.')
                  call dvout ( 1, h(i+1,i), ndigit,
     &                 '_napps: off diagonal element.')
               end if
               iend = i
               h(i+1,i) = zero
               go to 40
            end if
   30    continue
         iend = kplusp
   40    continue
c
c          if (msglvl .gt. 2) then
c              call ivout ( 1, istart, ndigit,
c      &                   '_napps: Start of current block ')
c              call ivout ( 1, iend, ndigit,
c      &                   '_napps: End of current block ')
c          end if
c
c        %------------------------------------------------%
c        | No reason to apply a shift to block of order 1 |
c        %------------------------------------------------%
c
         if ( istart .eq. iend ) go to 100
c
c        %------------------------------------------------------%
c        | If istart + 1 = iend then no reason to apply a       |
c        | complex conjugate pair of shifts on a 2 by 2 matrix. |
c        %------------------------------------------------------%
c
         if ( istart + 1 .eq. iend .and. abs( sigmai ) .gt. zero )
     &      go to 100
c
         h11 = h(istart,istart)
         h21 = h(istart+1,istart)
         if ( abs( sigmai ) .le. zero ) then
c
c           %---------------------------------------------%
c           | Real-valued shift ==> apply single shift QR |
c           %---------------------------------------------%
c
            f = h11 - sigmar
            g = h21
c
            do 80 i = istart, iend-1
c
c              %-----------------------------------------------------%
c              | Contruct the plane rotation G to zero out the bulge |
c              %-----------------------------------------------------%
c
               call dlartg (f, g, c, s, r)
               if (i .gt. istart) then
c
c                 %-------------------------------------------%
c                 | The following ensures that h(1:iend-1,1), |
c                 | the first iend-2 off diagonal of elements |
c                 | H, remain non negative.                   |
c                 %-------------------------------------------%
c
                  if (r .lt. zero) then
                     r = -r
                     c = -c
                     s = -s
                  end if
                  h(i,i-1) = r
                  h(i+1,i-1) = zero
               end if
c
c              %---------------------------------------------%
c              | Apply rotation to the left of H;  H &lt;- G'*H |
c              %---------------------------------------------%
c
               do 50 j = i, kplusp
                  t        =  c*h(i,j) + s*h(i+1,j)
                  h(i+1,j) = -s*h(i,j) + c*h(i+1,j)
                  h(i,j)   = t
   50          continue
c
c              %---------------------------------------------%
c              | Apply rotation to the right of H;  H &lt;- H*G |
c              %---------------------------------------------%
c
               do 60 j = 1, min(i+2,iend)
                  t        =  c*h(j,i) + s*h(j,i+1)
                  h(j,i+1) = -s*h(j,i) + c*h(j,i+1)
                  h(j,i)   = t
   60          continue
c
c              %----------------------------------------------------%
c              | Accumulate the rotation in the matrix Q;  Q &lt;- Q*G |
c              %----------------------------------------------------%
c
               do 70 j = 1, min( i+jj, kplusp )
                  t        =   c*q(j,i) + s*q(j,i+1)
                  q(j,i+1) = - s*q(j,i) + c*q(j,i+1)
                  q(j,i)   = t
   70          continue
c
c              %---------------------------%
c              | Prepare for next rotation |
c              %---------------------------%
c
               if (i .lt. iend-1) then
                  f = h(i+1,i)
                  g = h(i+2,i)
               end if
   80       continue
c
c           %-----------------------------------%
c           | Finished applying the real shift. |
c           %-----------------------------------%
c
         else
c
c           %----------------------------------------------------%
c           | Complex conjugate shifts ==> apply double shift QR |
c           %----------------------------------------------------%
c
            h12 = h(istart,istart+1)
            h22 = h(istart+1,istart+1)
            h32 = h(istart+2,istart+1)
c
c           %---------------------------------------------------------%
c           | Compute 1st column of (H - shift*I)*(H - conj(shift)*I) |
c           %---------------------------------------------------------%
c
            s    = 2.0*sigmar
            t = dlapy2 ( sigmar, sigmai )
            u(1) = ( h11 * (h11 - s) + t * t ) / h21 + h12
            u(2) = h11 + h22 - s
            u(3) = h32
c
            do 90 i = istart, iend-1
c
               nr = min ( 3, iend-i+1 )
c
c              %-----------------------------------------------------%
c              | Construct Householder reflector G to zero out u(1). |
c              | G is of the form I - tau*( 1 u )' * ( 1 u' ).       |
c              %-----------------------------------------------------%
c
               call dlarfg ( nr, u(1), u(2), 1, tau )
c
               if (i .gt. istart) then
                  h(i,i-1)   = u(1)
                  h(i+1,i-1) = zero
                  if (i .lt. iend-1) h(i+2,i-1) = zero
               end if
               u(1) = one
c
c              %--------------------------------------%
c              | Apply the reflector to the left of H |
c              %--------------------------------------%
c
               call dlarf ('Left', nr, kplusp-i+1, u, 1, tau,
     &                     h(i,i), ldh, workl)
c
c              %---------------------------------------%
c              | Apply the reflector to the right of H |
c              %---------------------------------------%
c
               ir = min ( i+3, iend )
               call dlarf ('Right', ir, nr, u, 1, tau,
     &                     h(1,i), ldh, workl)
c
c              %-----------------------------------------------------%
c              | Accumulate the reflector in the matrix Q;  Q &lt;- Q*G |
c              %-----------------------------------------------------%
c
               call dlarf ('Right', kplusp, nr, u, 1, tau,
     &                     q(1,i), ldq, workl)
c
c              %----------------------------%
c              | Prepare for next reflector |
c              %----------------------------%
c
               if (i .lt. iend-1) then
                  u(1) = h(i+1,i)
                  u(2) = h(i+2,i)
                  if (i .lt. iend-2) u(3) = h(i+3,i)
               end if
c
   90       continue
c
c           %--------------------------------------------%
c           | Finished applying a complex pair of shifts |
c           | to the current block                       |
c           %--------------------------------------------%
c
         end if
c
  100    continue
c
c        %---------------------------------------------------------%
c        | Apply the same shift to the next block if there is any. |
c        %---------------------------------------------------------%
c
         istart = iend + 1
         if (iend .lt. kplusp) go to 20
c
c        %---------------------------------------------%
c        | Loop back to the top to get the next shift. |
c        %---------------------------------------------%
c
  110 continue
c
c     %--------------------------------------------------%
c     | Perform a similarity transformation that makes   |
c     | sure that H will have non negative sub diagonals |
c     %--------------------------------------------------%
c
      do 120 j=1,kev
         if ( h(j+1,j) .lt. zero ) then
              call dscal( kplusp-j+1, -one, h(j+1,j), ldh )
              call dscal( min(j+2, kplusp), -one, h(1,j+1), 1 )
              call dscal( min(j+np+1,kplusp), -one, q(1,j+1), 1 )
         end if
 120  continue
c
      do 130 i = 1, kev
c
c        %--------------------------------------------%
c        | Final check for splitting and deflation.   |
c        | Use a standard test as in the QR algorithm |
c        | REFERENCE: LAPACK subroutine dlahqr        |
c        %--------------------------------------------%
c
         tst1 = abs( h( i, i ) ) + abs( h( i+1, i+1 ) )
         if( tst1.eq.zero )
     &       tst1 = dlanhs( '1', kev, h, ldh, workl )
         if( h( i+1,i ) .le. max( ulp*tst1, smlnum ) )
     &       h(i+1,i) = zero
 130  continue
c
c     %-------------------------------------------------%
c     | Compute the (kev+1)-st column of (V*Q) and      |
c     | temporarily store the result in WORKD(N+1:2*N). |
c     | This is needed in the residual update since we  |
c     | cannot GUARANTEE that the corresponding entry   |
c     | of H would be zero as in exact arithmetic.      |
c     %-------------------------------------------------%
c
      if (h(kev+1,kev) .gt. zero)
     &    call dgemv ('N', n, kplusp, one, v, ldv, q(1,kev+1), 1, zero,
     &                workd(n+1), 1)
c
c     %----------------------------------------------------------%
c     | Compute column 1 to kev of (V*Q) in backward order       |
c     | taking advantage of the upper Hessenberg structure of Q. |
c     %----------------------------------------------------------%
c
      do 140 i = 1, kev
         call dgemv ('N', n, kplusp-i+1, one, v, ldv,
     &               q(1,kev-i+1), 1, zero, workd, 1)
         call dcopy (n, workd, 1, v(1,kplusp-i+1), 1)
  140 continue
c
c     %-------------------------------------------------%
c     |  Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). |
c     %-------------------------------------------------%
c
      call dlacpy ('A', n, kev, v(1,kplusp-kev+1), ldv, v, ldv)
c
c     %--------------------------------------------------------------%
c     | Copy the (kev+1)-st column of (V*Q) in the appropriate place |
c     %--------------------------------------------------------------%
c
      if (h(kev+1,kev) .gt. zero)
     &   call dcopy (n, workd(n+1), 1, v(1,kev+1), 1)
c
c     %-------------------------------------%
c     | Update the residual vector:         |
c     |    r &lt;- sigmak*r + betak*v(:,kev+1) |
c     | where                               |
c     |    sigmak = (e_{kplusp}'*Q)*e_{kev} |
c     |    betak = e_{kev+1}'*H*e_{kev}     |
c     %-------------------------------------%
c
      call dscal (n, q(kplusp,kev), resid, 1)
      if (h(kev+1,kev) .gt. zero)
     &   call daxpy (n, h(kev+1,kev), v(1,kev+1), 1, resid, 1)
c
c       if (msglvl .gt. 1) then
c          call dvout ( 1, q(kplusp,kev), ndigit,
c      &        '_napps: sigmak = (e_{kev+p}^T*Q)*e_{kev}')
c          call dvout ( 1, h(kev+1,kev), ndigit,
c      &        '_napps: betak = e_{kev+1}^T*H*e_{kev}')
c          call ivout ( 1, kev, ndigit,
c      &               '_napps: Order of the final Hessenberg matrix ')
c          if (msglvl .gt. 2) then
c             call dmout (logfil, kev, kev, h, ldh, ndigit,
c      &      '_napps: updated Hessenberg matrix H for next iteration')
c          end if
c c
c       end if
c
 9000 continue
*      call arscnd (t1)
c       tnapps = tnapps + (t1 - t0)
c
      return
c
c     %---------------%
c     | End of dnapps |
c     %---------------%
c
      end