Numérotation des lignes :

C DSAITR    SOURCE    BP208322  20/02/06    21:15:26     10512          c-----------------------------------------------------------------------c\BeginDoccc\Name: dsaitrcc\Description:c  Reverse communication interface for applying NP additional steps toc  a K step symmetric 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 dsaupd.  The B-norm of r_{k+p} is alsoc  computed and returned.cc\Usage:c  call dsaitrc     ( IDO, BMAT, N, K, NP, MODE, 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 does not need to bec          recomputed in forming OP * Q.cc  BMAT    Character*1.  (INPUT)c          BMAT specifies the type of matrix B that defines thec          semi-inner product for the operator OP.  See dsaupd.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 order of H and the number of columns of V.cc  NP      Integer.  (INPUT)c          Number of additional Arnoldi steps to take.cc  MODE    Integer.  (INPUT)c          Signifies which form for "OP". If MODE=2 thenc          a reduction in the number of B matrix vector multipliesc          is possible since the B-norm of OP*x is equivalent toc          the inv(B)-norm of A*x.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          On INPUT the B-norm of r_{k}.c          On OUTPUT the B-norm of the updated residual r_{k+p}.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 2 array.  (INPUT/OUTPUT)c          H is used to store the generated symmetric tridiagonal matrixc          with the subdiagonal in the first column starting at H(2,1)c          and the main diagonal in the second column.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 where RESID is associatedc          with the K step Arnoldi factorization. Used to save somec          computation at the first step.c          On OUTPUT, WORKD(1:N) = B*RESID where RESID is associatedc          with the K+NP step Arnoldi factorization.cc  INFO    Integer.  (OUTPUT)c          = 0: Normal exit.c          > 0: Size of an invariant subspace of OP is found that isc               less than K + NP.cc\EndDoccc-----------------------------------------------------------------------cc\BeginLibcc\Local variables:c     xxxxxx  realcc\Routines called:c     dgetv0  ARPACK routine to generate the initial vector.c     ivout   ARPACK utility routine that prints integers.c     dmout   ARPACK utility routine that prints matrices.c     dvout   ARPACK utility routine that prints vectors.c     dlamch  LAPACK routine that determines machine constants.c     dlascl  LAPACK routine for careful scaling 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/93: Version ' 2.4'cc\SCCS Information: @(#)c FILE: saitr.F   SID: 2.6   DATE OF SID: 8/28/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        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 dsaupdc        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        alphaj &lt;- j-th component of w_{j}c        rnorm = || r_{j} ||c        betaj+1 = rnormc        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 dsaitr     &   (ido, bmat, n, k, np, mode, 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, mode, np      REAL*8     &           rnorm      real*8     T0,T1,T2,T3,T4,T5cc     %-----------------%c     | Array Arguments |c     %-----------------%c      integer    ipntr(3)      INTEGER    ITRAK(5)      REAL*8     &           h(ldh,2), 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    i, ierr, ipj, irj, ivj, iter, itry, j, msglvl,     &           infol, jj      REAL*8     &           rnorm1, wnorm, safmin, temp1      save       orth1, orth2, rstart, step3, step4,c      &           ierr, ipj, irj, ivj, iter, itry, j, msglvl,     &           ierr, ipj, irj, ivj, iter, itry, j,      &           rnorm1, safmin, 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, dvout, dmout,     &           dlascl, ivoutcc     %--------------------%c     | External Functions |c     %--------------------%c      REAL*8     &           ddot, dnrm2, dlamch      external   ddot, dnrm2, dlamchcc     %-----------------%c     | Data statements |c     %-----------------%c      data      first / .true. /cc     %-----------------------%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) then         first = .false.cc        %--------------------------------%c        | safmin = safe minimum is such  |c        | that 1/sfmin does not overflow |c        %--------------------------------%c         safmin = dlamch('safmin')      end ifc      if (ido .eq. 0) thencc        %-------------------------------%c        | Initialize timing statistics  |c        | & message level for debugging |c        %-------------------------------%c*         call arscnd (t0)c          msglvl = msaitrcc        %------------------------------%c        | Initial call to this routine |c        %------------------------------%c         info   = 0         step3  = .false.         step4  = .false.         rstart = .false.         orth1  = .false.         orth2  = .false.cc        %--------------------------------%c        | Pointer to the current step of |c        | the factorization to build     |c        %--------------------------------%c         j      = k + 1cc        %------------------------------------------%c        | Pointers used for reverse communication  |c        | when using WORKD.                        |c        %------------------------------------------%c         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.                                  |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     %--------------------------------------------------------------%c 1000 continuecc          if (msglvl .gt. 2) thenc             call ivout ( 1, j, ndigit,c      &                  '_saitr: generating Arnoldi vector no.')c             call dvout ( 1, rnorm, ndigit,c      &                  '_saitr: B-norm of the current residual =')c          end ifcc        %---------------------------------------------------------%c        | Check for exact zero. Equivalent to determing whether a |c        | j-step Arnoldi factorization is present.                |c        %---------------------------------------------------------%c         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,     &                     '_saitr: ****** 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            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, nbx,nopx, 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                tsaitr = tsaitr + (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. safmin) 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        %-----------------------------------%c*         call arscnd (t3)c          tmvopx = tmvopx + (t3 - t2)c         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 symmetric   |c        |          Arnoldi to length j. If MODE = 2 |c        |          then B*OP = B*inv(B)*A = A and   |c        |          we don't need to compute B*OP.   |c        | NOTE: If MODE = 2 WORKD(IVJ:IVJ+N-1) is   |c        | assumed to have A*v_{j}.                  |c        %-------------------------------------------%c         if (mode .eq. 2) go to 65*         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        %-----------------------------------%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   65    continue         if (mode .eq. 2) thencc           %----------------------------------%c           | Note that the B-norm of OP*v_{j} |c           | is the inv(B)-norm of A*v_{j}.   |c           %----------------------------------%c            wnorm = ddot (n, resid, 1, workd(ivj), 1)            wnorm = sqrt(abs(wnorm))         else 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         if (mode .ne. 2 ) then            call dgemv('T', n, j, one, v, ldv, workd(ipj), 1, zero,     &                  workd(irj), 1)         else if (mode .eq. 2) then            call dgemv('T', n, j, one, v, ldv, workd(ivj), 1, zero,     &                  workd(irj), 1)         end ifcc        %--------------------------------------%c        | Orthgonalize r_{j} against V_{j}.    |c        | RESID contains OP*v_{j}. See STEP 3. |c        %--------------------------------------%c         call dgemv('N', n, j, -one, v, ldv, workd(irj), 1, one,     &               resid, 1)cc        %--------------------------------------%c        | Extend H to have j rows and columns. |c        %--------------------------------------%c         h(j,2) = workd(irj + j - 1)         if (j .eq. 1  .or.  rstart) then            h(j,1) = zero         else            h(j,1) = rnorm         end if*         call arscnd (t4)c         orth1 = .true.         iter  = 0c*         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        %-----------------------------------------------------------%c         if (rnorm .gt. 0.717*wnorm) go to 100         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      &           '_saitr: re-orthonalization ; wnorm and rnorm are')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, but only   |c        | H(j,j) is updated.                           |c        %----------------------------------------------%c         call dgemv ('N', n, j, -one, v, ldv, workd(irj), 1,     &               one, resid, 1)c         if (j .eq. 1  .or.  rstart) h(j,1) = zero         h(j,2) = h(j,2) + workd(irj + 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,     &           '_saitr: Iterative refinement for Arnoldi residual')c             if (msglvl .gt. 2) thenc                 xtemp(1) = rnormc                 xtemp(2) = rnorm1c                 call dvout ( 2, xtemp, ndigit,c      &           '_saitr: 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           %--------------------------------%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        | Make sure the last off-diagonal element is non negative  |c        | If not perform a similarity transformation on H(1:j,1:j) |c        | and scale v(:,j) by -1.                                  |c        %----------------------------------------------------------%c         if (h(j,1) .lt. zero) then            h(j,1) = -h(j,1)            if ( j .lt. k+np) then               call dscal(n, -one, v(1,j+1), 1)            else               call dscal(n, -one, resid, 1)            end if         end ifcc        %------------------------------------%c        | STEP 6: Update  j = j+1;  Continue |c        %------------------------------------%c         j = j + 1         if (j .gt. k+np) then*            call arscnd (t1)c             tsaitr = tsaitr + (t1 - t0)            ido = 99cc             if (msglvl .gt. 1) thenc                call dvout ( k+np, h(1,2), ndigit,c      &         '_saitr: main diagonal of matrix H of step K+NP.')c                if (k+np .gt. 1) thenc                call dvout ( k+np-1, h(2,1), ndigit,c      &         '_saitr: sub diagonal of matrix H of step K+NP.')c                end ifc             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 dsaitr |c     %---------------%c      end

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