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.