Since the nodes of the dynamic fish model are not fully connected by
viscoelastic units, the effective stiffness matrix and,
consequently (since is diagonal), the system matrix , are sparse matrices. Sparsity can be exploited to increase the
computational efficiency of the factorization algorithms
[Press et al.1986]. We have employed the *Skyline Storage Scheme* for
storing , which is popular in numerical methods for finite
element analysis [Bathe and Wilson1976]. This storage scheme is most efficient
for storing sparse matrices that are symmetric and banded. We will now
explain our strategy in more detail.

The system matrix and effective load matrix are computed according to Eq. () and Eq. (), respectively. To assemble , we first assemble the instantaneous, internal effective viscoelastic stiffnesses of the system (see Eq. ()). Denoting the entries of as and the rows of as , the following procedure assembles and :

- Initialize and .
- For each viscoelastic unit , assemble the term of Eq. () into :
and assemble the term of Eq. () into :

- Assemble the remaining terms of and . For
:

Fig. shows the structure of the matrix
which depends solely on the connectivity of the *n*=23
nodes in the biomechanical fish model. Note that the matrix is
symmetric, sparse and banded. Bandedness means that all off-diagonal
entries beyond the bandwidth of the matrix are zero. Because
is symmetric, we can state the bandedness condition as

where is the *bandwidth* of and
is called the *half-bandwidth* of . For
example, if , then is a diagonal
matrix. For our dynamic fish model, (see
Fig. ).

The skyline storage scheme exploits the fact that the many zeros outside the skyline of are unnecessary in the factorization solution of Eq. . Hence, an effective way of storing the symmetric would be to store only the upper triangular half of the matrix below the skyline, including the diagonal. The skyline storage scheme stores into a one-dimensional array . In addition, it uses an auxiliary array to properly address the entries of stored in .

We explain the matrix storage and addressing scheme via a small
example system matrix of order 4 shown in
Fig. . Let be the row number of the
first nonzero entry in column *i* of . The variables ,
*i*=0, 1, 2, ..., define the *skyline* of the matrix and
are the *column heights*. Furthermore, the half-bandwidth of the
matrix, , equals . With the column heights
defined, we can now store all entries of below its skyline
consecutively in . Fig. shows how
the entries of are stored in . In addition to , we also define an array *MAX *, which stores the
indexes of the diagonal entries of , , in . Note that *MAX *(*i*) is equal to the sum of column
heights up to the (*i*-1)th column plus *i*. Hence, the number of
stored entries in the *i*th column of is equal to *
MAX *(*i*+1) - *MAX *(*i*) and the entry indexes
are *MAX *(*i*), *MAX *(*i*)+1, *
MAX *(*i*)+2, ..., *MAX *(*i*+1) - 1. Therefore,
storing in the form of and *MAX *
allows the entries of to be easily addressed.

A factorization solution of Eq. using skyline storage
is substantially more efficient compared to a Gaussian elimination
solution which has complexity in terms of the number of
operations required (an operation is defined as a multiplication and
an addition), where *n* is the order of . For factorization
with skyline storage, the number of operations required is
in the worst case of constant column heights
(i.e. for all *i*). Clearly, the larger the
difference between *n* and , the greater the advantage. In
simulating the dynamics of the artificial fish, we have *n*=23 and
which represents a significant improvement in
computational efficiency. Note that the column heights of a
system matrix and, hence, the effectiveness of the storage
scheme, depends on how the nodes in the discrete dynamic model are
numbered. For example, using our current numbering scheme, (Fig. ), however if we number node 22 as node
1 (referring to Fig. ), . Thus
in general, one should number the nodes in such a way that all
connections are as local as possible. We can frequently determine a
reasonable nodal point numbering by inspection. However, this
numbering may not be particularly easy to generate in certain cases
[Bathe and Wilson1976].

Xiaoyuan Tu | January 1996 |