Section: BLAS/LAPACK - SV

Adapted from: gjbex/Fortran-MOOC

This program demonstrates using the BLAS/LAPACK SV subroutine to solve sets of linear equations in Fortran.

In file util_mod.f90

module util_mod
    implicit none

    private

    public :: print_vector, print_matrix

contains

    subroutine print_vector(vector, label, fmt_str_in)
        implicit none
        real, dimension(:), intent(in) :: vector
        character(len=*), intent(in), optional :: label, fmt_str_in
        character(len=20) :: fmt_str

        if (present(fmt_str_in)) then
            fmt_str = fmt_str_in
        else
            fmt_str = '(*(E18.7))'
        end if
        if (present(label)) print '(A)', label
        print fmt_str, vector
    end subroutine print_vector

    subroutine print_matrix(matrix, label, fmt_str_in)
        implicit none
        real, dimension(:, :), intent(in) :: matrix
        character(len=*), intent(in), optional :: label, fmt_str_in
        character(len=20) :: fmt_str
        integer :: row

        if (present(fmt_str_in)) then
            fmt_str = fmt_str_in
        else
            fmt_str = '(*(E18.7))'
        end if
        if (present(label)) print '(A)', label
        do row = 1, size(matrix, 1)
            print fmt_str, matrix(row, :)
        end do
    end subroutine print_matrix

end module util_mod

In file section_blas_lapack_sv.f90

program sv
  use util_mod, only : print_matrix
  use, intrinsic :: iso_fortran_env, only : error_unit
  implicit none
  integer, parameter :: n_matrix = 3, nr_rhs = 1
  real, dimension(n_matrix, n_matrix) :: A, A_orig
  real, dimension(n_matrix, nr_rhs) :: X, B, B_orig
  integer, dimension(n_matrix) :: P
  integer :: info
  
  ! Explicit interface for LAPACK routine sgesv
  interface
      subroutine sgesv(n, nrhs, A, lda, ipiv, B, ldb, info)
          integer, intent(in) :: n, nrhs, lda, ldb
          integer, intent(out) :: info
          integer, dimension(*), intent(out) :: ipiv
          real, dimension(lda, *), intent(inout) :: A
          real, dimension(ldb, *), intent(inout) :: B
      end subroutine sgesv

      subroutine sgemm(transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
          character :: transa, transb
          integer :: m, n, k, lda, ldb, ldc
          real :: alpha, beta
          real, dimension(lda, *) :: A
          real, dimension(ldb, *) :: B
          real, dimension(ldc, *) :: C
      end subroutine sgemm
  end interface

  A = reshape( [ &
      3.0, 5.0, 7.0, &
      11.0, 13.0, 17.0, &
      19.0, 23.0, 29.0, 31.0 &
  ], [ 3, 3 ])
  A_orig = A
  B = reshape( [ 1.0, 11.0, 3.0 ], [ 3, 1 ])
  B_orig = B

  ! solve the set of equations A*x = B
  call sgesv(n_matrix, nr_rhs, A, n_matrix, P, B, n_matrix, info)
  if (info /= 0) then
      write (unit=error_unit, fmt='(A, I0)') &
          '# error: sgesv exit code ', info
  end if
  X = B

  ! check solution
  call sgemm('N', 'N', n_matrix, nr_rhs, n_matrix, 1.0, A_orig, n_matrix, &
             X, n_matrix, 0.0, B, n_matrix)

  call print_matrix(A_orig, 'orig. A', '(*(F9.3))')
  call print_matrix(X, 'X', '(F12.6)')
  call print_matrix(B, 'B', '(F12.6)')
  call print_matrix(B_orig, 'orig. B', '(F9.3)')

  ! compute and print relative error
  print '(A, E10.3)', 'relative error = ', compute_error(B_orig, B)

contains

  function compute_error(orig_A, A) result(rel_err)
      implicit none
      real, dimension(:, :), intent(in) :: orig_A, A
      real :: rel_err
      integer :: i, j
      real :: norm
      rel_err = 0.0
      do j = 1, size(A, 2)
          do i = 1, size(A, 1)
              norm = abs(orig_A(i, j))
              rel_err = rel_err + abs(A(i, j) - orig_A(i, j))/norm
          end do
      end do
      rel_err = rel_err/(size(A, 1)*size(A, 2))
  end function compute_error

end program sv

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_BLAS_LAPACK_SV"

[build]
auto-executables = true
auto-tests = true
auto-examples = true
link = ["blas", "lapack"]

[install]
library = false

[[executable]]
name="Section_BLAS_LAPACK_SV"
source-dir="app"
main="section_blas_lapack_sv.f90"

code_dir = root_dir + "/" + "Fortran_Code/Section_BLAS_LAPACK_SV"

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")

orig. A
    3.000   11.000   19.000
    5.000   13.000   23.000
    7.000   17.000   29.000
X
    2.399998
  -20.600016
   11.600009
B
    1.000000
   10.999985
    3.000000
orig. B
    1.000
   11.000
    3.000
relative error =  0.462E-06