Shift and invert on a Generalized Eigen-problem

  

Figure 3.1: Reverse communication interface for Shift-Invert.

     
c     %-----------------------------------------------------%
c     | Call a routine FAC to factor the matrix (A-sigma*M) |
c     |   into L*U.                                         |
c     |                                                     |
c     | A routine MV is called repeatedly below to          |
c     |   form z = Mv.                                      |
c     |                                                     |
c     | A routine SOLVE is used repeatedly below to solve   |
c     |   (A-sigma*M) w = z  using the single LU            |
c     |   factorization provided by FAC.                    |
c     %-----------------------------------------------------%
c
      call fac( (A-sigma*M), L, U)
c
c     %-------------------------------------------%
c     | M A I N   L O O P (Reverse communication) |
c     %-------------------------------------------%
c
 10   continue
c
c        %---------------------------------------------%
c        | Repeatedly call the routine DNAUPD and take |
c        | actions indicated by parameter IDO .        |
c        %---------------------------------------------%
c        
         call dnaupd ( ido, bmat, n, which, nev, tol, resid,
     &                       ncv, v, ldv, iparam, ipntr, workd,

         if (ido .eq. -1) then
c     
c           %--------------------------------------------%   
c           | Perform  y <-- OP*x = inv[A-sigma*M]*M*x   |   
c           | to force the starting vector into the      |   
c           | range of OP.                               |   
c           |               x = workd(ipntr(1))          |   
c           |               y = workd(ipntr(2))          |   
c           %--------------------------------------------%   
c
            call mv (workd(ipntr(1)), workd(ipntr(2)))
c

  

Figure 3.2: Reverse communication interface for Shift-Invert contd.

            call solve(L,U, workd(ipntr(2)))
c
c           %---  L O O P   B A C K to call DSAUPD again. ---%
c
            go to 10
c           
         else if (ido .eq. 1) then
c           
c           %-----------------------------------------%      
c           | Perform y <-- OP*x = inv[A-sigma*M]*M*x |      
c           | M*x has been saved in workd(ipntr(3)).  |      
c           |      M*x = workd(ipntr(3)               |   
c           |        y = workd(ipntr(2)).             |
c           %-----------------------------------------%
c
            call dcopy (n, workd(ipntr(3)), 1,
     &                                    workd(ipntr(2)), 1)
            call solve(L,U, workd(ipntr(2)))
c
c           %---  L O O P   B A C K to call DSAUPD again. ---%
c
            go to 10
c
        else if (ido .eq. 2) then
c
c           %-----------------------------------------%
c           |  Perform  y <--- M*x                    |
c           |               x = workd(ipntr(1))       |
c           |               y = workd(ipntr(2))       |
c           %-----------------------------------------%
c
            call mv (workd(ipntr(1)), workd(ipntr(2)))
c           
c           %---  L O O P   B A C K to call DSAUPD again. ---%
c
            go to 10
        end if

The above code segments shown in Figures 3.1-3.2  illustrate the reverse communication loop for dnaupd in shift-invert mode for a generalized non-symmetric eigenvalue problem. This loop structure will be identical for the symmetric problem. The only change needed is to replace dnaupd with dsaupd. The loop structure is also identical for the complex arithmetic subroutine znaupd.

In the example shown in Figures 3.1-3.2, the matrix ${\bf M}$ is assumed to be symmetric and positive semi-definite. In the structure above, the user will have to supply some routine such as fac to obtain a matrix factorization of that may repeatedly be used to solve linear systems. Moreover, a routine needs to be provided in place of mv to perform the matrix-vector product and a routine in place of solve is required to solve linear systems of the form as needed using the previously computed factorization.

When convergence has taken place (indicated by ido = 99), the reverse communication loop will be exited. Then, post-processing using the ARPACK subroutine dneupd must be done to recover the eigenvalues and corresponding eigenvectors of the original problem. When operating in Shift-invert mode, the eigenvalue selection parameter which is normally set to which = 'LM'. The routine dneupd is then used to convert the converged eigenvalues of OP to eigenvalues of the original problem (3.2.1). Also, when ${\bf M}$ is singular or ill-conditioned, the routine dneupd takes steps to purify the eigenvectors and rid them of numerical corruption from eigenvectors corresponding to near-infinite eigenvalues. These procedures are done automatically by the routine dneupd when operating in any one of the computational modes described above and later in this chapter.

The user may wish to construct alternative computational modes using spectral transformations that are not addressed by any of the modes specified in this chapter. The reverse communication interface will easily accommodate these modifications. However, it will most likely be necessary to construct explicit transformations of the eigenvalues of OP to eigenvalues of the original problem in these situations.