next up previous contents index
Next: Application to Reduced-Order Modeling Up: Band Lanczos Method   Previous: Basic Properties   Contents   Index

Algorithm

A complete statement of the band Lanczos algorithm is as follows.


\begin{algorithm}{Band Lanczos Method for NHEP
}
{
\begin{tabbing}
(nr)ss\=ijk...
...nvergence \\
\textup{(18)}
\> \> {\bf end for}
\end{tabbing}}
\end{algorithm}

We stress that the non-Hermitian band Lanczos method can be implemented by directly following the statement in Algorithm 7.16, together with the further details of some of the steps given below. To keep the length of the statement to one page, in Algorithm 7.16, all potentially nonzero entries of the banded parts of $T_j$ and $\tilde{T}_j$ are computed as inner products.

However, only roughly half of these entries need to be explicitly computed as inner products. The remaining entries can be obtained via the relation

\begin{displaymath}
t_{ik}=\tilde{t}_{ki}\, \delta_k / \delta_i,
\end{displaymath} (185)

which follows from the connection (7.71) of $T_j^{\rm (pr)}$ and $\tilde{T}_j^{\rm (pr)}$. In particular, in the following discussion of steps (10), (12), (13), and (14), we give formulas for the entries of $T_j$ and $\tilde{T}_j$ that use (7.79) whenever possible. Employing these formulas would minimize the number of inner products, but it would sacrifice some numerical stability.

Next, we discuss some of the steps of Algorithm 7.16 in more detail.

(4)
Generally, the decision if $\hat{v}_j$ needs to be deflated should be based on checking if
\begin{displaymath}
\Vert\hat{v}_j\Vert _2 \leq {\tt dtol},
\end{displaymath} (186)

where ${\tt dtol}$ is a suitably chosen small deflation tolerance. The vector $\hat{v}_j$ is then deflated if (7.80) is satisfied. In this case, the value $j-m_c$, if positive, is added to the index set ${\cal I}_w$, which contains the indices of the nonzero rows in the part $\tilde{T}_j^{\rm {(d)}}$ of $\tilde{T}_j$, and the current right block size is updated to $m_c=m_c-1$. If $m_c=0$, then the right block Krylov sequence (7.60) is exhausted, and the algorithm terminates naturally. If $m_c>0$, the vector $\hat{v_j}$ is deleted, the index $k+1$ of each of the remaining right candidate vectors $\hat{v}_{k+1}$ is reset to $k$, and, finally, the algorithm returns to step (3). If (7.80) is not satisfied, then no deflation is necessary and the algorithm proceeds with step (5).

(6)
In analogy to (7.80), the decision if $\hat{w}_j$ needs to be deflated is based on checking if
\begin{displaymath}
\Vert\hat{w}_j\Vert _2 \leq {\tt dtol}.
\end{displaymath} (187)

The vector $\hat{w}_j$ is then deflated if (7.81) is satisfied. In this case, the value $j-p_c$, if positive, is added to the index set ${\cal I}_v$, which contains the indices of the nonzero rows in the part $T_j^{\rm {(d)}}$ of $T_j$, and the current left block size is updated to $p_c=p_c-1$. If $p_c=0$, then the left block Krylov sequence (7.61) is exhausted, and the algorithm terminates naturally. If $p_c>0$, the vector $\hat{w}_j$ is deleted, the index $k+1$ of each of the remaining left candidate vectors $\hat{w}_{k+1}$ is reset to $k$, and finally, the algorithm returns to step (5). If (7.81) is not satisfied, then no deflation is necessary and the algorithm proceeds with step (7).

(7)
Both vectors $\hat{v}_j$ and $\hat{w}_j$ have passed the deflation check and are now normalized to become the next right and left Lanczos vectors $v_j$ and $w_j$, respectively. The normalization is such that

\begin{displaymath}
\left\Vert v_j \right\Vert _2 = \left\Vert w_j \right\Vert _2 = 1
\quad \mbox{for all}\quad j.
\end{displaymath}

(8)
Here, we compute $\delta_j$ and check for breakdown. If $\delta_j=0$, then look-ahead would be needed in order to continue the algorithm.

(9)
In this step, we advance the right block Krylov sequence by computing a new candidate vector, $\hat{v}_{j+m_c}$, as the $A$-multiple of the latest right Lanczos vector $v_j$.

(10)
The vector $\hat{v}_{j+m_c}$ is biorthogonalized against the left Lanczos vectors $w_k$, $k\in {\cal I}_v^{(\rm e)}$. Note that a TSMGS procedure is used to do this biorthogonalization.
One of the inner products required in the computation of $t_{kj}$ can be saved by employing (7.79). More precisely, the following formula for $t_{kj}$ should be used:
\begin{displaymath}
t_{kj}= \cases{
\displaystyle{\tilde{t}_{jk}\, \delta_j/\d...
...playstyle{{w_k^T \hat{v}_{j+m_c}}/{\delta_k}}
& \ otherwise.}
\end{displaymath} (188)

Note that the necessary entry $\tilde{t}_{jk}$ in (7.82) is available from step (7).

(11)
In this step, we advance the left block Krylov sequence by computing a new candidate vector, $\hat{w}_{j+p_c}$, as the $A^T$-multiple of the latest left Lanczos vector $w_j$.

(12)
The vector $\hat{w}_{j+p_c}$ is biorthogonalized against the right Lanczos vectors $v_k$, $k\in {\cal I}_w^{(\rm e)}$. Note that this biorthogonalization is done by means of a TSMGS procedure. Again, to save one inner product, the following formula for $\tilde{t}_{kj}$ should be used:

\begin{displaymath}
\tilde{t}_{kj} = \cases{
\displaystyle{t_{jk}\, \delta_j/\...
...playstyle{{\hat{w}_{j+p_c}^T v_k}/{\delta_k}}
& \ otherwise.}
\end{displaymath}

(13)
The right candidate vectors, $\hat{v}_{j+1},\hat{v}_{j+2},\ldots,\hat{v}_{j+m_c}$, are biorthogonalized against the latest left Lanczos vector $w_j$. To save inner products, the following formula for $t_{j,k-m_c}$ could be used:

\begin{displaymath}
t_{j,k-m_c} = \cases{
\displaystyle{{w_j^T \hat{v}_{k}}/{\d...
...lde{t}_{k-m_c,j}\, \delta_{k-m_c}/{\delta_j}}
& \ otherwise.}
\end{displaymath}

However, the use of this formula sacrifices some numerical stability.[*]

(14)
The left candidate vectors, $\hat{w}_{j+1},\hat{w}_{j+2},\ldots,\hat{w}_{j+p_c}$, are biorthogonalized against the latest right Lanczos vector $v_j$. To save inner products, the following formula for $\tilde{t}_{j,k-p_c}$ could be used:

\begin{displaymath}
\tilde{t}_{j,k-p_c} = \cases{
\displaystyle{{\hat{w}_{k}^T ...
...tyle{t_{k-p_c,j}\, \delta_{k-p_c}/{\delta_j}}
& \ otherwise.}
\end{displaymath}

However, the use of this formula sacrifices some numerical stability.

(15)
The entries $t_{jk}$ computed in this step are the potentially nonzero entries in the $j$th row of $T_j^{\rm {(pr)}}$ due to the vertical spikes caused by the deflation of $\hat{v}_k$ vectors. Note that again formula (7.79) is employed.

(16)
All the potentially nonzero elements in the $j$th row and the $j$th column of $T_j^{\rm {(pr)}}$ have now been computed and they are added to the matrix $T_{j-1}^{(\rm pr)}$ of the previous iteration $j-1$, to yield the current matrix $T_{j}^{(\rm pr)}$. Here, we use the convention that entries $t_{ik}$ that are not explicitly defined in Algorithm 7.16 are set to be zero. We remark that in the initial iterations, i.e., as long as $j\leq m_c$, respectively $j\leq p_c$, Algorithm 7.16 also produces entries $t_{jk}$, respectively $\tilde{t}_{jk}$, with negative indices $k\leq 0$. These entries arise from the biorthogonalization of the starting vectors, and they are not part of the matrix $T_j^{\rm {(pr)}}$. In particular, these entries are not needed when Algorithm 7.16 is used for eigenvalue computations. However, they are crucial when Algorithm 7.16 is applied to reduced-order modeling, as we will discuss in §7.10.4 below.

(17)
For eigenvalue computations, one tests for convergence by computing the eigenvalues $\theta_i^{(j)}$, $i=1,2,\ldots,j$, of the $j\times j$ matrix $T_j^{\rm {(pr)}}$. The algorithm is stopped if some of the $\theta_i^{(j)}$'s are good enough approximations to the desired eigenvalues of $A$. For reduced-order modeling, the algorithm is stopped if the $j$th-order model generated by the algorithm is a good enough approximation of the original linear dynamical system.


next up previous contents index
Next: Application to Reduced-Order Modeling Up: Band Lanczos Method   Previous: Basic Properties   Contents   Index
Susan Blackford 2000-11-20