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:

  1. 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.

  2. 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.

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

  4. 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.

  1. 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

  2. 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

  3. 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:

  1. 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.

  2. 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

  • rspatialDiscrVhdescribes 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.