matrix Module

Matrix Module programmed by Son, Sang-Kil in Feb. 2004 for PHSX815 : Computational Physics

Ref) Press et al., "Numerical Recipes in Fortran 77" Press et al., "Numerical Recipes in Fortran 90"

Main data types - LU_matrix : LU-decomposed matrix - F LU_decompose(A(N,N)) : create LU matrix by LU-decomposition of A - F LU_solve(LU, b(N)) : solve a linear equation for b with LU - F LU_inverse(LU) : inverse matrix of LU - F LU_determinant(LU) : determinant of LU - F LU_lower(LU) : return L matrix of LU - F LU_upper(LU) : return U matrix of LU

Auxiliary functions and subroutines - SR print_matrix(A(N,N)) : print out a matrix - SR print_vector(x(N)) : print out a vector


Uses


Interfaces

public interface operator(.m.)

Operators for efficient matrix multiplication with LAPACK, but better readibility. (It is assumed that the dimensionality fits.)

multiply two real square matrices or a square matrix times a vector AB= A .m. B or Av= A .m. v

  • private function mult_real_matrix_matrix(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:,:)

    Return Value real(kind=long), (size(a,dim=1),size(b,dim=2))

  • private function mult_real_matrix_vector(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:)

    Return Value real(kind=long), (size(a,dim=1))

public interface operator(.mt.)

multiply two real square matrices, the second is transposed A(B*T)= A .mt. B

  • private function mult_real_matrix_trans_matrix(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:,:)

    Return Value real(kind=long), (size(a,dim=1),size(b,dim=1))

public interface operator(.tm.)

multiply two real square matrices, the first is transposed (A*T)B= A .tm. B

  • private function mult_real_trans_matrix_matrix(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:,:)

    Return Value real(kind=long), (size(a,dim=2),size(b,dim=2))

  • private function mult_real_trans_matrix_vector(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:)

    Return Value real(kind=long), (size(a,dim=2))

public interface operator(.sub.)

substract two real square matrices element wise A-B = A .sub. B This functionality is provided by f90. Element wise, arrays (of same dimensionality) can be added, substracted, multiplied and divided. sub_real_matrix_matrix is not a LAPACK feature. (Michael)

  • private function sub_real_matrix_matrix(a, b) result(c)

    Arguments

    Type IntentOptional Attributes Name
    real(kind=long), intent(in) :: a(:,:)
    real(kind=long), intent(in) :: b(:,:)

    Return Value real(kind=long), (size(a,dim=1),size(a,dim=1))


Derived Types

type, public ::  lu_matrix

Data type: LU_matrix it contains the dimension of the matrix, N and LU-decomposed matrix, LU(N,N) also, permutation(N) and the sign for determinant, d

Components

Type Visibility Attributes Name Initial
integer, public :: n
real(kind=long), public, dimension(:, :), pointer :: lu
integer, public, dimension(:), pointer :: permutation
integer, public :: d

Functions

public function lu_decompose(a) result(m)

Carry out LU-decomposition for A and put the decomposed matrix into M implementation of Crout's method with partial pivoting.

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in), dimension(:, :) :: a

n*N matrix

Return Value type(lu_matrix)

public function lu_solve(m, b) result(x)

Solve a linear equation for b with LU-decomposed matrix M Return the solution for b as a new vector x b : row vector, not column vector, that is, just 1D array of n

Arguments

Type IntentOptional Attributes Name
type(lu_matrix), intent(in) :: m
real(kind=long), intent(in), dimension(:) :: b

Return Value real(kind=long), dimension(m%n)

public function lu_inverse(m) result(x)

Return the inverse matrix of LU-decomposed matrix M using F solve(M, b)

Arguments

Type IntentOptional Attributes Name
type(lu_matrix), intent(in) :: m

Return Value real(kind=long), dimension(m%n, m%n)

public function lu_determinant(m) result(d)

returns the determinant value of LU-decomposed matrix M

Arguments

Type IntentOptional Attributes Name
type(lu_matrix), intent(in) :: m

Return Value real(kind=long)

public function matrix_determinant(a) result(det)

returns the determinant of square matrix A slightly redundant to LU_determiant but allows also to return zero as a determinant

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in) :: a(:,:)

Return Value real(kind=long)

public function mattrace(a) result(b)

returns the trace of matrix A

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in) :: a(:,:)

Return Value real(kind=long)

public function mattraceproduct(a, b) result(trace)

computes the trace of the matrix A * B first dimension of A has to match second dimension of B second dimension of A has to match first dimension of B

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in) :: a(:,:)
real(kind=long), intent(in) :: b(:,:)

Return Value real(kind=long)


Subroutines

public subroutine lu_matrix_purge(m)

purges LU_matrix structure

Arguments

Type IntentOptional Attributes Name
type(lu_matrix) :: m

public subroutine print_matrix(a)

Print arbitrary N*M matrix

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in), dimension(:, :) :: a

public subroutine solve_axb(a, b, d, ierr, yes_verbose)

solves AxB=D for D

Arguments

Type IntentOptional Attributes Name
real(kind=long), intent(in) :: a(:,:)

matrix A of shape NxN
real(kind=long), intent(in) :: b(:)

vector B of len N
real(kind=long), intent(out) :: d(size(b,1))

result vector of len N
integer, intent(out) :: ierr

output status -1 is solution could not be found (A is singular)
logical, intent(in), optional :: yes_verbose

print message if A is singular