Numérotation des lignes :

dnaitr
C DNAITR    SOURCE    FANDEUR   22/05/02    21:15:10     11359          c-----------------------------------------------------------------------c\BeginDoccc\Name: dnaitrcc\Description:c  Reverse communication interface for applying NP additional steps toc  a K step nonsymmetric Arnoldi factorization.cc  Input:  OP*V_{k}  -  V_{k}*H = r_{k}*e_{k}^Tcc          with (V_{k}^T)*B*V_{k} = I, (V_{k}^T)*B*r_{k} = 0.cc  Output: OP*V_{k+p}  -  V_{k+p}*H = r_{k+p}*e_{k+p}^Tcc          with (V_{k+p}^T)*B*V_{k+p} = I, (V_{k+p}^T)*B*r_{k+p} = 0.cc  where OP and B are as in dnaupd.  The B-norm of r_{k+p} is alsoc  computed and returned.cc\Usage:c  call dnaitrc     ( IDO, BMAT, N, K, NP, NB, RESID, RNORM, V, LDV, H, LDH,c       IPNTR, WORKD, INFO )cc\Argumentsc  IDO     Integer.  (INPUT/OUTPUT)c          Reverse communication flag.c          -------------------------------------------------------------c          IDO =  0: first call to the reverse communication interfacec          IDO = -1: compute  Y = OP * X  wherec                    IPNTR(1) is the pointer into WORK for X,c                    IPNTR(2) is the pointer into WORK for Y.c                    This is for the restart phase to force the newc                    starting vector into the range of OP.c          IDO =  1: compute  Y = OP * X  wherec                    IPNTR(1) is the pointer into WORK for X,c                    IPNTR(2) is the pointer into WORK for Y,c                    IPNTR(3) is the pointer into WORK for B * X.c          IDO =  2: compute  Y = B * X  wherec                    IPNTR(1) is the pointer into WORK for X,c                    IPNTR(2) is the pointer into WORK for Y.c          IDO = 99: donec          -------------------------------------------------------------c          When the routine is used in the "shift-and-invert" mode, thec          vector B * Q is already available and do not need to bec          recompute in forming OP * Q.cc  BMAT    Character*1.  (INPUT)c          BMAT specifies the type of the matrix B that defines thec          semi-inner product for the operator OP.  See dnaupd.c          B = 'I' -> standard eigenvalue problem A*x = lambda*xc          B = 'G' -> generalized eigenvalue problem A*x = lambda*M**xcc  N       Integer.  (INPUT)c          Dimension of the eigenproblem.cc  K       Integer.  (INPUT)c          Current size of V and H.cc  NP      Integer.  (INPUT)c          Number of additional Arnoldi steps to take.cc  NB      Integer.  (INPUT)c          Blocksize to be used in the recurrence.c          Only work for NB = 1 right now.  The goal is to have ac          program that implement both the block and non-block method.cc  RESID   Double precision array of length N.  (INPUT/OUTPUT)c          On INPUT:  RESID contains the residual vector r_{k}.c          On OUTPUT: RESID contains the residual vector r_{k+p}.cc  RNORM   Double precision scalar.  (INPUT/OUTPUT)c          B-norm of the starting residual on input.c          B-norm of the updated residual r_{k+p} on output.cc  V       REAL*8 N by K+NP array.  (INPUT/OUTPUT)c          On INPUT:  V contains the Arnoldi vectors in the first Kc          columns.c          On OUTPUT: V contains the new NP Arnoldi vectors in the nextc          NP columns.  The first K columns are unchanged.cc  LDV     Integer.  (INPUT)c          Leading dimension of V exactly as declared in the callingc          program.cc  H       REAL*8 (K+NP) by (K+NP) array.  (INPUT/OUTPUT)c          H is used to store the generated upper Hessenberg matrix.cc  LDH     Integer.  (INPUT)c          Leading dimension of H exactly as declared in the callingc          program.cc  IPNTR   Integer array of length 3.  (OUTPUT)c          Pointer to mark the starting locations in the WORK forc          vectors used by the Arnoldi iteration.c          -------------------------------------------------------------c          IPNTR(1): pointer to the current operand vector X.c          IPNTR(2): pointer to the current result vector Y.c          IPNTR(3): pointer to the vector B * X when used in thec                    shift-and-invert mode.  X is the current operand.c          -------------------------------------------------------------cc  WORKD   Double precision work array of length 3*N.  (REVERSE COMMUNICATION)c          Distributed array to be used in the basic Arnoldi iterationc          for reverse communication.  The calling program should notc          use WORKD as temporary workspace during the iteration !!!!!!c          On input, WORKD(1:N) = B*RESID and is used to save somec          computation at the first step.cc  INFO    Integer.  (OUTPUT)c          = 0: Normal exit.c          > 0: Size of the spanning invariant subspace of OP found.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     dgetv0  ARPACK routine to generate the initial vector.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 matricesc     dvout   ARPACK utility routine that prints vectors.c     dlabad  LAPACK routine that computes machine constants.c     dlamch  LAPACK routine that determines machine constants.c     dlascl  LAPACK routine for careful scaling of a matrix.c     dlanhs  LAPACK routine that computes various norms of a matrix.c     dgemv   Level 2 BLAS routine for matrix vector multiplication.c     daxpy   Level 1 BLAS that computes a vector triad.c     dscal   Level 1 BLAS that scales a vector.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\Revision history:c     xx/xx/92: Version ' 2.4'cc\SCCS Information: @(#)c FILE: naitr.F   SID: 2.4   DATE OF SID: 8/27/96   RELEASE: 2cc\Remarksc  The algorithm implemented is:cc  restart = .false.c  Given V_{k} = [v_{1}, ..., v_{k}], r_{k};c  r_{k} contains the initial residual vector even for k = 0;c  Also assume that rnorm = || B*r_{k} || and B*r_{k} are alreadyc  computed by the calling program.cc  betaj = rnorm ; p_{k+1} = B*r_{k} ;c  For  j = k+1, ..., k+np  Doc     1) if ( betaj &lt; tol ) stop or restart depending on j.c        ( At present tol is zero )c        if ( restart ) generate a new starting vector.c     2) v_{j} = r(j-1)/betaj;  V_{j} = [V_{j-1}, v_{j}];c        p_{j} = p_{j}/betajc     3) r_{j} = OP*v_{j} where OP is defined as in dnaupdc        For shift-invert mode p_{j} = B*v_{j} is already available.c        wnorm = || OP*v_{j} ||c     4) Compute the j-th step residual vector.c        w_{j} =  V_{j}^T * B * OP * v_{j}c        r_{j} =  OP*v_{j} - V_{j} * w_{j}c        H(:,j) = w_{j};c        H(j,j-1) = rnormc        rnorm = || r_(j) ||c        If (rnorm > 0.717*wnorm) accept step and go back to 1)c     5) Re-orthogonalization step:c        s = V_{j}'*B*r_{j}c        r_{j} = r_{j} - V_{j}*s;  rnorm1 = || r_{j} ||c        alphaj = alphaj + s_{j};c     6) Iterative refinement step:c        If (rnorm1 > 0.717*rnorm) thenc           rnorm = rnorm1c           accept step and go back to 1)c        Elsec           rnorm = rnorm1c           If this is the first time in step 6), go to 5)c           Else r_{j} lies in the span of V_{j} numerically.c              Set r_{j} = 0 and rnorm = 0; go to 1)c        EndIfc  End Docc\EndLibcc-----------------------------------------------------------------------c      subroutine dnaitr     &   (ido, bmat, n, k, np, nb, resid, rnorm, v, ldv, h, ldh,     &    ipntr, workd, ITRAK, info)cc     %----------------------------------------------------%c     | Include files for debugging and timing information |c -INC TARTRAKc     %----------------------------------------------------%ccc     %------------------%c     | Scalar Arguments |c     %------------------%c      character  bmat*1      integer    ido, info, k, ldh, ldv, n, nb, np      REAL*8     &           rnormc       REAL*8     T0,T1,T2,T3,T4,T5cc     %-----------------%c     | Array Arguments |c     %-----------------%c      integer    ipntr(3)      INTEGER    ITRAK(5)      REAL*8     &           h(ldh,k+np), resid(n), v(ldv,k+np), workd(3*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      logical    first, orth1, orth2, rstart, step3, step4      integer    ierr, i, infol, ipj, irj, ivj, iter, itry, j, msglvl,     &           jj      REAL*8     &           betaj, ovfl, temp1, rnorm1, smlnum, tst1, ulp, unfl,     &           wnorm      save       first, orth1, orth2, rstart, step3, step4,c      &           ierr, ipj, irj, ivj, iter, itry, j, msglvl, ovfl,     &           ierr, ipj, irj, ivj, iter, itry, j, ovfl,     &           betaj, rnorm1, smlnum, ulp, unfl, wnorm      parameter (msglvl=0)cc     %-----------------------%c     | Local Array Arguments |c     %-----------------------%c      REAL*8     &           xtemp(2)cc     %----------------------%c     | External Subroutines |c     %----------------------%c      external   daxpy, dcopy, dscal, dgemv,     &           dgetv0, dlabad,     &           dvout, dmout, ivoutcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           ddot, dnrm2, dlanhs, dlamch      external   ddot, dnrm2, dlanhs, dlamchcc     %---------------------%**c     | Intrinsic Functions |**c     %---------------------%**c**      intrinsic    abs, sqrt**c**c     %-----------------%**c     | Data statements |**c     %-----------------%**c      data      first / .true. /**c**c     %-----------------------%**c     | Executable Statements |c     %-----------------------%c       T0 = 0.D0c       T1 = 0.D0c       T2 = 0.D0c       T3 = 0.D0c       T4 = 0.D0c       T5 = 0.D0        NOPX  =ITRAK(1)        NBX   =ITRAK(2)        NRORTH=ITRAK(3)        NITREF=ITRAK(4)        NRSTRT=ITRAK(5)c      if (first) thencc        %-----------------------------------------%c        | Set machine-dependent constants for the |c        | the splitting and deflation criterion.  |c        | 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 ifc      if (ido .eq. 0) thencc        %-------------------------------%c        | Initialize timing statistics  |c        | & message level for debugging |c        %-------------------------------%c*         call arscnd (t0)c          msglvl = mnaitrcc        %------------------------------%c        | Initial call to this routine |c        %------------------------------%c         info   = 0         step3  = .false.         step4  = .false.         rstart = .false.         orth1  = .false.         orth2  = .false.         j      = k + 1         ipj    = 1         irj    = ipj   + n         ivj    = irj   + n      end ifcc     %-------------------------------------------------%c     | When in reverse communication mode one of:      |c     | STEP3, STEP4, ORTH1, ORTH2, RSTART              |c     | will be .true. when ....                        |c     | STEP3: return from computing OP*v_{j}.          |c     | STEP4: return from computing B-norm of OP*v_{j} |c     | ORTH1: return from computing B-norm of r_{j+1}  |c     | ORTH2: return from computing B-norm of          |c     |        correction to the residual vector.       |c     | RSTART: return from OP computations needed by   |c     |         dgetv0.                                 |c     %-------------------------------------------------%c      if (step3)  go to 50      if (step4)  go to 60      if (orth1)  go to 70      if (orth2)  go to 90      if (rstart) go to 30cc     %-----------------------------%c     | Else this is the first step |c     %-----------------------------%cc     %--------------------------------------------------------------%c     |                                                              |c     |        A R N O L D I     I T E R A T I O N     L O O P       |c     |                                                              |c     | Note:  B*r_{j-1} is already in WORKD(1:N)=WORKD(IPJ:IPJ+N-1) |c     %--------------------------------------------------------------%  1000 continuecc          if (msglvl .gt. 1) thenc             call ivout ( 1, j, ndigit,c      &                  '_naitr: generating Arnoldi vector number')c             call dvout ( 1, rnorm, ndigit,c      &                  '_naitr: B-norm of the current residual is')c          end ifcc        %---------------------------------------------------%c        | STEP 1: Check if the B norm of j-th residual      |c        | vector is zero. Equivalent to determing whether   |c        | an exact j-step Arnoldi factorization is present. |c        %---------------------------------------------------%c         betaj = rnorm         if (rnorm .gt. zero) go to 40cc           %---------------------------------------------------%c           | Invariant subspace found, generate a new starting |c           | vector which is orthogonal to the current Arnoldi |c           | basis and continue the iteration.                 |c           %---------------------------------------------------%c            if (msglvl .gt. 0) then               call ivout ( 1, j, ndigit,     &                     '_naitr: ****** RESTART AT STEP ******')            end ifcc           %---------------------------------------------%c           | ITRY is the loop variable that controls the |c           | maximum amount of times that a restart is   |c           %---------------------------------------------%c            betaj  = zero            nrstrt = nrstrt + 1            itry   = 1   20       continue            rstart = .true.            ido    = 0   30       continuecc           %--------------------------------------%c           | If in reverse communication mode and |c           | RSTART = .true. flow returns here.   |c           %--------------------------------------%c            call dgetv0 (ido, bmat, itry, .false., n, j, v, ldv,     &                   resid, rnorm, ipntr, workd, nopx,nbx, ierr)            if (ido .ne. 99) go to 9000            if (ierr .lt. 0) then               itry = itry + 1               if (itry .le. 3) go to 20cc              %------------------------------------------------%c              | Give up after several restart attempts.        |c              | Set INFO to the size of the invariant subspace |c              | which spans OP and exit.                       |c              %------------------------------------------------%c               info = j - 1*               call arscnd (t1)c                tnaitr = tnaitr + (t1 - t0)               ido = 99               go to 9000            end ifc   40    continuecc        %---------------------------------------------------------%c        | STEP 2:  v_{j} = r_{j-1}/rnorm and p_{j} = p_{j}/rnorm  |c        | Note that p_{j} = B*r_{j-1}. In order to avoid overflow |c        | when reciprocating a small RNORM, test against lower    |c        | machine bound.                                          |c        %---------------------------------------------------------%c         call dcopy (n, resid, 1, v(1,j), 1)         if (rnorm .ge. unfl) then             temp1 = one / rnorm             call dscal (n, temp1, v(1,j), 1)             call dscal (n, temp1, workd(ipj), 1)         elsecc            %-----------------------------------------%c            | To scale both v_{j} and p_{j} carefully |c            | use LAPACK routine SLASCL               |c            %-----------------------------------------%c             call dlascl ('General', i, i, rnorm, one, n, 1,     &                    v(1,j), n, infol)             call dlascl ('General', i, i, rnorm, one, n, 1,     &                    workd(ipj), n, infol)         end ifcc        %------------------------------------------------------%c        | STEP 3:  r_{j} = OP*v_{j}; Note that p_{j} = B*v_{j} |c        | Note that this is not quite yet r_{j}. See STEP 4    |c        %------------------------------------------------------%c         step3 = .true.         nopx  = nopx + 1*         call arscnd (t2)         call dcopy (n, v(1,j), 1, workd(ivj), 1)         ipntr(1) = ivj         ipntr(2) = irj         ipntr(3) = ipj         ido = 1cc        %-----------------------------------%c        | Exit in order to compute OP*v_{j} |c        %-----------------------------------%c         go to 9000   50    continuecc        %----------------------------------%c        | Back from reverse communication; |c        | WORKD(IRJ:IRJ+N-1) := OP*v_{j}   |c        | if step3 = .true.                |c        %----------------------------------%c*         call arscnd (t3)c          tmvopx = tmvopx + (t3 - t2)          step3 = .false.cc        %------------------------------------------%c        | Put another copy of OP*v_{j} into RESID. |c        %------------------------------------------%c         call dcopy (n, workd(irj), 1, resid, 1)cc        %---------------------------------------%c        | STEP 4:  Finish extending the Arnoldi |c        |          factorization to length j.   |c        %---------------------------------------%c*         call arscnd (t2)         if (bmat .eq. 'G') then            nbx = nbx + 1            step4 = .true.            ipntr(1) = irj            ipntr(2) = ipj            ido = 2cc           %-------------------------------------%c           | Exit in order to compute B*OP*v_{j} |c           %-------------------------------------%c            go to 9000         else if (bmat .eq. 'I') then            call dcopy (n, resid, 1, workd(ipj), 1)         end if   60    continuecc        %----------------------------------%c        | Back from reverse communication; |c        | WORKD(IPJ:IPJ+N-1) := B*OP*v_{j} |c        | if step4 = .true.                |c        %----------------------------------%cc          if (bmat .eq. 'G') then*            call arscnd (t3)c             tmvbx = tmvbx + (t3 - t2)c          end ifc         step4 = .false.cc        %-------------------------------------%c        | The following is needed for STEP 5. |c        | Compute the B-norm of OP*v_{j}.     |c        %-------------------------------------%c         if (bmat .eq. 'G') then             wnorm = ddot (n, resid, 1, workd(ipj), 1)             wnorm = sqrt(abs(wnorm))         else if (bmat .eq. 'I') then            wnorm = dnrm2(n, resid, 1)         end ifcc        %-----------------------------------------%c        | Compute the j-th residual corresponding |c        | to the j step factorization.            |c        | Use Classical Gram Schmidt and compute: |c        | w_{j} &lt;-  V_{j}^T * B * OP * v_{j}      |c        | r_{j} &lt;-  OP*v_{j} - V_{j} * w_{j}      |c        %-----------------------------------------%ccc        %------------------------------------------%c        | Compute the j Fourier coefficients w_{j} |c        | WORKD(IPJ:IPJ+N-1) contains B*OP*v_{j}.  |c        %------------------------------------------%c         call dgemv ('T', n, j, one, v, ldv, workd(ipj), 1,     &               zero, h(1,j), 1)cc        %--------------------------------------%c        | Orthogonalize r_{j} against V_{j}.   |c        | RESID contains OP*v_{j}. See STEP 3. |c        %--------------------------------------%c         call dgemv ('N', n, j, -one, v, ldv, h(1,j), 1,     &               one, resid, 1)c         if (j .gt. 1) h(j,j-1) = betajc*         call arscnd (t4)c         orth1 = .true.c*         call arscnd (t2)         if (bmat .eq. 'G') then            nbx = nbx + 1            call dcopy (n, resid, 1, workd(irj), 1)            ipntr(1) = irj            ipntr(2) = ipj            ido = 2cc           %----------------------------------%c           | Exit in order to compute B*r_{j} |c           %----------------------------------%c            go to 9000         else if (bmat .eq. 'I') then            call dcopy (n, resid, 1, workd(ipj), 1)         end if   70    continuecc        %---------------------------------------------------%c        | Back from reverse communication if ORTH1 = .true. |c        | WORKD(IPJ:IPJ+N-1) := B*r_{j}.                    |c        %---------------------------------------------------%cc          if (bmat .eq. 'G') then*            call arscnd (t3)c             tmvbx = tmvbx + (t3 - t2)c          end ifc         orth1 = .false.cc        %------------------------------%c        | Compute the B-norm of r_{j}. |c        %------------------------------%c         if (bmat .eq. 'G') then            rnorm = ddot (n, resid, 1, workd(ipj), 1)            rnorm = sqrt(abs(rnorm))         else if (bmat .eq. 'I') then            rnorm = dnrm2(n, resid, 1)         end ifcc        %-----------------------------------------------------------%c        | STEP 5: Re-orthogonalization / Iterative refinement phase |c        | Maximum NITER_ITREF tries.                                |c        |                                                           |c        |          s      = V_{j}^T * B * r_{j}                     |c        |          r_{j}  = r_{j} - V_{j}*s                         |c        |          alphaj = alphaj + s_{j}                          |c        |                                                           |c        | The stopping criteria used for iterative refinement is    |c        | discussed in Parlett's book SEP, page 107 and in Gragg &  |c        | Reichel ACM TOMS paper; Algorithm 686, Dec. 1990.         |c        | Determine if we need to correct the residual. The goal is |c        | to enforce ||v(:,1:j)^T * r_{j}|| .le. eps * || r_{j} ||  |c        | The following test determines whether the sine of the     |c        | angle between  OP*x and the computed residual is less     |c        | than or equal to 0.717.                                   |c        %-----------------------------------------------------------%c         if (rnorm .gt. 0.717*wnorm) go to 100         iter  = 0         nrorth = nrorth + 1cc        %---------------------------------------------------%c        | Enter the Iterative refinement phase. If further  |c        | refinement is necessary, loop back here. The loop |c        | variable is ITER. Perform a step of Classical     |c        | Gram-Schmidt using all the Arnoldi vectors V_{j}  |c        %---------------------------------------------------%c   80    continuecc          if (msglvl .gt. 2) thenc             xtemp(1) = wnormc             xtemp(2) = rnormc             call dvout ( 2, xtemp, ndigit,c      &           '_naitr: re-orthonalization; wnorm and rnorm are')c             call dvout ( j, h(1,j), ndigit,c      &                  '_naitr: j-th column of H')c          end ifcc        %----------------------------------------------------%c        | Compute V_{j}^T * B * r_{j}.                       |c        | WORKD(IRJ:IRJ+J-1) = v(:,1:J)'*WORKD(IPJ:IPJ+N-1). |c        %----------------------------------------------------%c         call dgemv ('T', n, j, one, v, ldv, workd(ipj), 1,     &               zero, workd(irj), 1)cc        %---------------------------------------------%c        | Compute the correction to the residual:     |c        | r_{j} = r_{j} - V_{j} * WORKD(IRJ:IRJ+J-1). |c        | The correction to H is v(:,1:J)*H(1:J,1:J)  |c        | + v(:,1:J)*WORKD(IRJ:IRJ+J-1)*e'_j.         |c        %---------------------------------------------%c         call dgemv ('N', n, j, -one, v, ldv, workd(irj), 1,     &               one, resid, 1)         call daxpy (j, one, workd(irj), 1, h(1,j), 1)c         orth2 = .true.*         call arscnd (t2)         if (bmat .eq. 'G') then            nbx = nbx + 1            call dcopy (n, resid, 1, workd(irj), 1)            ipntr(1) = irj            ipntr(2) = ipj            ido = 2cc           %-----------------------------------%c           | Exit in order to compute B*r_{j}. |c           | r_{j} is the corrected residual.  |c           %-----------------------------------%c            go to 9000         else if (bmat .eq. 'I') then            call dcopy (n, resid, 1, workd(ipj), 1)         end if   90    continuecc        %---------------------------------------------------%c        | Back from reverse communication if ORTH2 = .true. |c        %---------------------------------------------------%cc          if (bmat .eq. 'G') then*            call arscnd (t3)c             tmvbx = tmvbx + (t3 - t2)c          end ifcc        %-----------------------------------------------------%c        | Compute the B-norm of the corrected residual r_{j}. |c        %-----------------------------------------------------%c         if (bmat .eq. 'G') then             rnorm1 = ddot (n, resid, 1, workd(ipj), 1)             rnorm1 = sqrt(abs(rnorm1))         else if (bmat .eq. 'I') then             rnorm1 = dnrm2(n, resid, 1)         end ifc         if (msglvl .gt. 0 .and. iter .gt. 0) then            call ivout ( 1, j, ndigit,     &           '_naitr: Iterative refinement for Arnoldi residual')c             if (msglvl .gt. 2) thenc                 xtemp(1) = rnormc                 xtemp(2) = rnorm1c                 call dvout ( 2, xtemp, ndigit,c      &           '_naitr: iterative refinement ; rnorm and rnorm1 are')c             end if         end ifcc        %-----------------------------------------%c        | Determine if we need to perform another |c        | step of re-orthogonalization.           |c        %-----------------------------------------%c         if (rnorm1 .gt. 0.717*rnorm) thencc           %---------------------------------------%c           | No need for further refinement.       |c           | The cosine of the angle between the   |c           | corrected residual vector and the old |c           | residual vector is greater than 0.717 |c           | In other words the corrected residual |c           | and the old residual vector share an  |c           | angle of less than arcCOS(0.717)      |c           %---------------------------------------%c            rnorm = rnorm1c         elsecc           %-------------------------------------------%c           | Another step of iterative refinement step |c           %-------------------------------------------%c            nitref = nitref + 1            rnorm  = rnorm1            iter   = iter + 1            if (iter .le. 1) go to 80cc           %-------------------------------------------------%c           | Otherwise RESID is numerically in the span of V |c           %-------------------------------------------------%c            do 95 jj = 1, n               resid(jj) = zero  95        continue            rnorm = zero         end ifcc        %----------------------------------------------%c        | Branch here directly if iterative refinement |c        | wasn't necessary or after at most NITER_REF  |c        | steps of iterative refinement.               |c        %----------------------------------------------%c  100    continuec         rstart = .false.         orth2  = .false.c*         call arscnd (t5)c          titref = titref + (t5 - t4)cc        %------------------------------------%c        | STEP 6: Update  j = j+1;  Continue |c        %------------------------------------%c         j = j + 1         if (j .gt. k+np) then*            call arscnd (t1)c             tnaitr = tnaitr + (t1 - t0)            ido = 99            do 110 i = max(1,k), k+np-1cc              %--------------------------------------------%c              | 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', k+np, h, ldh, workd(n+1) )               if( abs( h( i+1,i ) ).le.max( ulp*tst1, smlnum ) )     &              h(i+1,i) = zero 110        continuecc             if (msglvl .gt. 2) thenc                call dmout (logfil, k+np, k+np, h, ldh, ndigit,c      &          '_naitr: Final upper Hessenberg matrix H of order K+NP')c             end ifc            go to 9000         end ifcc        %--------------------------------------------------------%c        | Loop back to extend the factorization by another step. |c        %--------------------------------------------------------%c      go to 1000cc     %---------------------------------------------------------------%c     |                                                               |c     |  E N D     O F     M A I N     I T E R A T I O N     L O O P  |c     |                                                               |c     %---------------------------------------------------------------%c 9000 continue      ITRAK(1)=NOPX        ITRAK(2)=NBX         ITRAK(3)=NRORTH      ITRAK(4)=NITREF      ITRAK(5)=NRSTRTcc     %---------------%c     | End of dnaitr |c     %---------------%c      end      

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