Section: Block Matrices

Section: Block Matrices#

This program demonstrates using the Forall construct operate on a block matrix in Fortran.#

The Fortran program block_matrix is designed to generate a block diagonal matrix, where each block along the diagonal is filled with ones, and the rest of the matrix is filled with zeros. The size of the matrix and the blocks within it can be specified by command-line arguments or default to a predefined size if no arguments are provided. Here’s an in-depth explanation of the components and functionalities of the program:

Key Components of the Program#

Variable Declarations
- matrix_size: The size of the square matrix.
- block_size: The size of the blocks that will be filled with ones along the diagonal.
- data: A dynamically allocated 2D integer array that holds the matrix data.
Flow of Execution
- The program starts by calling get_parameters to set the matrix_size and block_size, either from command-line arguments or default values.
- It then creates the block matrix using the create_block_matrix function.
- The matrix is printed using print_matrix.
- Finally, the allocated matrix is deallocated.

Detailed Explanation of Subroutines and Functions#

`get_parameters(matrix_size, block_size)`#

Purpose: To fetch and validate the size of the matrix and blocks from command-line arguments.
Input Handling: If two command-line arguments are provided, they are read and assigned to matrix_size and block_size. If the reading fails (detected by iostat), an error message is displayed.
Defaults: If not exactly two command-line arguments are given, defaults to a matrix size of 9 and a block size of 3.

`create_block_matrix(matrix_size, block_size)`#

Purpose: To generate a block diagonal matrix based on specified sizes.
Process:
- The matrix B is allocated with dimensions matrix_size x matrix_size.
- A nested forall construct is used to fill the diagonal blocks with ones. The construct iterates over starting points of blocks (variable i) and then fills a block_size x block_size block from each starting point (i, i).
- The conditional within the forall ensures that the filling does not exceed the matrix dimensions, avoiding out-of-bound errors.
Error Handling: If allocation fails, an error message is printed and the program is stopped.

`print_matrix(A)`#

Purpose: To print the matrix to the standard output.
Details: Iterates over each row of the matrix A and prints the elements formatted as integers.

Special Features and Techniques Used#

Command Argument Processing: Demonstrates basic usage of get_command_argument for customizable program behavior.
Error Handling in Allocation: Uses stat and error_unit to manage and respond to allocation failures.
Matrix Construction with forall: Efficiently constructs block matrices using a forall statement, a parallel loop construct suited for operations on arrays.

Sample Output Explanation#

Given a matrix_size of 9 and a block_size of 3, the matrix B would look like this:

Each 1 represents a filled element within a 3x3 block along the diagonal in a 9x9 matrix.

Usage Scenario#

This program could be useful in numerical simulations or mathematical modeling where block diagonal matrices are required, such as in block iterative methods for solving large sparse systems of equations.

Program Code#

In file section_block_matrices.f90

program block_matrix
    implicit none
    integer :: matrix_size, block_size
    integer, dimension(:, :), allocatable :: data

    call get_parameters(matrix_size, block_size)
    data = create_block_matrix(matrix_size, block_size) 
    call print_matrix(data)
    deallocate (data)

contains

    subroutine get_parameters(matrix_size, block_size)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        integer, intent(out) :: matrix_size, block_size
        character(len=1024) :: buffer, msg
        integer :: status

        if (command_argument_count() == 2) then
            call get_command_argument(1, buffer)
            read (buffer, fmt=*, iostat=status, iomsg=msg) matrix_size
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') 'error: ', msg
            end if
            call get_command_argument(2, buffer)
            read (buffer, fmt=*, iostat=status, iomsg=msg) block_size
            if (status /= 0) then
                write (unit=error_unit, fmt='(2A)') 'error: ', msg
            end if
        else
            matrix_size = 9
            block_size = 3
        end if
    end subroutine get_parameters

    function create_block_matrix(matrix_size, block_size) result(B)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        integer, value :: matrix_size, block_size
        integer, dimension(:, :), allocatable :: B
        integer :: i, j, k, status

        allocate (B(matrix_size, matrix_size), stat=status)
        if (status /= 0) then
            write (unit=error_unit, fmt='(2A)') &
                'error: can not allocate matrix'
            stop 2
        end if
        B = 0
        forall (i = 1:size(B, 1):block_size, &
                j = 0:block_size - 1, &
                k = 0:block_size - 1, &
                i + j <= size(B, 1) .and. i + k <= size(B, 2))
            B(i + j, i + k) = 1
        end forall
    end function create_block_matrix

    subroutine print_matrix(A)
        implicit none
        integer, dimension(:, :), intent(in) :: A
        integer :: i

        do i = 1, size(A, 1)
            print '(*(I3))', A(i, :)
        end do
    end subroutine print_matrix

end program block_matrix

The above program is compiled and run using Fortran Package Manager (fpm):

Build the Program using FPM (Fortran Package Manager)#

import os
root_dir = ""
root_dir = os.getcwd()

Since the code makes use of the LAPACK library, the following FPM configuration file (fpm.toml) was used:

name = "Section_Block_Matrices"

[build]
auto-executables = true
auto-tests = true
auto-examples = true
module-naming = false

[install]
library = false

[fortran]
implicit-typing = false
implicit-external = false
source-form = "free"

[[executable]]
name="Section_Block_Matrices"
source-dir="app"
main="section_block_matrices.f90"

code_dir = root_dir + "/" + "Fortran_Code/Section_Block_Matrices/"

os.chdir(code_dir)

build_status = os.system("fpm build 2>/dev/null")

Run the Program using FPM (Fortran Package Manager)#

exec_status = \
    os.system("fpm run 2>/dev/null")