Example

This example program demonstrates how to use the functions and subroutines contained in the matrices subpackage.

program example_all

  use linalg_kinds, only : wp
  use linalg_matrices_diagonalize, only : diagonalizerealsymmetric
  use linalg_matrices_diagonalize, only : diagonalizehermitian
  use linalg_matrices_diagonalize, only : svd_square
  use linalg_matrices_inverse, only : inversereal
  use linalg_matrices_inverse, only : inversecomplex

  implicit none

  integer, parameter :: dimen = 3
  real(wp) :: A(dimen,dimen)
  complex(wp) :: B(dimen,dimen)
  real(wp) :: C(dimen,dimen)
  real(wp) :: eigenvaluesA(dimen)
  real(wp) :: eigenvectorsA(dimen,dimen)
  real(wp) :: eigenvaluesB(dimen)
  real(wp) :: eigenvaluesC(dimen)
  complex(wp) :: eigenvectorsB(dimen,dimen)
  real(wp) :: zelC(dimen,dimen)
  real(wp) :: zerC(dimen,dimen)
  real(wp) :: Ainv(dimen,dimen)
  complex(wp) :: Binv(dimen,dimen)
  integer :: kont
  integer :: order
  complex(wp) :: test(dimen,dimen)
  real(wp) :: test2(dimen,dimen)
  real(wp) :: test3(dimen,dimen)
  integer :: i

  A(1,1) = 1.0E0_wp
  A(1,2) = 0.0E0_wp
  A(1,3) = 2.0E0_wp
  A(2,1) = 0.0E0_wp
  A(2,2) = 4.0E0_wp
  A(2,3) = 5.0E0_wp
  A(3,1) = 2.0E0_wp
  A(3,2) = 5.0E0_wp
  A(3,3) = 3.0E0_wp

  B = (1,0) * A
  B(1,3) = B(1,3) + (0,1) * A(2,2)
  B(3,1) = B(3,1) - (0,1) * A(2,2)
  B(1,2) = B(1,2) + (0,1) * A(3,3)
  B(2,1) = B(2,1) - (0,1) * A(3,3)
  B(2,3) = B(2,3) + (0,1) * A(1,1)
  B(3,2) = B(3,2) - (0,1) * A(1,1)

  C = 2 * B
  C(3,3) = 1.0E0_wp
  C(1,2) = 0.5E0_wp
  C(3,1) = 2.5E0_wp

  order = 1
  kont = 0

  call diagonalizerealsymmetric(  &
    A,  &
    eigenvaluesA,  &
    eigenvectorsA,  &
    order)

  call diagonalizehermitian(  &
    B,  &
    eigenvaluesB,  &
    eigenvectorsB,  &
    kont)

  call svd_square(  &
    dimen,  &
    C,  &
    eigenvaluesC,  &
    zelC,  &
    zerC,  &
    kont,  &
    order)


  write(*,*)
  write(*,*) "==========================================="
  write(*,*)
  write(*,*) "Diagonalize real symmetric matrix:"
  write(*,*)
  write(*,*) "Matrix:"
  do i=1,dimen
    write(*,*) "(", A(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Eigenvalues:"
  write(*,*) "    ", eigenvaluesA
  write(*,*)
  write(*,*) "Eigenvectors:"
  write(*,*) "    1st: ", eigenvectorsA(1,:)
  write(*,*) "    2nd: ", eigenvectorsA(2,:)
  write(*,*) "    3rd: ", eigenvectorsA(3,:)
  write(*,*)
  write(*,*) "Checks:"
  do i=1,dimen
    write(*,*) "    ", matmul(A,eigenvectorsA(i,:))
    write(*,*) " =? ", eigenvectorsA(i,:) * eigenvaluesA(i)
  enddo
  write(*,*)
  write(*,*) "==========================================="
  write(*,*)
  write(*,*) "Diagonalize hermitian matrix:"
  write(*,*)
  write(*,*) "Matrix:"
  do i=1,dimen
    write(*,*) "(", B(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Eigenvalues:"
  write(*,*) "    ", eigenvaluesB
  write(*,*)
  write(*,*) "Eigenvectors:"
  write(*,*) "    1st: ", eigenvectorsB(1,:)
  write(*,*) "    2nd: ", eigenvectorsB(2,:)
  write(*,*) "    3rd: ", eigenvectorsB(3,:)
  write(*,*)
  write(*,*) "Checks:"
  do i=1,dimen
    write(*,*) "    ", matmul(B, eigenvectorsB(i,:))
    write(*,*) " =? ", eigenvectorsB(i,:) * eigenvaluesB(i)
  enddo
  write(*,*)
  write(*,*) "==========================================="
  write(*,*)
  write(*,*) "Singular value decomposition of real matrix:"
  write(*,*)
  write(*,*) "Matrix:"
  do i=1,dimen
    write(*,*) "(", C(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Eigenvalues:"
  write(*,*) "    ", eigenvaluesC
  write(*,*)
  write(*,*) "Left transformation matrix"
  do i=1,dimen
    write(*,*) "    (", zelC(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Right transformation matrix"
  do i=1,dimen
    write(*,*) "    (", zerC(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Checks:"
  write(*,*) "    Left transformation matrix is orthogonal:"
  test = matmul(zelC, transpose(zelC))
  do i=1,dimen
    write(*,*) "    ", real(test(i,:))
  enddo
  write(*,*) "    Right transformation matrix is orthogonal:"
  test = matmul(zerC, transpose(zerC))
  do i=1,dimen
    write(*,*) "    ", real(test(i,:))
  enddo
  write(*,*) "    Check relation zel^T diag zer = orignal matrix:"
  test2 = 0.0E0_wp
  do i=1,dimen
    test2(i,i) = eigenvaluesC(i)
  enddo
  test3 = matmul(transpose(zelC), matmul(test2, zerC))
  do i=1,dimen
    write(*,*) "    ", test3(i,:)
  enddo
  write(*,*) "    Check relation zel C zer^T = diagonal matrix:"
  test3 = matmul(zelC, matmul(C, transpose(zerC)))
  do i=1,dimen
    write(*,*) "    ", test3(i,:)
  enddo
  write(*,*)
  write(*,*) "==========================================="
  write(*,*)
  write(*,*) "Inverse of real matrix:"
  A(3,2) = 5 * A(1,2)
  Ainv = inversereal(A)
  write(*,*)
  write(*,*) "Matrix:"
  do i=1,dimen
    write(*,*) "(", A(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Inverse matrix:"
  do i=1,dimen
    write(*,*) "(", Ainv(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Check:"
  write(*,*) "    Check that A Ainv = 1:"
  A = matmul(Ainv, A)
  do i=1,dimen
    write(*,*) "    (", A(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "==========================================="
  write(*,*)
  write(*,*) "Inverse of complex matrix:"
  B = (1,0) * A - (0,1) * C
  B(1,3) = 12.0E0_wp
  B(3,2) = -0.05E0_wp
  Binv = inversecomplex(B)
  write(*,*)
  write(*,*) "Matrix:"
  do i=1,dimen
    write(*,*) "(", B(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Inverse matrix:"
  do i=1,dimen
    write(*,*) "(", Binv(i,:), ")"
  enddo
  write(*,*)
  write(*,*) "Check:"
  write(*,*) "    Check that B Binv = 1:"
  B = matmul(Binv, B)
  do i=1,dimen
    write(*,*) "    (", B(i,:), ")"
  enddo
  write(*,*)

end program example_all