Suppose that a matrix 



M


{\displaystyle M}

 contains complex numbers. Then the conjugate transpose of 



M


{\displaystyle M}

 is a matrix 




M

H




{\displaystyle M^{H}}

 containing the complex conjugates of the matrix transposition of 



M


{\displaystyle M}

.

 



(

M

H



)

j
i


=



M

i
j


¯




{\displaystyle (M^{H})_{ji}={\overline {M_{ij}}}}




This means that row 



j


{\displaystyle j}

, column 



i


{\displaystyle i}

 of the conjugate transpose equals the

complex conjugate of row 



i


{\displaystyle i}

, column 



j


{\displaystyle j}

 of the original matrix.



In the next list, 



M


{\displaystyle M}

 must also be a square matrix.

 A Hermitian matrix equals its own conjugate transpose: 




M

H


=
M


{\displaystyle M^{H}=M}

.
 A normal matrix is commutative in multiplication with its conjugate transpose: 




M

H


M
=
M

M

H




{\displaystyle M^{H}M=MM^{H}}

.
 A unitary matrix has its inverse equal to its conjugate transpose: 




M

H


=

M

−
1




{\displaystyle M^{H}=M^{-1}}

. 
 This is true iff 




M

H


M
=

I

n




{\displaystyle M^{H}M=I_{n}}

 and iff 



M

M

H


=

I

n




{\displaystyle MM^{H}=I_{n}}

, where 




I

n




{\displaystyle I_{n}}

 is the identity matrix.


Task

Given some matrix of complex numbers, find its conjugate transpose.

Also determine if the matrix is a:

 Hermitian matrix,
 normal matrix, or 
 unitary matrix.


See also

 MathWorld entry: conjugate transpose
 MathWorld entry: Hermitian matrix
 MathWorld entry: normal matrix
 MathWorld entry: unitary matrix


import complex, strformat
 
type Matrix[M, N: static Positive] = array[M, array[N, Complex[float]]]
 
const Eps = 1e-10   # Tolerance used for float comparisons.
 
 
####################################################################################################
# Templates.
 
template `[]`(m: Matrix; i, j: Natural): Complex[float] =
  ## Allow to get value of an element using m[i, j] syntax.
  m[i][j]
 
template `[]=`(m: var Matrix; i, j: Natural; val: Complex[float]) =
  ## Allow to set value of an element using m[i, j] syntax.
  m[i][j] = val
 
 
####################################################################################################
# General operations.
 
func `$`(m: Matrix): string =
  ## Return the string representation of a matrix using one line per row.
 
  for i, row in m:
    result.add(if i == 0: '[' else: ' ')
    for j, val in row:
      if j != 0: result.add(' ')
      result.add(&"({val.re:7.4f}, {val.im:7.4f})")
    result.add(if i == m.high: ']' else: '\n')
 
#---------------------------------------------------------------------------------------------------
 
func conjugateTransposed[M, N: static int](m: Matrix[M, N]): Matrix[N, M] =
  ## Return the conjugate transpose of a matrix.
 
  for i in 0..<m.M:
    for j in 0..<m.N:
      result[j, i] = m[i, j].conjugate()
 
#---------------------------------------------------------------------------------------------------
 
func `*`[M, K, N: static int](m1: Matrix[M, K]; m2: Matrix[K, N]): Matrix[M, N] =
  # Compute the product of two matrices.
 
  for i in 0..<M:
    for j in 0..<N:
      for k in 0..<K:
        result[i, j] = result[i, j] + m1[i, k] * m2[k, j]
 
 
####################################################################################################
# Properties.
 
func isHermitian(m: Matrix): bool =
  ## Check if a matrix is hermitian.
 
  when m.M != m.N:
    {.error: "hermitian test only allowed for square matrices".}
  else:
    for i in 0..<m.M:
      for j in i..<m.N:
        if m[i, j] != m[j, i].conjugate:
          return false
    result = true
 
#---------------------------------------------------------------------------------------------------
 
func isNormal(m: Matrix): bool =
  ## Check if a matrix is normal.
 
  when m.M != m.N:
    {.error: "normal test only allowed for square matrices".}
  else:
    let h = m.conjugateTransposed
    result = m * h == h * m
 
#---------------------------------------------------------------------------------------------------
 
func isIdentity(m: Matrix): bool =
  ## Check if a matrix is the identity matrix.
 
  when m.M != m.N:
    {.error: "identity test only allowed for square matrices".}
  else:
    for i in 0..<m.M:
      for j in 0..<m.N:
        if i == j:
          if abs(m[i, j] - 1.0) > Eps:
            return false
        else:
          if abs(m[i, j]) > Eps:
            return false
    result = true
 
#---------------------------------------------------------------------------------------------------
 
func isUnitary(m: Matrix): bool =
  ## Check if a matrix is unitary.
 
  when m.M != m.N:
    {.error: "unitary test only allowed for square matrices".}
  else:
    let h = m.conjugateTransposed
    result = (m * h).isIdentity and (h * m).isIdentity
 
#———————————————————————————————————————————————————————————————————————————————————————————————————
 
when isMainModule:
 
  import math
 
  proc test(m: Matrix) =
    echo "\n"
    echo "Matrix"
    echo "------"
    echo m
    echo ""
    echo "Conjugate transposed"
    echo "--------------------"
    echo m.conjugateTransposed
 
    when m.M == m.N:
      # Only for squares matrices.
      echo ""
      echo "Hermitian: ", m.isHermitian
      echo "Normal: ", m.isNormal
      echo "Unitary: ", m.isUnitary
 
  #-------------------------------------------------------------------------------------------------
 
  # Non square matrix.
  const M1: Matrix[2, 3] = [[1.0 + im 2.0, 3.0 + im 0.0, 2.0 + im 5.0],
                            [3.0 - im 1.0, 2.0 + im 0.0, 0.0 + im 3.0]]
 
    # Square matrices.
  const M2: Matrix[2, 2] = [[3.0 + im 0.0, 2.0 + im 1.0],
                            [2.0 - im 1.0, 1.0 + im 0.0]]
 
  const M3: Matrix[3, 3] = [[1.0 + im 0.0, 1.0 + im 0.0, 0.0 + im 0.0],
                            [0.0 + im 0.0, 1.0 + im 0.0, 1.0 + im 0.0],
                            [1.0 + im 0.0, 0.0 + im 0.0, 1.0 + im 0.0]]
 
  const SR2 = 1 / sqrt(2.0)
  const M4: Matrix[3, 3] = [[SR2 + im 0.0, SR2 + im 0.0, 0.0 + im 0.0],
                            [0.0 + im SR2, 0.0 - im SR2, 0.0 + im 0.0],
                            [0.0 + im 0.0, 0.0 + im 0.0, 0.0 + im 1.0]]
 
  test(M1)
  test(M2)
  test(M3)
  test(M4)