# Matrices/Vectors

There exist two types of vectors and matrices, "scalar" vectors/matrices and "block" vectors/matrices. "Scalar" matrices/vectors correspond to `t_spatialDiscretisation`

structures, while "block" matrices correspond to `t_blockDiscretisation`

structures:

`t_vectorScalar`

(linearsystemscalar.f90)- Encapsules all degrees of freedom of a finite element space $V_h$
- Is associated to a finite element space, identified by a
`t_scalarDiscretisation`

structure.

`t_matrixScalar`

(linearsystemscalar.f90)- Describes a linear operator $A_h: V_h \to W_h$, where $V_h$ and $W_h$ are associated finite element spaces.
- Is associated to two finite element spaces $V_h$ and $W_h$, identified by corresponding
`t_scalarDiscretisation`

structures. $V_h$ is referred to as "trial space" and $W_h$ as "test space".

**Note:**$A_h$ is associated to the trial space $V_h$ and the test space $W_h$. The target space is actually $W_h^*$ which is, however, associated to $W_h$ because the space is finite dimensional.`t_vectorBlock`

(linearsystemblock.f90)- Encapsules all degrees of freedom of a tensor product space $X_h$, e.g., $X_h = V_h \times V_h$.
- Is associated to a tensor product space, identified by a
`t_blockDiscretisation`

structure. - Contains a list of
`t_vectorScalar`

structures, each describing a subvector in the block vector. - Subvectors can be reached by
`(t_vectorBlock)%RvectorBlock(:)`

. - By default, memory for a block vector is allocated as a continuous 1D array. The subvectors are realised as parts of the complete vector.

`t_matrixBlock`

(linearsystemblock.f90)- Describes an operator $G_h : X_h \to Y_h$, where $X_h$ and $Y_h$ are tensor products of finite element spaces.
- Is associated to two tensor product spaces $X_h$ and $Y_h$, identified by corresponding
`t_blockDiscretisation`

structures. $X_h$ is referred to as "trial space" and $Y_h$ as "test space".

Note:$G_h$ is associated to the trial space $X_h$ and the test space $Y_h$. The target space is actually $Y_h^*$ which is, however, associated to $Y_h$ because the space is finite dimensional.

Contents | |
---|---|

1 | Working with "scalar" vectors |

2 | Working with "block" vectors |

3 | Working with "scalar" matrices |

4 | Working with "block" matrices |

## Working with "scalar" vectors (linearsystemscalar.f90)

### Essential properties

The following properties are the essential standard properties of a "scalar" vector structure `t_vectorScalar`

from the viewpoint of the user:

Property | Content |
---|---|

NEQ | Number of equations |

p_rspatialDiscr | Pointer an underlying scalar discretisation structure defining an associated finite element space |

There is a 1D data array associated to a vector. By default, this array is of type *double precision*. If one wants to manually modify the content of a vector, a pointer to the content can be obtained using the `lsyssc_getbase_XXXX`

routine, e.g.,

```
type(t_vectorScalar) :: rx
real(DP), dimension(:), pointer :: p_Dx
...
call lsyssc_getbase_double (rx,p_Dx)
! Do something with p_Dx...
```

### Creation/destruction

Vectors can be created based on a discretisation structure. "Scalar" vectors are created based on a "scalar" discretisation structure using `lsyssc_createVector`

and destroyed using `lsyssc_releaseVector`

:

```
type(t_scalarDiscretisation) :: rdiscrVh
type(t_vectorScalar) :: rx
...
call lsyssc_createVector (rdiscrVh,rx)
...
call lsyssc_releaseVector (rx)
```

By default, the vector is not initialised. To initialise a vector with zero upon creation, the additional optional flag `bclear`

can be set to `.true.`

:

```
...
call lsyssc_createVector (rdiscrVh,rx,.true.)
```

### Linear algebra

For the linear algebra, the module `linearsystemscalar.f90`

provides all necessary functionality. `linearsystemscalar.f90`

is an extension of the module `linearalgebra.f90`

which provides the link to high performance BLAS routines. The following table lists a subset of important subroutines for the linear algebra of scalar vectors:

Subroutine | Purpose |
---|---|

lsyssc_clearVector | Clear a vector, i.e. overwrites all entries with 0.0 or with a defined value |

lsyssc_copyVector | Copy a vector to another one |

lsyssc_scaleVector | Scale a vector by a constant |

lsyssc_addConstant | Adds a constant to a vector |

lsyssc_vectorLinearComb | Linear combination of two vectors ("DAXPY") |

lsyssc_scalarProduct | Calculate the scalar product of two vectors |

lsyssc_vectorNorm | Calculate the norm of a vector |

## Working with "block" vectors (linearsystemblock.f90)

### Essential properties

The following properties are the essential standard properties of a "block" vector structure `t_vectorScalar`

from the viewpoint of the user:

Property | Content |
---|---|

NEQ | Total number of equations |

p_rblockDiscr | Pointer an underlying scalar discretisation structure defining an associated finite element space |

nblocks | Number of blocks |

RvectorBlock(:) | Array of subvectors |

There is a 1D data array associated to a vector. By default, this array is of type *double precision*. If one wants to manually modify the content of a vector, a pointer to the content can be obtained using the `lsysbl_getbase_XXXX`

routine, e.g.,

```
type(t_vectorBlock) :: rx
real(DP), dimension(:), pointer :: p_Dx
...
call lsysbl_getbase_double (rx,p_Dx)
! Do something with p_Dx...
```

The data array is of length `NEQ`

and contains the data of all subvectors. If the data of a subvector has to be accessed, this can be done via access to the subvector using `RvectorBlock`

. The following code demonstrates the access to the data of the 3rd subvector:

```
type(t_vectorBlock) :: rx
real(DP), dimension(:), pointer :: p_Dx
...
! Get a pointer to the data of the 3rd subvector
call lsyssc_getbase_double (rx%RvectorBlock(3),p_Dx)
```

### Creation/destruction

Similarly to "scalar" vectors, "block" vectors are created based on a "block" discretisation structure using `lsysbl_createVector`

and destroyed using `lsysbl_releaseVector`

:

```
type(t_blockDiscretisation) :: rdiscrXh
type(t_vectorBlock) :: rx
...
call lsysbl_createVector (rdiscrXh,rx)
...
call lsysbl_releaseVector (rx)
```

The vector is only initialised with zero if the additional flag `bclear=.true.`

is specified:

```
...
call lsysbl_createVector (rdiscrXh,rx,true.)
```

### Linear algebra

For the linear algebra, the module `linearsystemblock.f90`

provides all necessary functionality. `linearsystemblock.f90`

boils down all operations to operations on the subvectors or applies BLAS routines directly if possible. The following table lists a subset of important subroutines for the linear algebra of scalar vectors:

Subroutine | Purpose |
---|---|

lsysbl_clearVector | Clear a vector, i.e. overwrites all entries with 0.0 or with a defined value |

lsysbl_copyVector | Copy a vector to another one |

lsysbl_scaleVector | Scale a vector by a constant |

lsysbl_vectorLinearComb | Linear combination of two vectors ("DAXPY") |

lsysbl_scalarProduct | Calculate the scalar product of two vectors |

lsysbl_vectorNorm | Calculate the norm of a vector |

lsysbl_vectorNormBlock | Calculates the norms of all subvectors |

## Working with "scalar" matrices (linearsystemscalar.f90)

### Essential properties

The following properties are the essential standard properties of a "scalar" matrix structure `t_matrixScalar`

from the viewpoint of the user:

Property | Content |
---|---|

cmatrixFormat | Format of the matrix |

NEQ | Number of rows in the matrix |

NCOLS | Number of columns in the matrix |

dscaleFactor | A scaling factor of the matrix; default=1.0 |

p_rspatialDiscrTrial | Pointer to the trial space |

p_rspatialDiscrTest | Pointer to the test space |

Matrices can be saved in different matrix formats (identified by `cmatrixFormat`

in the structure). Matrix formats are referred to by appropriate IDs. The following table gives an overview about the most commonly used data formats in FEAT2:

Id | Matrix format |
---|---|

LSYSSC_MATRIX1 | Full matrix, saves all entries |

LSYSSC_MATRIX9 | Standard CSR matrix + diagonal pointer |

LSYSSC_MATRIXD | Diagonal matrix, saves only the main diagonal |

### Matrix formats

Each matrix format has its own properties realised in the structure `t_matrixScalar`

. The following paragraphs give a short overview about the matrix formats.

**LSYSSC_MATRIX1 - Full matrix**This matrix format saves all entries in the matrix in one large continuous array of size

`(t_matrixScalar)%NA=NEQ*NCOLS`

. By default, the matrix data is saved in double precision. A pointer to the matrix data can be obtained using the`lsyssc_getbase_XXXX`

routines. The data is stored "row-wise" (similar to C/C++) to speed up matrix-vector multiplication.The follwing example demonstrates how to fill row 50 with "ones":

`type(t_matrixScalar) :: rmatrix real(DP), dimension(:), pointer :: p_Da ... ! Get the data call lsyssc_getbase_double (rmatrix,p_Da) ! Full row 50 with one. do i=1,rmatrix%NCOLS p_Da( (50-1)*rmatrix%NCOLS + i ) = 1.0_DP end do`

**LSYSSC_MATRIX9 - CSR matrix**This matrix format saves only the nonzero entries in a matrix in one continuous array of size

`(t_matrixScalar)%NA`

. By default, the matrix data is saved in double precision. The following arrays are associated to such a matrix:Array name Content `DA`

data array of size `NA`

with the nonzero entries`KCOL`

integer array of size `NA`

; saves for every entry the column in the matrix`KLD`

integer array of size `NEQ+1`

; saves for every row the start index in`DA`

/`KCOL`

; furthermore, there is`KLD(NEQ+1)=NA+1`

`KDIAGONAL`

integer array of size NEQ; saves for every row an index to the diagonal entry in `DA`

/`KCOL`

(;this is an extension to the usual CSR structure); if there is no diagonal entry,`KDIAGONAL`

points to the beginning of the next row, similar to`KLD`

A pointer to the data arrays can be obtaind using the

`lsyssc_getbase_XXXX`

routines. The follwing example demonstrates how to set the entry at row 50 / column 40 to 1 (assuming that A(40,50) exists):`type(t_matrixScalar) :: rmatrix real(DP), dimension(:), pointer :: p_Da integer, dimension(:), pointer :: p_Kcol integer, dimension(:), pointer :: p_Kld ... ! Get the data call lsyssc_getbase_double (rmatrix,p_Da) call lsyssc_getbase_Kcol (rmatrix,p_Kcol) call lsyssc_getbase_Kld (rmatrix,p_Kld) ! Set A(50,40)=1. Search for column 40 in row 50. do i=p_Kld(50),p_Kld(50+1)-1 if (p_Kcol(i) == 40) then p_Da(i) = 1.0_DP exit end if end do`

**LSYSSC_MATRIXD - Diagonal matrix**This matrix format saves only the entries of the main diagonal in one continuous array of size

`(t_matrixScalar)%NA=NEQ`

. By default, the matrix data is saved in double precision. A pointer to the matrix data can be obtained using the`lsyssc_getbase_XXXX`

routines.The follwing example demonstrates how to fill the main diagonal with "ones":

`type(t_matrixScalar) :: rmatrix real(DP), dimension(:), pointer :: p_Da ... ! Get the data call lsyssc_getbase_double (rmatrix,p_Da) ! Fill the main diagonal with "(1,1,1,1,...)" do i=1,rmatrix%NEQ p_Da(i) = 1.0_DP end do`

### Content and structure

A matrix structure `t_matrixScalar`

encapsules two types of data:

**Matrix structure**Data of type "matrix structure" describe structural data of the matrix. Apart of simple data like NEQ, this term in particular refers to the column/row structure of CSR matrix (

`KCOL`

,`KLD`

), i.e., the pattern of zero/nonzero entries.**Matrix content**The term "matrix content" refers to the actual entries of the matrix which are calculated, e.g., during the integration.

FEAT2 differs between these two types of data as they are logically completely different and may be used independent of each other. Matrices always have a "structure" but not necessarily a content:

- For simple matrices like full matrices of diagonal matrices, the "structure" contains only simple data like
`NEQ`

and does not need memory. - For more complicated matrices like CSR, the "structure" is the nonzero pattern (which needs memory) while the "content" is the actual matrix data.
- "Graph matrices" which realise a graph only contain structural data and no content.
- This includes, in particular, the structure of finite element matrices. Different matrices for the same finite element space $V_h$ all have the same structure, so the structure can be calculated once in advance.

To work with the structure and the content, the module `linearsystemscalar.f90`

provides a couple of subroutines:

Subroutine | Purpose |
---|---|

lsyssc_convertMatrix | Convert a matrix into another matrix structure |

lsyssc_hasMatrixStructure | Check if a matrix has a structure in memory or not |

lsyssc_hasMatrixContent | Check if a matrix has a content in memory or not |

lsyssc_releaseMatrixContent | Releases the content of the matrix, the structure will stay unchanged |

### Creation/destruction of scalar matrices

Depending on the type of the matrix, a "scalar" matrix is created in one or two steps, depending on the type.

**Full matrix**Full matrices are always created in matrix format LSYSSC_MATRIX1. This type of matrix is created using the subroutine

`lsyssc_createFullMatrix`

from`linearsystemscalar.f90`

. Example:`type(t_spatialDiscretisation) :: rspatialDiscr type(t_matrixScalar) :: rmatrix ... call lsyssc_createFullMatrix (rspatialDiscr,rmatrix [,bclear=.true.])`

where the optional parameter

`bclear=.true.`

can be specified to initialise the matrix with zero upon creation. FEAT2 reserves memory for the entries (content) on the heap.Alternatively, the matrix can be created by directly specifying the size of the matrix (without discretisation structure attached). The following example creates an

`NEQ x NEQ`

matrix without attached discretisation structure:`type(t_matrixScalar) :: rmatrix ... call lsyssc_createFullMatrix (rmatrix,NEQ [,bclear=.true.])`

**Diagonal matrix**Diagonal matrices can be created in all matrix formats using the subroutine

`lsyssc_createDiagMatrixStruc`

from`linearsystemscalar.f90`

. The matrix format has to be specified as parameter. Example:`type(t_spatialDiscretisation) :: rspatialDiscr type(t_matrixScalar) :: rmatrix ... call lsyssc_createDiagMatrixStruc (& rspatialDiscr,LSYSSC_MATRIXD,rmatrix [,bclear=.true.])`

where the optional parameter

`bclear=.true.`

can be specified to initialise the matrix with zero upon creation. FEAT2 reserves memory for the entries (content) on the heap.Alternatively, if the matrix is square, it can be created by directly specifying the size of the matrix (without discretisation structure attached). The following example creates a diagonal

`NEQ x NEQ`

matrix without attached discretisation structure:`type(t_matrixScalar) :: rmatrix ... call lsyssc_createDiagMatrixStruc (& rmatrix,LSYSSC_MATRIXD,NEQ [,bclear=.true.])`

**CSR matrices**For the creation of CSR matrices, special routines from

`bilinearformevaluation.f90`

are designed to determine the matrix pattern based on the underlying finite element space. The creation of a CSR matrix requires two steps: Creation of the structure and creation of the content. The main entry point is the subroutine`bilf_createMatrixStructure`

. The following example demonstrates how to create a CSR matrix ("format 9") based on a discretisation:`type(t_spatialDiscretisation) :: rspatialDiscr type(t_matrixScalar) :: rmatrix ... call bilf_createMatrixStructure (& rspatialDiscr,LSYSSC_MATRIX9,rmatrix)`

CSR matrices handle structure and content independent of each other. After a call to

`lsyssc_createDiagMatrixStruc`

, the "content" of the matrix is empty, no data arrays are attached to the matrix. To create a data array for the entries, one has to use the additional call`lsyssc_allocEmptyMatrix`

:`call lsyssc_allocEmptyMatrix (rmatrix [,bclear=.true.])`

The optional parameter

`bclear=.true.`

can be specified to initialise the entries with zero upon creation.

Independent of the creation, a matrix is destroyed using `lsyssc_releaseMatrix`

:

```
call lsyssc_releaseMatrix (rmatrix)
```

### Different trial and test spaces

The matrix creation routines usually allow to specify different trial and test spaces which usually leads to a rectangular matrix with more or less columns than rows. Such situations usually appear during the creation of submatrices in a block matrix, see below. The alternative test space can be specified as additional parameter during the creation of the matrix structure.

The following example creates a CSR matrix with different trial and test spaces. The trial space is identified by `rspatialDiscrVh`

and the test space as `rspatialDiscrWh`

:

```
type(t_spatialDiscretisation) :: rspatialDiscrVh
type(t_spatialDiscretisation) :: rspatialDiscrWh
type(t_matrixScalar) :: rmatrix
...
call bilf_createMatrixStructure (&
rspatialDiscrVh,LSYSSC_MATRIX9,rmatrix,rspatialDiscrWh)
```

### Linear algebra

For the linear algebra, the module `linearsystemscalar.f90`

provides all necessary functionality. The following table lists a subset of important subroutines for the linear algebra of scalar matrices:

Subroutine | Purpose |
---|---|

lsyssc_matVec | Matrix-vector multiplication |

lsyssc_invertedDiagMatVec | Multiply a vector with the inverse of the diagonal of a scalar matrix (-> used, e.g., for the Jacobi iteration) |

lsysbl_allocEmptyMatrix | Allocates memory for the entries of a matrix |

lsyssc_clearMatrix | Clear a matrix |

lsyssc_scaleMatrix | Scale a matrix by a constant |

lsyssc_matrixLinearComb | Linear combination of two matrices |

lsyssc_calcDeterminant | Calculates the determinant of a square matrix |

lsyssc_calcGerschgorin | Apply the theorem of Gerschgorin to get approximations for the min/max eigenvalues of a matrix |

lsyssc_clearOffdiags | Clear all offdiagonal entries in a matrix |

lsyssc_multMatMat | Multiplies two matrices |

lsyssc_transposeMatrix | Calculate the transpose of a matrix |

lsyssc_lumpMatrix | Lumping of a matrix |

lsyssc_initialiseIdentityMatrix | Initialises the content of a matrix to an identity matrix |

lsyssc_copyMatrix | Copies matrix structure and content |

### Sharing/copying of matrix data

Depending on the underlying problem, finite element space etc., finite element matrices can be tremendeously large. FEAT2 implements the special functionality of "shared" matrix data to save memory. Every matrix can "share" its structure and/or content with another "parent" matrix. Memory is allocated and maintained only in this "parent" matrix, while the "child" matrix does not need additional memory for the shared data. This technique allows to work with so-called "template matrices" and is a good technique to ensure that matrices have a common nonzero pattern. Whereever one wants to "derive" one matrix or another, it should be taken into account whether the data should actually be copied or shared.

There are three main routines in `linearsystemscalar.f90`

responsible for the implementation of shared matrices:

Subroutine | Purpose |
---|---|

lsyssc_duplicateMatrix | Extended version of `lsyssc_copyMatrix` |

lsyssc_isMatrixStructureShared | Tests if the structure of a matrix is shared with another matrix |

lsyssc_isMatrixContentShared | Tests if the content of a matrix is shared with another matrix |

The subroutine `lsyssc_duplicateMatrix`

is an extended and much cheaper variant of `lsyssc_copyMatrix`

and should be used whereever possible. It allows to selectively copy or share the data of one matrix with another. The routine accepts four parameters:

- A source matrix.
- A destination matrix.
- What to do with the structure of the source.
- What to do with the content of the source.

Let us demonstrate this on an example with template matrices. The following example at first create a template finite element matrix in CSR format. This template matrix only has a structure defining the nonzero pattern of the entries. No content is attached. We create it only for the purpose to create other matrices later:

```
type(t_spatialDiscretisation) :: rspatialDiscr
type(t_matrixScalar) :: rmatrixTemplate
type(t_matrixScalar) :: rmatrixLaplace, rmatrixMass
...
call bilf_createMatrixStructure (rspatialDiscr,LSYSSC_MATRIX9,rmatrixTemplate)
```

In a second step, we use `lsyssc_duplicateMatrix`

to create a matrix structure for a Laplace and a mass matrix. Both matrices will "share" the matrix structure (so `KCOL`

, `KLD`

,...) with the template matrix, this does not need additional memory. The content of the two matrices is different, of course, so we tell `lsyssc_duplicateMatrix`

to allocate new memory for the content, which we initialise by zero in a second step:

```
....
call lsyssc_duplicateMatrix (rmatrixTemplate,rmatrixLaplace,&
LSYSSC_DUP_SHARE,LSYSSC_DUP_EMPTY)
call lsyssc_duplicateMatrix (rmatrixTemplate,rmatrixMass,&
LSYSSC_DUP_SHARE,LSYSSC_DUP_EMPTY)
call lsyssc_clearMatrix (rmatrixLaplace)
call lsyssc_clearMatrix (rmatrixMass)
...
```

As can be seen, the constants `LSYSSC_DUP_xxxx`

define what to do with the structure (here: "share") and the content (here: allocate a new, "empty" content). The follwing list shows the effect of the different constants if applied to structure and/or content. For this table, assume a call

```
call lsyssc_duplicateMatrix (rA,rB,LSYSSC_DUP_...,LSYSSC_DUP_...)
```

Constant | structure/content(rB) |
---|---|

LSYSSC_DUP_SHARE | shared with structure/content(rA) |

LSYSSC_DUP_EMPTY | allocate new memory, no initialisation |

LSYSSC_DUP_COPYOVERWRITE | overwrite with structure/content(rA) |

LSYSSC_DUP_COPY | copy from structure/content(rA), allocate new memory |

LSYSSC_DUP_ASIS | copy if structure/content(rB) is not shared, ignore if shared |

LSYSSC_DUP_REMOVE | remove any old structure/content(rB), deallocated if necessary |

LSYSSC_DUP_DISMISS | drop any old structure/content(rB), no deallocation |

Warning:If a "child" matrix shares its date with a parent matrix, the user is responsible to take care of the data of the parent matrix. Changing data arrays of the parent matrix would also change the data arrays of the child matrix.

Child matrices have to be released with `call lsyssc_releaseMartix`

in the same way as other matrices. FEAT2 takes care that releasing a child matrix will not to a destruction of the data in the parent matrix.

Warning:The opposite is not true. If the parent matrix is destroyed, the child matrix gets invalid. If the parent matrix has been released, the user shall never use a child matrix for anything except for destruction.

## Working with "block" matrices (linearsystemblock.f90)

In contrast to "scalar" matrices, "block" matrices do not contain any actual data. A "block" matrix is realised as a (m x n) array of scalar matrices that describe the actual operators between the finite element spaces. It is a logical representation of an operator between the corresponding tensor product spaces.

### Essential properties

The following properties are the essential standard properties of a "block" matrix structure `t_matrixBlock`

from the viewpoint of the user:

Property | Content |
---|---|

NEQ | Total number of rows |

NCOLS | Total number of columns |

nblocksPerRow | Number of "block columns", i.e., number of blocks in one "block row" of the block matrix |

nblocksPerCol | Number of "block rows", i.e., number of blocks in one "block column" of the block matrix |

RmatrixBlock(:,:) | 2D array of all submatrices |

p_rblockDiscrTest | Pointer to the trial space |

p_rblockDiscrTrial | Pointer to the test space |

### Creation/destruction

A block matrix is created upon a "block discretisation" by a simple call to `lsysbl_createMatrix`

and released after use with `lsysbl_releaseMatrix`

:

```
type(t_blockDiscretisation) :: rdiscr
type(t_matrixScalar) :: rmatrix
...
call lsysbl_createMatrix (rdiscr,rmatrix)
...
call lsysbl_releaseMatrix (rmatrix)
```

The command `lsysbl_createMatrix`

creates an "empty" matrix with the number of block rows/columns matching the number of finite element spaces in rdiscr. No memory is allocated for any submatrix, and no submatrix is initialised.

### Working with submatrices

The user can access the submatrices via `(t_matrixBlock)%RmatrixBlock(i,j)`

. This array provides empty memory for all the submatrices, and the user has to create submatrices where necessary with the usual techniques for the structure `t_matrixScalar`

. The user is responsible to initialise the submatrices with the correct "scalar" discretisation structure.

In the following example, it is assumed that

`rspatialDiscrVh`

describes a finite element space $V_h$,`rspatialDiscrWh`

describes a finite element space $W_h$,`rblockDiscr`

describes the tensor product space $X_h$ := $V_h$ x $W_h$, created by`rspatialDiscrVh`

and`rspatialDiscrWh`

.

The example creates a 2x2 block matrix, realising an operator $G_h : X_h \to X_h$. In the diagonal blocks (1,1) and (2,2), empty finite element matrices are created in CSR format and initialised with zero.

```
type(t_scalarDiscretisation) :: rspatialDiscrVh,rspatialDiscrWh
type(t_blockDiscretisation) :: rblockDiscr
type(t_matrixBlock) :: rmatrix
...
call lsysbl_createMatrix (rblockDiscr,rmatrix)
! Create block (1,1)
call bilf_createMatrixStructure (&
rspatialDiscrVh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(1,1))
! Create block (2,2)
call bilf_createMatrixStructure (&
rspatialDiscrWh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(2,2))
! Allocate all created submatrices and initialise with zero
call lsysbl_allocEmptyMatrix(rmatrix,.true.)
```

If, additionally, the matrix at position (1,2) and/or (2,1) should be created, an extended syntax for `bilf_createMatrixStructure`

has to be used. For these matrices trial and test spaces differ, and this has to be taken into account in the call to `bilf_createMatrixStructure`

. The following example creates submatrices at all four blocks:

```
type(t_scalarDiscretisation) :: rspatialDiscrVh,rspatialDiscrWh
type(t_blockDiscretisation) :: rblockDiscr
type(t_matrixBlock) :: rmatrix
...
call lsysbl_createMatrix (rblockDiscr,rmatrix)
! Create block (1,1)
call bilf_createMatrixStructure (&
rspatialDiscrVh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(1,1))
! Create block (2,2)
call bilf_createMatrixStructure (&
rspatialDiscrWh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(2,2))
! Create block (1,2). Offdiagonal block.
call bilf_createMatrixStructure (&
rspatialDiscrVh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(1,2),&
rspatialDiscrWh)
! Create block (2,1). Offdiagonal block.
call bilf_createMatrixStructure (&
rspatialDiscrWh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(2,1),&
rspatialDiscrVh)
! Allocate all created submatrices and initialise with zero
call lsysbl_allocEmptyMatrix(rmatrix,.true.)
```

If some finite element spaces are the same, it is advisable to use template matrices in order to save memory. In the following example, it is assumed that

`rspatialDiscrVh`

describes a finite element space $V_h$,`rblockDiscr`

describes the tensor product space $X_h$ := $V_h$ x $V_h$, created by`rspatialDiscr`

.

The example creates a 2x2 block matrix, realising an operator $G_h : X_h \to X_h$. In the diagonal blocks (1,1) and (2,2), empty finite element matrices are created in CSR format and initialised with zero. The submatrices are created based on a "template" matrix for the finite element space. For this purpose, the routine `lsysbl_duplicateMatrix`

is used which copies a scalar submatrix into a block of a block matrix.

```
type(t_scalarDiscretisation) :: rspatialDiscrVh
type(t_blockDiscretisation) :: rblockDiscr
type(t_matrixBlock) :: rmatrixTemplate
type(t_matrixBlock) :: rmatrix
...
! Create a template matrix for the space $V_h$
call bilf_createMatrixStructure (&
rspatialDiscrVh,LSYSSC_MATRIX9,rmatrix%RmatrixBlock(1,1))
! Create the block matrix
call lsysbl_createMatrix (rblockDiscr,rmatrix)
! Create block (1,1) / (2,2) using the template matrix
call lsysbl_duplicateMatrix (rmatrixTemplate,rmatrix,1,1,&
LSYSSC_DUP_SHARE,LSYSSC_DUP_EMPTY)
call lsysbl_duplicateMatrix (rmatrixTemplate,rmatrix,2,2,&
LSYSSC_DUP_SHARE,LSYSSC_DUP_EMPTY)
! Initialise with zero
call lsysbl_clearMatrix(rmatrix)
```

The following table gives a brief overview about the most important subroutines in `linearsystemblock.f90`

which are used to deal with submatrices:

Subroutine | Purpose |
---|---|

lsysbl_duplicateMatrix | Duplicates a scalar matrix into a submatrix of a block matrix or duplicates a block matrix to another one |

lsysbl_isSubmatrixPresent | Checks if a submatrix at a position (i,j) is present |

lsysbl_releaseMatrixRow | Releases a row of submatrices from a block matrix |

lsysbl_releaseMatrixColumn | Releases a column of submatrices from a block matrix |

### Linear algebra

For the linear algebra with block matrices, the module `linearsystemblock.f90`

provides a set of subroutines. The most important subroutines can be found in the following table:

Subroutine | Purpose |
---|---|

lsysbl_matVec | Matrix-vector multiplication |

lsysbl_invertedDiagMatVec | Multiply a vector with the inverse of the diagonal of a matrix |

lsysbl_allocEmptyMatrix | Allocates memory for the entries of a matrix |

lsysbl_clearMatrix | Clear a matrix |

lsysbl_scaleMatrix | Scale a matrix by a constant |

lsysbl_copyMatrix | Copies matrix structural data and content |

lsysbl_duplicateMatrix | Extended version of `lsysbl_copyMatrix` . Applies `lsyssc_duplicateMatrix` to every submatrix or duplicates a scalar matrix into a submatrix of a block matrix. |

Other functionality can be emulated by the user by applying the corresponding linear algebra routines to the scalar subroutines.