Next: 1.3 Solution by an iterative method Up: 1 Solution of linear systems Prev: 1.1 Introduction Contents

1.2 Solution by a direct method

1.2.1 Introduction

The methods presented in this chapter are the three classical solution methods for a linear system by matrix factorization . We can therefore distinguish between:

• the Cholesky  method, is used when matrix A of the linear system to be solved is symmetric positive definite    : , where L is the lower triangular matrix.
• the Crout method , is used when the matrix is symmetric  : , where L is the lower triangular matrix with unit diagonal, and D is the diagonal matrix.
• the Gauss  method, is used for all non-singular  matrices : A=L U, where L is the lower triangular matrix with unit diagonal, and U is the upper triangular matrix.

Several versions of these three methods are implemented in the MODULEF code, according to the memory space necessary to store matrix A.

To start with, we assume that matrix A can be stored entirely in main memory, so that the storage of the matrix and the right-hand-side vector can be described exactly.

1.2.2 The MUA data structure

For more details consult [MODULEF User Guide - 2].

Contents

This DS stores, for each line, the coefficients lying between the first column with a a priori non-zero coefficient and the diagonal of a sparse finite element matrix. This type of storage is called profile or skyline storage.

DS MUA consists of 6 arrays of predefined order.

Array MUA0:

General information.

This integer array contains 32 variables, consisting of a general description of the job (title, date, name), of DS MUA (type, level, ...), and indicates the presence or absence of array MUA1.

1:20 TITRE

the job title in 20 words of 4 characters,

21:22 DATE

the date of creation in 2 words of 4 characters,

23:28 NOMCRE

the creator's name in 6 words of 4 characters,

29 'MUA '

the DS type,

30 NIVEAU

the DS level,

31 ETAT

a reserved parameter,

32 NTACM

the number of supplementary arrays associated with the DS (they are described in array MUA1).

Array MUA1:

Description of any supplementary arrays.

This array is analogous to array B1 of DS B (see this DS).

Array MUA2:

General description of the matrix.

This integer array contains 12 values.

1 NTYP

the type of matrix coefficients,

2 NBLOC

the number of blocks or pages in the matrix,

3 NTCOL

the number of columns of the largest block,

4 NMOPB

the number of words required in main memory to store the largest block,

5 NMATA

the number of matrices on file,

6 NCODSA

the type of matrix storage:

• 1: symmetric matrix, only the lower triangle is stored, line by line,
• 0: diagonal matrix, or
• -1: non-symmetric matrix, ordered line after line, with the diagonal coefficient placed at the end.
  

7 LBDP

the largest difference + 1 between the 2 node numbers in the same element if NCODSA is non-zero, 1 if NCODSA=0 i.e. the half band-width in terms of nodes,

8 LBDPDL

the largest difference + 1 between the degree of freedom numbers of 2 nodes belonging to the same element if NCODSA is non-zero, 1 if NCODSA=0 i.e. the half band-width in terms of degrees of freedom,

9 ND

the number of degrees of freedom per node if it is constant, 0 otherwise,

10 NOE

the number of nodes,

11 NTDL

the matrix order,

12 MUA5AR

the assembly is already performed if 1,
the assembly is not yet performed if 0.

Array MUA3:

Pointer to the last degree of freedom of each block.

This integer array , of length NBLOC+1, contains:

• MUA3(1) = 0,
• MUA3(i+1) = the number of the last degree of freedom of block i in the matrix.

Array MUA4:

Pointer to the diagonal coefficient of each line.

This integer array, of length NTDL+1 or 2 depending on NCODSA, contains:

• NCODSA = 0
  • MUA4(1) = 0
  • MUA4(2) = NTDL,
• NCODSA 0
  • MUA4(1) = 0
  • MUA4(i+1) = the address of diagonal coefficient i.

Array MUA5:

The matrix coefficients.

This NTYP type array contains the matrix coefficients:

• MUA5(k) = value of the k-th matrix coefficient (if it is not paginated),
• or MUA5(k-MUA3(i)) the k-th coefficient of block i (if it is paginated).

Example 1: Symmetric matrix

Consider the 12-th order matrix given below (the values correspond to the the rows in the storage and not to the matrix coefficients) constituting 3 pages:

Array MUA5, of length MUA4(13), is subdivided into 3 pages, line 1 to 5, 6 to 10 and 11 to 12. Each page contains a maximum of 32 words. Thus, we have:

• NCODSA = 1
• array MUA3:
  

• array MUA4:
  

Example 2: Non-symmetric matrix

Consider the 12-th order matrix given below (the values correspond to the the rows in the storage and not to the matrix coefficients) constituting 3 pages:

Here, each page contains a maximum of 12 words. Thus, we have:

• NCODSA = -1
• array MUA3:
  

• array MUA4:
  

For the symmetric case we store the lower triangular part of matrix A in a one-dimensional array MUA5 , as follows: for each variable i,we store all the coefficients of the row, for , (where mc(i) is the smallest column index such that ). These coefficients (zero or not) are ordered by increasing column index. The rows are then ordered in order of increasing variable in array MUA5. To maintain the correspondence between MUA5 and A, a pointer MUA4  is used which gives the address of the diagonal coefficient in MUA5 for each row i. The address ia in MUA5 of an arbitrary coefficient is therefore calculated by the formula ( With the convention MUA4(1)=0)

This storage mode, called profile or skyline storage, is particularly useful as it is also used for matrix L. In fact, the matrices A and L have the same profile, and in practice they are often stored in the same array).

For the non-symmetric case , the same storage is used by adding the columns in the following manner: for each variable i, we store all the coefficients of the row, for , (where mc(i) is the smallest column index such that ). These coefficients (zero or not) are ordered by increasing column index; then we store all the coefficients of the column, for , ending with the diagonal coefficient. This information is then ordered in order of increasing variable in array MUA5. To maintain the correspondence between MUA5 and A, a pointer MUA4 is used which gives the address of the diagonal coefficient in MUA5 for each row i. The address ia in MUA5 of an arbitrary coefficient is therefore calculated by the formula

1.2.3 The B data structure

For more details consult [MODULEF User Guide - 2].

Contents

This DS stores a coefficient array in the form of a vector (right-hand-side (RHS) of a system, solution, etc.).

In practice, matrix B has the form B(NDSM,NTDL) where NDSM designates, for example, the number of load cases of RHS and NTDL is the number of degrees of freedom.

DS B consists of five arrays of predefined order.

Array B0:

General information.

This integer array contains 32 variables, consisting of a general description of the job (title, date, name), of DS B (type, level, ...), and indicates the presence or absence of array B1.

1:20 TITRE

the job title in 20 words of 4 characters,

21:22 DATE

the date of creation in 2 words of 4 characters,

23:28 NOMCRE

the creator's name in 6 words of 4 characters,

29 'B '

the DS type,

30 NIVEAU

the DS level,

31 ETAT

a reserved parameter,

32 NTACM

the number of supplementary arrays associated with the DS (they are described in array B1).

Array B1:

Description of any supplementary arrays.

This integer array exists if NTACM (see B0) is non-zero. It contains, for each array to be described, 22 variables. In general, there are no associated arrays, except for certain applications where it is necessary to store values in the DS which do not enter its standard definition.

If array B1 exists, it contains sequentially:

• Loop i from 1 to NTACM
  • array name in 4 characters,
  • array address in the super-array M,
  • number of words of this array,
  • array type (1 = integer, 2 = real single precision, 4 = character, 5 = real double precision, ...)
  • comments corresponding to the contents of the array in 18 words of 4 characters.
• End of loop.

Array B2:

General description of the coefficient array.

This integer array contains 11 values.

1 NTYP

the type of array B4,

2 NIND

the number of indices (fixed to 2),

3 M1

the maximum value of its first index (for example, NDSM, the number of load cases),

4 M2

the maximum value of its second index (for example, NTDL if B4 only has one page),

5 NCOD

the code for splitting array B4 into blocks, with:

• NCOD = 1 if all the blocks have the same size (except the last),
  
• NCOD = 0 if not,
  

6 NBLOC

the number of blocks,

7 NOE

the number of nodes,

8 NBBLOC

the number of B4 arrays when the DS is on file,

9 ND

the number of degrees of freedom per node if it is constant, or 0,

10 NTDL

the number of degrees of freedom or the maximum value of the second index of array B4 when the NBLOC blocks are considered,

11 NCODSB

the storage code for array B4 with:

• NCODSB = -1 for a storage by degrees of freedom (generally, this storage is used),
• NCODSB = -2 for a storage by nodes.

Array B3:

Pointer to the blocks.

This integer array contains NBLOC+1 variables indicating the segmentation of array B4.

• B3(1) = 0
• B3(i+1) = number of the last node in block i if NCODSB = -2
  B3(i+1) = number of the last degree of freedom in block i if NCODSB = -1

Array B4:

the array coefficients.

This NTYP type array contains the NDSM*M2 values of each page.

B4(1,l+1) to B4(NDSM,l+1) then

B4(1,l+2) to B4(NDSM,l+2) then

....

B4(1,l') to B4(NDSM,l')

where
l = NDL(B3(i))+1, l' = NDL(B3(i+1)) where NDL is a function such that NDL(1) = 0 and NDL(i+1) is the last degree of freedom of node i.

This set of 5 arrays (B0, B1, B2, B3 and B4) managed by a few subroutines, contain the NDSM right-hand-sides (or solutions) of the linear system A x = b. The coefficients of array B4(NDSM,NTDL) represent the NTDL  components of the NDSM  right-hand-sides of the linear system to be solved. After the solution procedure, the solutions of the linear system are also stored in array B4 (the same, or one of a distinct structure B).