Numérotation des lignes :

dsapps
C DSAPPS    SOURCE    FANDEUR   22/05/02    21:15:14     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dsappscc\Description:c  Given the Arnoldi factorizationcc     A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,cc  apply NP shifts implicitly 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 of order KEV+NP. Q is the product ofc  rotations resulting from the NP bulge chasing sweeps.  The updated Arnoldic  factorization becomes:cc     A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.cc\Usage:c  call dsappsc     ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ, WORKD )cc\Argumentsc  N       Integer.  (INPUT)c          Problem size, i.e. dimension of matrix A.cc  KEV     Integer.  (INPUT)c          INPUT: KEV+NP is the size of the input matrix H.c          OUTPUT: KEV is the size of the updated matrix HNEW.cc  NP      Integer.  (INPUT)c          Number of implicit shifts to be applied.cc  SHIFT   Double precision array of length NP.  (INPUT)c          The shifts to be applied.cc  V       REAL*8 N by (KEV+NP) array.  (INPUT/OUTPUT)c          INPUT: V contains the current KEV+NP Arnoldi vectors.c          OUTPUT: VNEW = V(1:n,1:KEV); the updated Arnoldi vectorsc          are 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 2 array.  (INPUT/OUTPUT)c          INPUT: H contains the symmetric tridiagonal matrix of thec          Arnoldi factorization with the subdiagonal in the 1st columnc          starting at H(2,1) and the main diagonal in the 2nd column.c          OUTPUT: H contains the updated tridiagonal matrix in thec          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          INPUT: RESID contains the the residual vector r_{k+p}.c          OUTPUT: RESID is the updated residual vector rnew_{k}.cc  Q       REAL*8 KEV+NP by KEV+NP work array.  (WORKSPACE)c          Work array used to accumulate the rotations during the bulgec          chase sweep.cc  LDQ     Integer.  (INPUT)c          Leading dimension of Q exactly as declared in the callingc          program.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.c  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitlyc     Restarted Arnoldi Iteration", Rice University Technical Reportc     TR95-13, Department of Computational and Applied Mathematics.cc\Routines called:c     ivout   ARPACK utility routine that prints integers.c     arscnd  ARPACK utility routine for timing. -> deleted by BP in 2020c     dvout   ARPACK utility routine that prints vectors.c     dlamch  LAPACK routine that determines machine constants.c     dlartg  LAPACK Givens rotation construction routine.c     dlacpy  LAPACK matrix copy 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     12/16/93: Version ' 2.4'cc\SCCS Information: @(#)c FILE: sapps.F   SID: 2.6   DATE OF SID: 3/28/97   RELEASE: 2cc\Remarksc  1. In this version, each shift is applied to all the subblocks ofc     the tridiagonal matrix H and not just to the submatrix that itc     comes from. This routine assumes that the subdiagonal elementsc     of H that are stored in h(1:kev+np,1) are nonegative upon inputc     and enforce this condition upon output. This version incorporatesc     deflation. See code for documentation.cc\EndLibcc-----------------------------------------------------------------------c      subroutine dsapps     &   ( n, kev, np, shift, v, ldv, h, ldh, resid, q, ldq, 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,2), q(ldq,kev+np), resid(n), shift(np),     &           v(ldv,kev+np), workd(2*n)cc     %------------%c     | Parameters |c     %------------%c      REAL*8     &           one, zero      parameter (one = 1.0D+0, zero = 0.0D+0)cc     %---------------%c     | Local Scalars |c     %---------------%c      integer    i, iend, istart, itop, j, jj, kplusp, msglvl      parameter (msglvl=0)      logical    first      REAL*8     &           a1, a2, a3, a4, big, c, epsmch, f, g, r, s      save       epsmch, firstccc     %----------------------%c     | External Subroutines |c     %----------------------%c      external   daxpy, dcopy, dscal, dlacpy,     &           dlartg, dlaset, dvout,     &           ivout, dgemvcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           dlamch      external   dlamchcc     %----------------------%**c     | Intrinsics Functions |**c     %----------------------%**c**      intrinsic  abs**c**c     %----------------%**c     | Data statments |**c     %----------------%**c      data       first / .true. /**c**c     %-----------------------%**c     | Executable Statements |c     %-----------------------%c      if (first) then         epsmch = dlamch('Epsilon-Machine')         first = .false.      end if      itop = 1cc     %-------------------------------%c     | Initialize timing statistics  |c     | & message level for debugging |c     %-------------------------------%c*      call arscnd (t0)c       msglvl = msappsc      kplusp = kev + npcc     %----------------------------------------------%c     | Initialize Q to the identity matrix of order |c     | kplusp used to accumulate the rotations.     |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     | Apply the np shifts implicitly. Apply each shift to the  |c     | whole matrix and not just to the submatrix from which it |c     | comes.                                                   |c     %----------------------------------------------------------%c      do 90 jj = 1, npc         istart = itopcc        %----------------------------------------------------------%c        | Check for splitting and deflation. Currently we consider |c        | an off-diagonal element h(i+1,1) negligible if           |c        |         h(i+1,1) .le. epsmch*( |h(i,2)| + |h(i+1,2)| )   |c        | for i=1:KEV+NP-1.                                        |c        | If above condition tests true then we set h(i+1,1) = 0.  |c        | Note that h(1:KEV+NP,1) are assumed to be non negative.  |c        %----------------------------------------------------------%c   20    continuecc        %------------------------------------------------%c        | The following loop exits early if we encounter |c        | a negligible off diagonal element.             |c        %------------------------------------------------%c         do 30 i = istart, kplusp-1            big   = abs(h(i,2)) + abs(h(i+1,2))            if (h(i+1,1) .le. epsmch*big) then               if (msglvl .gt. 0) then                  call ivout ( 1, i, ndigit,     &                 '_sapps: deflation at row/column no.')                  call ivout ( 1, jj, ndigit,     &                 '_sapps: occured before shift number.')                  call dvout ( 1, h(i+1,1), ndigit,     &                 '_sapps: the corresponding off diagonal element')               end if               h(i+1,1) = zero               iend = i               go to 40            end if   30    continue         iend = kplusp   40    continuec         if (istart .lt. iend) thencc           %--------------------------------------------------------%c           | Construct the plane rotation G'(istart,istart+1,theta) |c           | that attempts to drive h(istart+1,1) to zero.          |c           %--------------------------------------------------------%c             f = h(istart,2) - shift(jj)             g = h(istart+1,1)             call dlartg (f, g, c, s, r)cc            %-------------------------------------------------------%c            | Apply rotation to the left and right of H;            |c            | H &lt;- G' * H * G,  where G = G(istart,istart+1,theta). |c            | This will create a "bulge".                           |c            %-------------------------------------------------------%c             a1 = c*h(istart,2)   + s*h(istart+1,1)             a2 = c*h(istart+1,1) + s*h(istart+1,2)             a4 = c*h(istart+1,2) - s*h(istart+1,1)             a3 = c*h(istart+1,1) - s*h(istart,2)             h(istart,2)   = c*a1 + s*a2             h(istart+1,2) = c*a4 - s*a3             h(istart+1,1) = c*a3 + s*a4cc            %----------------------------------------------------%c            | Accumulate the rotation in the matrix Q;  Q &lt;- Q*G |c            %----------------------------------------------------%c             do 60 j = 1, min(istart+jj,kplusp)                a1            =   c*q(j,istart) + s*q(j,istart+1)                q(j,istart+1) = - s*q(j,istart) + c*q(j,istart+1)                q(j,istart)   = a1   60        continueccc            %----------------------------------------------%c            | The following loop chases the bulge created. |c            | Note that the previous rotation may also be  |c            | done within the following loop. But it is    |c            | kept separate to make the distinction among  |c            | the bulge chasing sweeps and the first plane |c            | rotation designed to drive h(istart+1,1) to  |c            | zero.                                        |c            %----------------------------------------------%c             do 70 i = istart+1, iend-1cc               %----------------------------------------------%c               | Construct the plane rotation G'(i,i+1,theta) |c               | that zeros the i-th bulge that was created   |c               | by G(i-1,i,theta). g represents the bulge.   |c               %----------------------------------------------%c                f = h(i,1)                g = s*h(i+1,1)cc               %----------------------------------%c               | Final update with G(i-1,i,theta) |c               %----------------------------------%c                h(i+1,1) = c*h(i+1,1)                call dlartg (f, g, c, s, r)cc               %-------------------------------------------%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 ifcc               %--------------------------------------------%c               | Apply rotation to the left and right of H; |c               | H &lt;- G * H * G',  where G = G(i,i+1,theta) |c               %--------------------------------------------%c                h(i,1) = rc                a1 = c*h(i,2)   + s*h(i+1,1)                a2 = c*h(i+1,1) + s*h(i+1,2)                a3 = c*h(i+1,1) - s*h(i,2)                a4 = c*h(i+1,2) - s*h(i+1,1)c                h(i,2)   = c*a1 + s*a2                h(i+1,2) = c*a4 - s*a3                h(i+1,1) = c*a3 + s*a4cc               %----------------------------------------------------%c               | Accumulate the rotation in the matrix Q;  Q &lt;- Q*G |c               %----------------------------------------------------%c                do 50 j = 1, min( i+jj, kplusp )                   a1       =   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)   = a1   50           continuec   70        continuec         end ifcc        %--------------------------%c        | Update the block pointer |c        %--------------------------%c         istart = iend + 1cc        %------------------------------------------%c        | Make sure that h(iend,1) is non-negative |c        | If not then set h(iend,1) &lt;-- -h(iend,1) |c        | and negate the last column of Q.         |c        | We have effectively carried out a        |c        | similarity on transformation H           |c        %------------------------------------------%c         if (h(iend,1) .lt. zero) then             h(iend,1) = -h(iend,1)             call dscal(kplusp, -one, q(1,iend), 1)         end ifcc        %--------------------------------------------------------%c        | Apply the same shift to the next block if there is any |c        %--------------------------------------------------------%c         if (iend .lt. kplusp) go to 20cc        %-----------------------------------------------------%c        | Check if we can increase the the start of the block |c        %-----------------------------------------------------%c         do 80 i = itop, kplusp-1            if (h(i+1,1) .gt. zero) go to 90            itop  = itop + 1   80    continuecc        %-----------------------------------%c        | Finished applying the jj-th shift |c        %-----------------------------------%c   90 continuecc     %------------------------------------------%c     | All shifts have been applied. Check for  |c     | more possible deflation that might occur |c     | after the last shift is applied.         |c     %------------------------------------------%c      do 100 i = itop, kplusp-1         big   = abs(h(i,2)) + abs(h(i+1,2))         if (h(i+1,1) .le. epsmch*big) then            if (msglvl .gt. 0) then               call ivout ( 1, i, ndigit,     &              '_sapps: deflation at row/column no.')               call dvout ( 1, h(i+1,1), ndigit,     &              '_sapps: the corresponding off diagonal element')            end if            h(i+1,1) = zero         end if 100  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 not necessary if h(kev+1,1) = 0.         |c     %-------------------------------------------------%c      if ( h(kev+1,1) .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 that Q is an upper triangular matrix |c     | with lower bandwidth np.                              |c     | Place results in v(:,kplusp-kev:kplusp) temporarily.  |c     %-------------------------------------------------------%c      do 130 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)  130 continuecc     %-------------------------------------------------%c     |  Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). |c     %-------------------------------------------------%c      call dlacpy ('All', n, kev, v(1,np+1), ldv, v, ldv)cc     %--------------------------------------------%c     | Copy the (kev+1)-st column of (V*Q) in the |c     | appropriate place if h(kev+1,1) .ne. zero. |c     %--------------------------------------------%c      if ( h(kev+1,1) .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_{kev+p}'*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,1) .gt. zero)     &   call daxpy (n, h(kev+1,1), v(1,kev+1), 1, resid, 1)cc       if (msglvl .gt. 1) thenc          call dvout ( 1, q(kplusp,kev), ndigit,c      &      '_sapps: sigmak of the updated residual vector')c          call dvout ( 1, h(kev+1,1), ndigit,c      &      '_sapps: betak of the updated residual vector')c          call dvout ( kev, h(1,2), ndigit,c      &      '_sapps: updated main diagonal of H for next iteration')c          if (kev .gt. 1) thenc          call dvout ( kev-1, h(2,1), ndigit,c      &      '_sapps: updated sub diagonal of H for next iteration')c          end ifc       end ifc*      call arscnd (t1)c       tsapps = tsapps + (t1 - t0)c 9000 continuec       returncc     %---------------%c     | End of dsapps |c     %---------------%c      end     

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