Numérotation des lignes :

dnapps
C DNAPPS    SOURCE    FANDEUR   22/05/02    21:15:10     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dnappscc\Description:c  Given the Arnoldi factorizationcc     A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,cc  apply NP implicit shifts resulting incc     A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Qcc  where Q is an orthogonal matrix which is the product of rotationsc  and reflections resulting from the NP bulge chage sweeps.c  The updated Arnoldi factorization becomes:cc     A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.cc\Usage:c  call dnappsc     ( N, KEV, NP, SHIFTR, SHIFTI, V, LDV, H, LDH, RESID, Q, LDQ,c       WORKL, WORKD )cc\Argumentsc  N       Integer.  (INPUT)c          Problem size, i.e. size of matrix A.cc  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 onlyc          updated on ouput when fewer than NP shifts are applied inc          order to keep the conjugate pair together.cc  NP      Integer.  (INPUT)c          Number of implicit shifts to be applied.cc  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 thec          conjugate pairs are in consecutive locations.cc  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 vectorsc          in the first KEV columns of V.cc  LDV     Integer.  (INPUT)c          Leading dimension of V exactly as declared in the callingc          program.cc  H       REAL*8 (KEV+NP) by (KEV+NP) array.  (INPUT/OUTPUT)c          On INPUT, H contains the current KEV+NP by KEV+NP upperc          Hessenber matrix of the Arnoldi factorization.c          On OUTPUT, H contains the updated KEV by KEV upper Hessenbergc          matrix in the KEV leading submatrix.cc  LDH     Integer.  (INPUT)c          Leading dimension of H exactly as declared in the callingc          program.cc  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.cc  Q       REAL*8 KEV+NP by KEV+NP work array.  (WORKSPACE)c          Work array used to accumulate the rotations and reflectionsc          during the bulge chase sweep.cc  LDQ     Integer.  (INPUT)c          Leading dimension of Q exactly as declared in the callingc          program.cc  WORKL   Double precision work array of length (KEV+NP).  (WORKSPACE)c          Private (replicated) array on each PE or array allocated onc          the front end.cc  WORKD   Double precision work array of length 2*N.  (WORKSPACE)c          Distributed array used in the application of the accumulatedc          orthogonal matrix Q.cc\EndDoccc-----------------------------------------------------------------------cc\BeginLibcc\Local variables:c     xxxxxx  realcc\References:c  1. D.C. Sorensen, "Implicit Application of Polynomial Filters inc     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),c     pp 357-385.cc\Routines called:c     ivout   ARPACK utility routine that prints integers.c     arscnd  ARPACK utility routine for timing. -> deleted by BP in 2020c     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 toc             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.cc\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.4'cc\SCCS Information: @(#)c FILE: napps.F   SID: 2.4   DATE OF SID: 3/28/97   RELEASE: 2cc\Remarksc  1. In this version, each shift is applied to all the sublocks ofc     the Hessenberg matrix H and not just to the submatrix that itc     comes from. Deflation as in LAPACK routine dlahqr (QR algorithmc     for upper Hessenberg matrices ) is used.c     The subdiagonals of H are enforced to be non-negative.cc\EndLibcc-----------------------------------------------------------------------c      subroutine dnapps     &   ( n, kev, np, shiftr, shifti, v, ldv, h, ldh, resid, q, ldq,     &     workl, workd )cc     %----------------------------------------------------%c     | Include files for debugging and timing information |c -INC TARTRAKc     %----------------------------------------------------%ccc     %------------------%c     | Scalar Arguments |c     %------------------%c      integer    kev, ldh, ldq, ldv, n, npcc     %-----------------%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)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      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)cc     %----------------------%c     | External Subroutines |c     %----------------------%c      external   daxpy, dcopy, dscal,     &           dlacpy, dlarfg, dlarf,     &           dlaset, dlabad, dlartgcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           dlamch, dlanhs, dlapy2      external   dlamch, dlanhs, dlapy2cc     %----------------------%**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) thencc        %-----------------------------------------------%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 ifcc     %-------------------------------%c     | Initialize timing statistics  |c     | & message level for debugging |c     %-------------------------------%c*      call arscnd (t0)c       msglvl = mnapps      kplusp = kev + npcc     %--------------------------------------------%c     | Initialize Q to the identity to accumulate |c     | the rotations and reflections              |c     %--------------------------------------------%c      call dlaset ('All', kplusp, kplusp, zero, one, q, ldq)cc     %----------------------------------------------%c     | Quick return if there are no shifts to apply |c     %----------------------------------------------%c      if (np .eq. 0) go to 9000cc     %----------------------------------------------%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)cc          if (msglvl .gt. 2 ) thenc             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 ifcc        %-------------------------------------------------%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 ) thencc           %-----------------------------------------%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 ) thencc           %------------------------------------%c           | Start of a complex conjugate pair. |c           %------------------------------------%c            cconj = .true.         else if ( jj .eq. np .and. abs( sigmai ) .gt. zero ) thencc           %----------------------------------------------%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    continuecc        %--------------------------------------------------%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-1cc           %----------------------------------------%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    continuecc          if (msglvl .gt. 2) thenc              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 ifcc        %------------------------------------------------%c        | No reason to apply a shift to block of order 1 |c        %------------------------------------------------%c         if ( istart .eq. iend ) go to 100cc        %------------------------------------------------------%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 100c         h11 = h(istart,istart)         h21 = h(istart+1,istart)         if ( abs( sigmai ) .le. zero ) thencc           %---------------------------------------------%c           | Real-valued shift ==> apply single shift QR |c           %---------------------------------------------%c            f = h11 - sigmar            g = h21c            do 80 i = istart, iend-1cc              %-----------------------------------------------------%c              | Contruct the plane rotation G to zero out the bulge |c              %-----------------------------------------------------%c               call dlartg (f, g, c, s, r)               if (i .gt. istart) thencc                 %-------------------------------------------%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 ifcc              %---------------------------------------------%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          continuecc              %---------------------------------------------%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          continuecc              %----------------------------------------------------%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          continuecc              %---------------------------%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       continuecc           %-----------------------------------%c           | Finished applying the real shift. |c           %-----------------------------------%c         elsecc           %----------------------------------------------------%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)cc           %---------------------------------------------------------%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) = h32c            do 90 i = istart, iend-1c               nr = min ( 3, iend-i+1 )cc              %-----------------------------------------------------%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) = onecc              %--------------------------------------%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)cc              %---------------------------------------%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)cc              %-----------------------------------------------------%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)cc              %----------------------------%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 ifc   90       continuecc           %--------------------------------------------%c           | Finished applying a complex pair of shifts |c           | to the current block                       |c           %--------------------------------------------%c         end ifc  100    continuecc        %---------------------------------------------------------%c        | Apply the same shift to the next block if there is any. |c        %---------------------------------------------------------%c         istart = iend + 1         if (iend .lt. kplusp) go to 20cc        %---------------------------------------------%c        | Loop back to the top to get the next shift. |c        %---------------------------------------------%c  110 continuecc     %--------------------------------------------------%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  continuec      do 130 i = 1, kevcc        %--------------------------------------------%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  continuecc     %-------------------------------------------------%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)cc     %----------------------------------------------------------%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 continuecc     %-------------------------------------------------%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)cc     %--------------------------------------------------------------%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)cc     %-------------------------------------%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)cc       if (msglvl .gt. 1) thenc          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) thenc             call dmout (logfil, kev, kev, h, ldh, ndigit,c      &      '_napps: updated Hessenberg matrix H for next iteration')c          end ifc cc       end ifc 9000 continue*      call arscnd (t1)c       tnapps = tnapps + (t1 - t0)c      returncc     %---------------%c     | End of dnapps |c     %---------------%c      end     

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