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 |

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 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 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.