Numérotation des lignes :

dgetv0
C DGETV0    SOURCE    FANDEUR   22/05/02    21:15:06     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dgetv0cc\Description:c  Generate a random initial residual vector for the Arnoldi process.c  Force the residual vector to be in the range of the operator OP.cc\Usage:c  call dgetv0c     ( IDO, BMAT, ITRY, INITV, N, J, V, LDV, RESID, RNORM,c       IPNTR, WORKD, IERR )cc\Argumentsc  IDO     Integer.  (INPUT/OUTPUT)c          Reverse communication flag.  IDO must be zero on the firstc          call to dgetv0.c          -------------------------------------------------------------c          IDO =  0: first call to the reverse communication interfacec          IDO = -1: compute  Y = OP * X  wherec                    IPNTR(1) is the pointer into WORKD for X,c                    IPNTR(2) is the pointer into WORKD for Y.c                    This is for the initialization phase to force thec                    starting vector into the range of OP.c          IDO =  2: compute  Y = B * X  wherec                    IPNTR(1) is the pointer into WORKD for X,c                    IPNTR(2) is the pointer into WORKD for Y.c          IDO = 99: donec          -------------------------------------------------------------cc  BMAT    Character*1.  (INPUT)c          BMAT specifies the type of the matrix B in the (generalized)c          eigenvalue problem A*x = lambda*B*x.c          B = 'I' -> standard eigenvalue problem A*x = lambda*xc          B = 'G' -> generalized eigenvalue problem A*x = lambda*B*xcc  ITRY    Integer.  (INPUT)c          ITRY counts the number of times that dgetv0 is called.c          It should be set to 1 on the initial call to dgetv0.cc  INITV   Logical variable.  (INPUT)c          .TRUE.  => the initial residual vector is given in RESID.c          .FALSE. => generate a random initial residual vector.cc  N       Integer.  (INPUT)c          Dimension of the problem.cc  J       Integer.  (INPUT)c          Index of the residual vector to be generated, with respect toc          the Arnoldi process.  J > 1 in case of a "restart".cc  V       REAL*8 N by J array.  (INPUT)c          The first J-1 columns of V contain the current Arnoldi basisc          if this is a "restart".cc  LDV     Integer.  (INPUT)c          Leading dimension of V exactly as declared in the callingc          program.cc  RESID   Double precision array of length N.  (INPUT/OUTPUT)c          Initial residual vector to be generated.  If RESID isc          provided, force RESID into the range of the operator OP.cc  RNORM   Double precision scalar.  (OUTPUT)c          B-norm of the generated residual.cc  IPNTR   Integer array of length 3.  (OUTPUT)cc  WORKD   Double precision work array of length 2*N.  (REVERSE COMMUNICATION).c          On exit, WORK(1:N) = B*RESID to be used in SSAITR.cc  IERR    Integer.  (OUTPUT)c          =  0: Normal exit.c          = -1: Cannot generate a nontrivial restarted residual vectorc                in the range of the operator OP.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     arscnd  ARPACK utility routine for timing. -> deleted by BP in 2020c     dvout   ARPACK utility routine for vector output.c     dlarnv  LAPACK routine for generating a random vector.c     dgemv   Level 2 BLAS routine for matrix vector multiplication.c     dcopy   Level 1 BLAS that copies one vector to another.c     ddot    Level 1 BLAS that computes the scalar product of two vectors.c     dnrm2   Level 1 BLAS that computes the norm of a vector.cc\Authorc     Danny Sorensen               Phuong Vuc     Richard Lehoucq              CRPC / Rice Universityc     Dept. of Computational &     Houston, Texasc     Applied Mathematicsc     Rice Universityc     Houston, Texascc\SCCS Information: @(#)c FILE: getv0.F   SID: 2.7   DATE OF SID: 04/07/99   RELEASE: 2cc\EndLibcc-----------------------------------------------------------------------c      subroutine dgetv0     &   ( ido, bmat, itry, initv, n, j, v, ldv, resid, rnorm,     &     ipntr, workd, nopx,nbx, ierr )cc     %----------------------------------------------------%c     | Include files for debugging and timing information |c -INC TARTRAKc     %----------------------------------------------------%ccc     %------------------%c     | Scalar Arguments |c     %------------------%cc       REAL*8 T0c       REAL*8 T1c       REAL*8 T2c       REAL*8 T3      character  bmat*1      logical    initv      integer    ido, ierr, itry, j, ldv, n      REAL*8     &           rnormcc     %-----------------%c     | Array Arguments |c     %-----------------%c      integer    ipntr(3)      REAL*8     &           resid(n), v(ldv,j), 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 & Arrays |c     %------------------------%c      logical    first, inits, orth      integer    idist, iseed(4), iter, msglvl, jj      REAL*8     &           rnorm0      save       first, iseed, inits, iter, msglvl, orth, rnorm0cc     %----------------------%c     | External Subroutines |c     %----------------------%c      external   dcopy, dgemvc      external   dlarnv, dvout, arsncdcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           ddot, dnrm2      external   ddot, dnrm2cc     %---------------------%**c     | Intrinsic Functions |**c     %---------------------%**c**      intrinsic    abs, sqrt**c**c     %-----------------%**c     | Data Statements |**c     %-----------------%**c      data       inits /.true./**c**c     %-----------------------%**c     | Executable Statements |c     %-----------------------%c       T0 = 0.D0c       T1 = 0.D0c       T2 = 0.D0c       T3 = 0.D0ccc     %-----------------------------------%c     | Initialize the seed of the LAPACK |c     | random number generator           |c     %-----------------------------------%c      if (inits) then          iseed(1) = 1          iseed(2) = 3          iseed(3) = 5          iseed(4) = 7          inits = .false.      end ifc      if (ido .eq.  0) thencc        %-------------------------------%c        | Initialize timing statistics  |c        | & message level for debugging |c        %-------------------------------%c*         call arscnd (t0)         msglvl = mgetv0c         ierr   = 0         iter   = 0         first  = .FALSE.         orth   = .FALSE.cc        %-----------------------------------------------------%c        | Possibly generate a random starting vector in RESID |c        | Use a LAPACK random number generator used by the    |c        | matrix generation routines.                         |c        |    idist = 1: uniform (0,1)  distribution;          |c        |    idist = 2: uniform (-1,1) distribution;          |c        |    idist = 3: normal  (0,1)  distribution;          |c        %-----------------------------------------------------%c         if (.not.initv) then            idist = 2            call dlarnv (idist, iseed, n, resid)         end ifcc        %----------------------------------------------------------%c        | Force the starting vector into the range of OP to handle |c        | the generalized problem when B is possibly (singular).   |c        %----------------------------------------------------------%c*         call arscnd (t2)         if (bmat .eq. 'G') then            nopx = nopx + 1            ipntr(1) = 1            ipntr(2) = n + 1            call dcopy (n, resid, 1, workd, 1)            ido = -1            go to 9000         end if      end ifcc     %-----------------------------------------%c     | Back from computing OP*(initial-vector) |c     %-----------------------------------------%c      if (first) go to 20cc     %-----------------------------------------------%c     | Back from computing B*(orthogonalized-vector) |c     %-----------------------------------------------%c      if (orth)  go to 40cc       if (bmat .eq. 'G') then*         call arscnd (t3)c          tmvopx = tmvopx + (t3 - t2)c       end ifcc     %------------------------------------------------------%c     | Starting vector is now in the range of OP; r = OP*r; |c     | Compute B-norm of starting vector.                   |c     %------------------------------------------------------%c*      call arscnd (t2)      first = .TRUE.      if (bmat .eq. 'G') then         nbx = nbx + 1         call dcopy (n, workd(n+1), 1, resid, 1)         ipntr(1) = n + 1         ipntr(2) = 1         ido = 2         go to 9000      else if (bmat .eq. 'I') then         call dcopy (n, resid, 1, workd, 1)      end ifc   20 continuecc       if (bmat .eq. 'G') then*         call arscnd (t3)c          tmvbx = tmvbx + (t3 - t2)c       end ifc      first = .FALSE.      if (bmat .eq. 'G') then          rnorm0 = ddot (n, resid, 1, workd, 1)          rnorm0 = sqrt(abs(rnorm0))      else if (bmat .eq. 'I') then           rnorm0 = dnrm2(n, resid, 1)      end if      rnorm  = rnorm0cc     %---------------------------------------------%c     | Exit if this is the very first Arnoldi step |c     %---------------------------------------------%c      if (j .eq. 1) go to 50cc     %----------------------------------------------------------------c     | Otherwise need to B-orthogonalize the starting vector against |c     | the current Arnoldi basis using Gram-Schmidt with iter. ref.  |c     | This is the case where an invariant subspace is encountered   |c     | in the middle of the Arnoldi factorization.                   |c     |                                                               |c     |       s = V^{T}*B*r;   r = r - V*s;                           |c     |                                                               |c     | Stopping criteria used for iter. ref. is discussed in         |c     | Parlett's book, page 107 and in Gragg & Reichel TOMS paper.   |c     %---------------------------------------------------------------%c      orth = .TRUE.   30 continue       call dgemv ('T', n, j-1, one, v, ldv, workd, 1,     &            zero, workd(n+1), 1)      call dgemv ('N', n, j-1, (-one), v, ldv, workd(n+1), 1,     &            one, resid, 1)cc     %----------------------------------------------------------%c     | Compute the B-norm of the orthogonalized starting vector |c     %----------------------------------------------------------%c*      call arscnd (t2)      if (bmat .eq. 'G') then         nbx = nbx + 1         call dcopy (n, resid, 1, workd(n+1), 1)         ipntr(1) = n + 1         ipntr(2) = 1         ido = 2         go to 9000      else if (bmat .eq. 'I') then         call dcopy (n, resid, 1, workd, 1)      end ifc   40 continuecc       if (bmat .eq. 'G') then*         call arscnd (t3)c          tmvbx = tmvbx + (t3 - t2)c       end ifc      if (bmat .eq. 'G') then         rnorm = ddot (n, resid, 1, workd, 1)         rnorm = sqrt(abs(rnorm))      else if (bmat .eq. 'I') then         rnorm = dnrm2(n, resid, 1)      end ifcc     %--------------------------------------%c     | Check for further orthogonalization. |c     %--------------------------------------%cc       if (msglvl .gt. 2) thenc           call dvout ( 1, rnorm0, ndigit,c      &                '_getv0: re-orthonalization ; rnorm0 is')c           call dvout ( 1, rnorm, ndigit,c      &                '_getv0: re-orthonalization ; rnorm is')c       end ifc      if (rnorm .gt. 0.717*rnorm0) go to 50c      iter = iter + 1      if (iter .le. 5) thencc        %-----------------------------------%c        | Perform iterative refinement step |c        %-----------------------------------%c         rnorm0 = rnorm         go to 30      elsecc        %------------------------------------%c        | Iterative refinement step "failed" |c        %------------------------------------%c         do 45 jj = 1, n            resid(jj) = zero   45    continue         rnorm = zero         ierr = -1      end ifc   50 continue c       if (msglvl .gt. 0) thenc          call dvout ( 1, rnorm, ndigit,c      &        '_getv0: B-norm of initial - restarted starting vector')c       end ifc       if (msglvl .gt. 3) thenc          call dvout ( n, resid, ndigit,c      &        '_getv0: initial - restarted starting vector')c       end if      ido = 99c*      call arscnd (t1)c       tgetv0 = tgetv0 + (t1 - t0)c 9000 continue      returncc     %---------------%c     | End of dgetv0 |c     %---------------%c      end      

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