Next: Algorithm Outline and Discussion Up: Numerical Solution Previous: Numerical Solution

## System Matrix Assembling and the Skyline Storage Scheme

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 :

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

and assemble the term of Eq. () into :

3. 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 ith 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].

Next: Algorithm Outline and Discussion Up: Numerical Solution Previous: Numerical Solution
 Xiaoyuan Tu January 1996