Conjugate transpose

From Rosetta Code
Jump to: navigation, search
Task
Conjugate transpose
You are encouraged to solve this task according to the task description, using any language you may know.

Suppose that a matrix M contains complex numbers. Then the conjugate transpose of M is a matrix MH containing the complex conjugates of the matrix transposition of M.

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

This means that row j, column i of the conjugate transpose equals the complex conjugate of row i, column j of the original matrix.

In the next list, M must also be a square matrix.

Given some matrix of complex numbers, find its conjugate transpose. Also determine if it is a Hermitian matrix, normal matrix, or a unitary matrix.

Contents

[edit] Ada

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
with Ada.Numerics.Complex_Arrays; use Ada.Numerics.Complex_Arrays;
procedure ConTrans is
subtype CM is Complex_Matrix;
S2O2 : constant Float := 0.7071067811865;
 
procedure Print (mat : CM) is begin
for row in mat'Range(1) loop for col in mat'Range(2) loop
Put(mat(row,col), Exp=>0, Aft=>4);
end loop; New_Line; end loop;
end Print;
 
function almostzero(mat : CM; tol : Float) return Boolean is begin
for row in mat'Range(1) loop for col in mat'Range(2) loop
if abs(mat(row,col)) > tol then return False; end if;
end loop; end loop;
return True;
end almostzero;
 
procedure Examine (mat : CM) is
CT : CM := Conjugate (Transpose(mat));
isherm, isnorm, isunit : Boolean;
begin
isherm := almostzero(mat-CT, 1.0e-6);
isnorm := almostzero(mat*CT-CT*mat, 1.0e-6);
isunit := almostzero(CT-Inverse(mat), 1.0e-6);
Print(mat);
Put_Line("Conjugate transpose:"); Print(CT);
Put_Line("Hermitian?: " & isherm'Img);
Put_Line("Normal?: " & isnorm'Img);
Put_Line("Unitary?: " & isunit'Img);
end Examine;
 
hmat : CM := ((3.0+0.0*i, 2.0+1.0*i), (2.0-1.0*i, 1.0+0.0*i));
nmat : CM := ((1.0+0.0*i, 1.0+0.0*i, 0.0+0.0*i),
(0.0+0.0*i, 1.0+0.0*i, 1.0+0.0*i),
(1.0+0.0*i, 0.0+0.0*i, 1.0+0.0*i));
umat : CM := ((S2O2+0.0*i, S2O2+0.0*i, 0.0+0.0*i),
(0.0+S2O2*i, 0.0-S2O2*i, 0.0+0.0*i),
(0.0+0.0*i, 0.0+0.0*i, 0.0+1.0*i));
begin
Put_Line("hmat:"); Examine(hmat); New_Line;
Put_Line("nmat:"); Examine(nmat); New_Line;
Put_Line("umat:"); Examine(umat);
end ConTrans;
Output:
hmat:
( 3.0000, 0.0000)( 2.0000, 1.0000)
( 2.0000,-1.0000)( 1.0000, 0.0000)
Conjugate transpose:
( 3.0000,-0.0000)( 2.0000, 1.0000)
( 2.0000,-1.0000)( 1.0000,-0.0000)
Hermitian?: TRUE
Normal?: TRUE
Unitary?: FALSE

nmat:
( 1.0000, 0.0000)( 1.0000, 0.0000)( 0.0000, 0.0000)
( 0.0000, 0.0000)( 1.0000, 0.0000)( 1.0000, 0.0000)
( 1.0000, 0.0000)( 0.0000, 0.0000)( 1.0000, 0.0000)
Conjugate transpose:
( 1.0000,-0.0000)( 0.0000,-0.0000)( 1.0000,-0.0000)
( 1.0000,-0.0000)( 1.0000,-0.0000)( 0.0000,-0.0000)
( 0.0000,-0.0000)( 1.0000,-0.0000)( 1.0000,-0.0000)
Hermitian?: FALSE
Normal?: TRUE
Unitary?: FALSE

umat:
( 0.7071, 0.0000)( 0.7071, 0.0000)( 0.0000, 0.0000)
( 0.0000, 0.7071)( 0.0000,-0.7071)( 0.0000, 0.0000)
( 0.0000, 0.0000)( 0.0000, 0.0000)( 0.0000, 1.0000)
Conjugate transpose:
( 0.7071,-0.0000)( 0.0000,-0.7071)( 0.0000,-0.0000)
( 0.7071,-0.0000)( 0.0000, 0.7071)( 0.0000,-0.0000)
( 0.0000,-0.0000)( 0.0000,-0.0000)( 0.0000,-1.0000)
Hermitian?: FALSE
Normal?: TRUE
Unitary?: TRUE

[edit] C

/*28th August, 2012
Abhishek Ghosh
 
Uses C99 specified complex.h, complex datatype has to be defined and operation provided if used on non-C99 compilers */

 
#include<stdlib.h>
#include<stdio.h>
#include<complex.h>
 
typedef struct
{
int rows, cols;
complex **z;
} matrix;
 
matrix
transpose (matrix a)
{
int i, j;
matrix b;
 
b.rows = a.cols;
b.cols = a.rows;
 
b.z = malloc (b.rows * sizeof (complex *));
 
for (i = 0; i < b.rows; i++)
{
b.z[i] = malloc (b.cols * sizeof (complex));
for (j = 0; j < b.cols; j++)
{
b.z[i][j] = conj (a.z[j][i]);
}
}
 
return b;
}
 
int
isHermitian (matrix a)
{
int i, j;
matrix b = transpose (a);
 
if (b.rows == a.rows && b.cols == a.cols)
{
for (i = 0; i < b.rows; i++)
{
for (j = 0; j < b.cols; j++)
{
if (b.z[i][j] != a.z[i][j])
return 0;
}
}
}
 
else
return 0;
 
return 1;
}
 
matrix
multiply (matrix a, matrix b)
{
matrix c;
int i, j;
 
if (a.cols == b.rows)
{
c.rows = a.rows;
c.cols = b.cols;
 
c.z = malloc (c.rows * (sizeof (complex *)));
 
for (i = 0; i < c.rows; i++)
{
c.z[i] = malloc (c.cols * sizeof (complex));
c.z[i][j] = 0 + 0 * I;
for (j = 0; j < b.cols; j++)
{
c.z[i][j] += a.z[i][j] * b.z[j][i];
}
}
 
}
 
return c;
}
 
int
isNormal (matrix a)
{
int i, j;
matrix a_ah, ah_a;
 
if (a.rows != a.cols)
return 0;
 
a_ah = multiply (a, transpose (a));
ah_a = multiply (transpose (a), a);
 
for (i = 0; i < a.rows; i++)
{
for (j = 0; j < a.cols; j++)
{
if (a_ah.z[i][j] != ah_a.z[i][j])
return 0;
}
}
 
return 1;
}
 
int
isUnitary (matrix a)
{
matrix b;
int i, j;
if (isNormal (a) == 1)
{
b = multiply (a, transpose(a));
 
for (i = 0; i < b.rows; i++)
{
for (j = 0; j < b.cols; j++)
{
if ((i == j && b.z[i][j] != 1) || (i != j && b.z[i][j] != 0))
return 0;
}
}
return 1;
}
return 0;
}
 
 
int
main ()
{
complex z = 3 + 4 * I;
matrix a, aT;
int i, j;
printf ("Enter rows and columns :");
scanf ("%d%d", &a.rows, &a.cols);
 
a.z = malloc (a.rows * sizeof (complex *));
printf ("Randomly Generated Complex Matrix A is : ");
for (i = 0; i < a.rows; i++)
{
printf ("\n");
a.z[i] = malloc (a.cols * sizeof (complex));
for (j = 0; j < a.cols; j++)
{
a.z[i][j] = rand () % 10 + rand () % 10 * I;
printf ("\t%f + %fi", creal (a.z[i][j]), cimag (a.z[i][j]));
}
}
 
aT = transpose (a);
 
printf ("\n\nTranspose of Complex Matrix A is : ");
for (i = 0; i < aT.rows; i++)
{
printf ("\n");
aT.z[i] = malloc (aT.cols * sizeof (complex));
for (j = 0; j < aT.cols; j++)
{
aT.z[i][j] = rand () % 10 + rand () % 10 * I;
printf ("\t%f + %fi", creal (aT.z[i][j]), cimag (aT.z[i][j]));
}
}
 
printf ("\n\nComplex Matrix A %s hermitian",
isHermitian (a) == 1 ? "is" : "is not");
printf ("\n\nComplex Matrix A %s unitary",
isUnitary (a) == 1 ? "is" : "is not");
printf ("\n\nComplex Matrix A %s normal",
isNormal (a) == 1 ? "is" : "is not");
 
 
 
return 0;
}
Output:
Enter rows and columns :3 3
Randomly Generated Complex Matrix A is :
        3.000000 + 6.000000i    7.000000 + 5.000000i    3.000000 + 5.000000i
        6.000000 + 2.000000i    9.000000 + 1.000000i    2.000000 + 7.000000i
        0.000000 + 9.000000i    3.000000 + 6.000000i    0.000000 + 6.000000i

Transpose of Complex Matrix A is :
        2.000000 + 6.000000i    1.000000 + 8.000000i    7.000000 + 9.000000i
        2.000000 + 0.000000i    2.000000 + 3.000000i    7.000000 + 5.000000i
        9.000000 + 2.000000i    2.000000 + 8.000000i    9.000000 + 7.000000i

Complex Matrix A is not hermitian

Complex Matrix A is not unitary

Complex Matrix A is not normal

[edit] D

Translation of: Python
A well typed and mostly imperative version:
import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric;
 
T[][] conjugateTranspose(T)(in T[][] m) pure nothrow {
auto r = new typeof(return)(m[0].length, m.length);
foreach (immutable nr, const row; m)
foreach (immutable nc, immutable c; row)
r[nc][nr] = c.conj;
return r;
}
 
bool isRectangular(T)(in T[][] M) pure nothrow @nogc {
return M.all!(row => row.length == M[0].length);
}
 
T[][] matMul(T)(in T[][] A, in T[][] B) pure nothrow
in {
assert(A.isRectangular && B.isRectangular &&
!A.empty && !B.empty && A[0].length == B.length);
} body {
auto result = new T[][](A.length, B[0].length);
auto aux = new T[B.length];
 
foreach (immutable j; 0 .. B[0].length) {
foreach (immutable k, const row; B)
aux[k] = row[j];
foreach (immutable i, const ai; A)
result[i][j] = dotProduct(ai, aux);
}
 
return result;
}
 
/// Check any number of complex matrices for equality within
/// some bits of mantissa.
bool areEqual(T)(in Complex!T[][][] matrices, in size_t nBits=20)
pure nothrow {
static bool allSame(U)(in U[] v) pure nothrow @nogc {
return v[1 .. $].all!(c => c == v[0]);
}
 
bool allNearSame(in Complex!T[] v) pure nothrow @nogc {
auto v0 = v[0].Complex!T; // To avoid another cast.
return v[1 .. $].all!(c=> feqrel(v0.re, cast()c.re) >= nBits &&
feqrel(v0.im, cast()c.im) >= nBits);
}
 
immutable x = matrices.map!(m => m.length).array;
if (!allSame(x))
return false;
immutable y = matrices.map!(m => m[0].length).array;
if (!allSame(y))
return false;
foreach (immutable s; 0 .. x[0])
foreach (immutable t; 0 .. y[0])
if (!allNearSame(matrices.map!(m => m[s][t]).array))
return false;
return true;
}
 
bool isHermitian(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow {
return [m, ct].areEqual;
}
 
bool isNormal(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow {
return [matMul(m, ct), matMul(ct, m)].areEqual;
}
 
auto complexIdentitymatrix(in size_t side) pure nothrow {
return side.iota
.map!((in r) => side.iota.map!(c => complex(r == c)).array)
.array;
}
 
bool isUnitary(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow {
immutable mct = matMul(m, ct);
immutable ident = mct.length.complexIdentitymatrix;
return [mct, matMul(ct, m), ident].areEqual;
}
 
void main() {
alias C = complex;
immutable x = 2 ^^ 0.5 / 2;
 
immutable data = [[[C(3.0, 0.0), C(2.0, 1.0)],
[C(2.0, -1.0), C(1.0, 0.0)]],
 
[[C(1.0, 0.0), C(1.0, 0.0), C(0.0, 0.0)],
[C(0.0, 0.0), C(1.0, 0.0), C(1.0, 0.0)],
[C(1.0, 0.0), C(0.0, 0.0), C(1.0, 0.0)]],
 
[[C(x, 0.0), C(x, 0.0), C(0.0, 0.0)],
[C(0.0, -x), C(0.0, x), C(0.0, 0.0)],
[C(0.0, 0.0), C(0.0, 0.0), C(0.0, 1.0)]]];
 
foreach (immutable mat; data) {
enum mFormat = "[%([%(%1.3f, %)],\n %)]]";
writefln("Matrix:\n" ~ mFormat, mat);
immutable ct = conjugateTranspose(mat);
"Its conjugate transpose:".writeln;
writefln(mFormat, ct);
writefln("Hermitian? %s.", isHermitian(mat, ct));
writefln("Normal?  %s.", isNormal(mat, ct));
writefln("Unitary?  %s.\n", isUnitary(mat, ct));
}
}
Output:
Matrix:
[[3.000+0.000i, 2.000+1.000i],
 [2.000-1.000i, 1.000+0.000i]]
Its conjugate transpose:
[[3.000-0.000i, 2.000+1.000i],
 [2.000-1.000i, 1.000-0.000i]]
Hermitian? true.
Normal?    true.
Unitary?   false.

Matrix:
[[1.000+0.000i, 1.000+0.000i, 0.000+0.000i],
 [0.000+0.000i, 1.000+0.000i, 1.000+0.000i],
 [1.000+0.000i, 0.000+0.000i, 1.000+0.000i]]
Its conjugate transpose:
[[1.000-0.000i, 0.000-0.000i, 1.000-0.000i],
 [1.000-0.000i, 1.000-0.000i, 0.000-0.000i],
 [0.000-0.000i, 1.000-0.000i, 1.000-0.000i]]
Hermitian? false.
Normal?    true.
Unitary?   false.

Matrix:
[[0.707+0.000i, 0.707+0.000i, 0.000+0.000i],
 [0.000-0.707i, 0.000+0.707i, 0.000+0.000i],
 [0.000+0.000i, 0.000+0.000i, 0.000+1.000i]]
Its conjugate transpose:
[[0.707-0.000i, 0.000+0.707i, 0.000-0.000i],
 [0.707-0.000i, 0.000-0.707i, 0.000-0.000i],
 [0.000-0.000i, 0.000-0.000i, 0.000-1.000i]]
Hermitian? false.
Normal?    true.
Unitary?   true.

[edit] Alternative Version

A more functional version that contains some typing problems (same output).

import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric, std.exception, std.traits;
 
// alias CM(T) = Complex!T[][]; // Not yet useful.
 
auto conjugateTranspose(T)(in Complex!T[][] m) /*pure*/ nothrow
if (!hasIndirections!T) {
return iota(m[0].length)
.map!(i => m.transversal(i).map!conj.array)
.array
.assumeUnique; //*
}
 
T[][] matMul(T)(immutable T[][] A, immutable T[][] B) pure nothrow {
immutable Bt = B[0].length.iota.map!(i=> B.transversal(i).array)
.array;
return A.map!((in a) => Bt.map!(b => a.dotProduct(b)).array).array;
}
 
/// Check any number of complex matrices for equality within
/// some bits of mantissa.
bool areEqual(T)(in Complex!T[][][] matrices, in size_t nBits=20)
pure nothrow {
static bool allSame(U)(in U[] v) pure nothrow @nogc {
return v[1 .. $].all!(c => c == v[0]);
}
 
bool allNearSame(in Complex!T[] v) pure nothrow @nogc {
auto v0 = v[0].Complex!T; // To avoid another cast.
return v[1 .. $].all!(c=> feqrel(v0.re, cast()c.re) >= nBits &&
feqrel(v0.im, cast()c.im) >= nBits);
}
 
immutable x = matrices.map!(m => m.length).array;
if (!allSame(x))
return false;
immutable y = matrices.map!(m => m[0].length).array;
if (!allSame(y))
return false;
foreach (immutable s; 0 .. x[0])
foreach (immutable t; 0 .. y[0])
if (!allNearSame(matrices.map!(m => m[s][t]).array))
return false;
return true;
}
 
bool isHermitian(T)(in Complex!T[][] m, in Complex!T[][] ct)
pure nothrow {
return [m, ct].areEqual;
}
 
bool isNormal(T)(immutable Complex!T[][] m, immutable Complex!T[][] ct)
pure nothrow {
return [matMul(m, ct), matMul(ct, m)].areEqual;
}
 
auto complexIdentitymatrix(in size_t side) pure nothrow {
return side.iota
.map!((in r) => side.iota.map!(c => complex(r == c)).array)
.array;
}
 
bool isUnitary(T)(immutable Complex!T[][] m, immutable Complex!T[][] ct)
pure nothrow {
immutable mct = matMul(m, ct);
immutable ident = mct.length.complexIdentitymatrix;
return [mct, matMul(ct, m), ident].areEqual;
}
 
void main() {
alias C = complex;
immutable x = 2 ^^ 0.5 / 2;
 
foreach (/*immutable*/ const matrix;
[[[C(3.0, 0.0), C(2.0, 1.0)],
[C(2.0, -1.0), C(1.0, 0.0)]],
 
[[C(1.0, 0.0), C(1.0, 0.0), C(0.0, 0.0)],
[C(0.0, 0.0), C(1.0, 0.0), C(1.0, 0.0)],
[C(1.0, 0.0), C(0.0, 0.0), C(1.0, 0.0)]],
 
[[C(x, 0.0), C(x, 0.0), C(0.0, 0.0)],
[C(0.0, -x), C(0.0, x), C(0.0, 0.0)],
[C(0.0, 0.0), C(0.0, 0.0), C(0.0, 1.0)]]]) {
immutable mat = matrix.assumeUnique; //*
 
enum mFormat = "[%([%(%1.3f, %)],\n %)]]";
writefln("Matrix:\n" ~ mFormat, mat);
immutable ct = conjugateTranspose(mat);
"Its conjugate transpose:".writeln;
writefln(mFormat, ct);
writefln("Hermitian? %s.", isHermitian(mat, ct));
writefln("Normal?  %s.", isNormal(mat, ct));
writefln("Unitary?  %s.\n", isUnitary(mat, ct));
}
}

[edit] Factor

Before the fix to Factor bug #484, m. gave the wrong answer and this code failed. Factor 0.94 is too old to work.

Works with: Factor version development (future 0.95)
USING: kernel math.functions math.matrices sequences ;
IN: rosetta.hermitian
 
: conj-t ( matrix -- conjugate-transpose )
flip [ [ conjugate ] map ] map ;
 
: hermitian-matrix? ( matrix -- ? )
dup conj-t = ;
 
: normal-matrix? ( matrix -- ? )
dup conj-t [ m. ] [ swap m. ] 2bi = ;
 
: unitary-matrix? ( matrix -- ? )
[ dup conj-t m. ] [ length identity-matrix ] bi = ;

Usage:

USE: rosetta.hermitian
IN: scratchpad { { C{ 1 2 } 0 }
                 { 0 C{ 3 4 } } }
               [ hermitian-matrix? . ]
               [ normal-matrix? . ]
               [ unitary-matrix? . ] tri
f
t
f

[edit] Fortran

The examples and algorithms are taken from the j solution, except for UnitaryQ. The j solution uses the matrix inverse verb. Compilation on linux, assuming the program is file f.f08 :
gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none f.f08 -o f
 
program conjugate_transpose
 
complex, dimension(3, 3) :: a
integer :: i
a = reshape((/ (i, i=1,9) /), shape(a))
call characterize(a)
a(:2,:2) = reshape((/cmplx(3,0),cmplx(2,-1),cmplx(2,1),cmplx(1,0)/),(/2,2/))
call characterize(a(:2,:2))
call characterize(cmplx(reshape((/1,0,1,1,1,0,0,1,1/),(/3,3/)),0))
a(3,:) = (/cmplx(0,0), cmplx(0,0), cmplx(0,1)/)*sqrt(2.0)
a(2,:) = (/cmplx(0,-1),cmplx(0,1),cmplx(0,0)/)
a(1,:) = (/1,1,0/)
a = a * sqrt(2.0)/2.0
call characterize(a)
 
contains
 
subroutine characterize(a)
complex, dimension(:,:), intent(in) :: a
integer :: i, j
do i=1, size(a,1)
print *,(a(i, j), j=1,size(a,1))
end do
print *,'Is Hermitian? ',HermitianQ(a)
print *,'Is normal? ',NormalQ(a)
print *,'Unitary? ',UnitaryQ(a)
print '(/)'
end subroutine characterize
 
function ct(a) result(b) ! return the conjugate transpose of a matrix
complex, dimension(:,:), intent(in) :: a
complex, dimension(size(a,1),size(a,1)) :: b
b = conjg(transpose(a))
end function ct
 
function identity(n) result(b) ! return identity matrix
integer, intent(in) :: n
real, dimension(n,n) :: b
integer :: i
b = 0
do i=1, n
b(i,i) = 1
end do
end function identity
 
logical function HermitianQ(a)
complex, dimension(:,:), intent(in) :: a
HermitianQ = all(a .eq. ct(a))
end function HermitianQ
 
logical function NormalQ(a)
complex, dimension(:,:), intent(in) :: a
NormalQ = all(matmul(ct(a),a) .eq. matmul(a,ct(a)))
end function NormalQ
 
logical function UnitaryQ(a)
! if A inverse equals A star
! then multiplying each side by A should result in the identity matrix
! Thus show that A times A star is sufficiently close to I .
complex, dimension(:,:), intent(in) :: a
UnitaryQ = all(abs(matmul(a,ct(a)) - identity(size(a,1))) .lt. 1e-6)
end function UnitaryQ
 
end program conjugate_transpose
 
-*- mode: compilation; default-directory: "/tmp/" -*-
Compilation started at Fri Jun  7 16:31:38

a=./f && make $a && time $a
gfortran -std=f2008 -Wall -fopenmp -ffree-form -fall-intrinsics -fimplicit-none f.f08 -o f
 (  1.00000000    ,  0.00000000    ) (  4.00000000    ,  0.00000000    ) (  7.00000000    ,  0.00000000    )
 (  2.00000000    ,  0.00000000    ) (  5.00000000    ,  0.00000000    ) (  8.00000000    ,  0.00000000    )
 (  3.00000000    ,  0.00000000    ) (  6.00000000    ,  0.00000000    ) (  9.00000000    ,  0.00000000    )
 Is Hermitian?   F
 Is normal?   F
 Unitary?   F


 (  3.00000000    ,  0.00000000    ) (  2.00000000    ,  1.00000000    )
 (  2.00000000    , -1.00000000    ) (  1.00000000    ,  0.00000000    )
 Is Hermitian?   T
 Is normal?   T
 Unitary?   F


 (  1.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    )
 (  1.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    ) (  1.00000000    ,  0.00000000    )
 Is Hermitian?   F
 Is normal?   T
 Unitary?   F


 ( 0.707106769    ,  0.00000000    ) ( 0.707106769    ,  0.00000000    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,-0.707106769    ) (  0.00000000    , 0.707106769    ) (  0.00000000    ,  0.00000000    )
 (  0.00000000    ,  0.00000000    ) (  0.00000000    ,  0.00000000    ) (  0.00000000    , 0.999999940    )
 Is Hermitian?   F
 Is normal?   T
 Unitary?   T



real	0m0.002s
user	0m0.000s
sys	0m0.000s

Compilation finished at Fri Jun  7 16:31:38

[edit] Go

package main
 
import (
"fmt"
"math"
"math/cmplx"
)
 
// a type to represent matrices
type matrix struct {
ele []complex128
cols int
}
 
// conjugate transpose, implemented here as a method on the matrix type.
func (m *matrix) conjTranspose() *matrix {
r := &matrix{make([]complex128, len(m.ele)), len(m.ele) / m.cols}
rx := 0
for _, e := range m.ele {
r.ele[rx] = cmplx.Conj(e)
rx += r.cols
if rx >= len(r.ele) {
rx -= len(r.ele) - 1
}
}
return r
}
 
// program to demonstrate capabilites on example matricies
func main() {
show("h", matrixFromRows([][]complex128{
{3, 2 + 1i},
{2 - 1i, 1}}))
 
show("n", matrixFromRows([][]complex128{
{1, 1, 0},
{0, 1, 1},
{1, 0, 1}}))
 
show("u", matrixFromRows([][]complex128{
{math.Sqrt2 / 2, math.Sqrt2 / 2, 0},
{math.Sqrt2 / -2i, math.Sqrt2 / 2i, 0},
{0, 0, 1i}}))
}
 
func show(name string, m *matrix) {
m.print(name)
ct := m.conjTranspose()
ct.print(name + "_ct")
 
fmt.Println("Hermitian:", m.equal(ct, 1e-14))
 
mct := m.mult(ct)
ctm := ct.mult(m)
fmt.Println("Normal:", mct.equal(ctm, 1e-14))
 
i := eye(m.cols)
fmt.Println("Unitary:", mct.equal(i, 1e-14) && ctm.equal(i, 1e-14))
}
 
// two constructors
func matrixFromRows(rows [][]complex128) *matrix {
m := &matrix{make([]complex128, len(rows)*len(rows[0])), len(rows[0])}
for rx, row := range rows {
copy(m.ele[rx*m.cols:(rx+1)*m.cols], row)
}
return m
}
 
func eye(n int) *matrix {
r := &matrix{make([]complex128, n*n), n}
n++
for x := 0; x < len(r.ele); x += n {
r.ele[x] = 1
}
return r
}
 
// print method outputs matrix to stdout
func (m *matrix) print(heading string) {
fmt.Print("\n", heading, "\n")
for e := 0; e < len(m.ele); e += m.cols {
fmt.Printf("%6.3f ", m.ele[e:e+m.cols])
fmt.Println()
}
}
 
// equal method uses ε to allow for floating point error.
func (a *matrix) equal(b *matrix, ε float64) bool {
for x, aEle := range a.ele {
if math.Abs(real(aEle)-real(b.ele[x])) > math.Abs(real(aEle))*ε ||
math.Abs(imag(aEle)-imag(b.ele[x])) > math.Abs(imag(aEle))*ε {
return false
}
}
return true
}
 
// mult method taken from matrix multiply task
func (m1 *matrix) mult(m2 *matrix) (m3 *matrix) {
m3 = &matrix{make([]complex128, (len(m1.ele)/m1.cols)*m2.cols), m2.cols}
for m1c0, m3x := 0, 0; m1c0 < len(m1.ele); m1c0 += m1.cols {
for m2r0 := 0; m2r0 < m2.cols; m2r0++ {
for m1x, m2x := m1c0, m2r0; m2x < len(m2.ele); m2x += m2.cols {
m3.ele[m3x] += m1.ele[m1x] * m2.ele[m2x]
m1x++
}
m3x++
}
}
return m3
}

Output:

h
[( 3.000+0.000i) (+2.000+1.000i)] 
[( 2.000-1.000i) (+1.000+0.000i)] 

h_ct
[( 3.000-0.000i) (+2.000+1.000i)] 
[( 2.000-1.000i) (+1.000-0.000i)] 
Hermitian: true
Normal: true
Unitary: false

n
[( 1.000+0.000i) (+1.000+0.000i) (+0.000+0.000i)] 
[( 0.000+0.000i) (+1.000+0.000i) (+1.000+0.000i)] 
[( 1.000+0.000i) (+0.000+0.000i) (+1.000+0.000i)] 

n_ct
[( 1.000-0.000i) (+0.000-0.000i) (+1.000-0.000i)] 
[( 1.000-0.000i) (+1.000-0.000i) (+0.000-0.000i)] 
[( 0.000-0.000i) (+1.000-0.000i) (+1.000-0.000i)] 
Hermitian: false
Normal: true
Unitary: false

u
[( 0.707+0.000i) (+0.707+0.000i) (+0.000+0.000i)] 
[( 0.000+0.707i) (+0.000-0.707i) (+0.000+0.000i)] 
[( 0.000+0.000i) (+0.000+0.000i) (+0.000+1.000i)] 

u_ct
[( 0.707-0.000i) (+0.000-0.707i) (+0.000-0.000i)] 
[( 0.707-0.000i) (+0.000+0.707i) (+0.000-0.000i)] 
[( 0.000-0.000i) (+0.000-0.000i) (+0.000-1.000i)] 
Hermitian: false
Normal: true
Unitary: true

[edit] J

Solution:
   ct =: +@|:                      NB.  Conjugate transpose (ct A is A_ct)
Examples:
   X         =: +/ . *             NB. Matrix Multiply (x)
 
HERMITIAN =: 3 2j1 ,: 2j_1 1
(-: ct) HERMITIAN NB. A_ct = A
1
 
NORMAL =: 1 1 0 , 0 1 1 ,: 1 0 1
((X~ -: X) ct) NORMAL NB. A_ct x A = A x A_ct
1
 
UNITARY =: (-:%:2) * 1 1 0 , 0j_1 0j1 0 ,: 0 0 0j1 * %:2
(ct -: %.) UNITARY NB. A_ct = A^-1
1
Reference (example matrices for other langs to use):
   HERMITIAN;NORMAL;UNITARY
+--------+-----+--------------------------+
| 3 2j1|1 1 0| 0.707107 0.707107 0|
|2j_1 1|0 1 1|0j_0.707107 0j0.707107 0|
| |1 0 1| 0 0 0j1|
+--------+-----+--------------------------+
NB. In J, PjQ is P + Q*i and the 0.7071... is sqrt(2)
 
hermitian=: -: ct
normal =: (X~ -: X) ct
unitary=: ct -: %.
 
(hermitian,normal,unitary)&.>HERMITIAN;NORMAL;UNITARY
+-----+-----+-----+
|1 1 0|0 1 0|0 1 1|
+-----+-----+-----+

[edit] Julia

Julia has a built-in matrix type, and the conjugate-transpose of a complex matrix A is simply:

A'

(similar to Matlab). You can check whether A is Hermitian via the built-in function

ishermitian(A)

Ignoring the possibility of roundoff errors for floating-point matrices (like most of the examples in the other languages), you can check whether a matrix is normal or unitary by the following functions

isnormal(A) = size(A,1) == size(A,2) && A'*A == A*A'
isunitary(A) = size(A,1) == size(A,2) && A'*A == eye(A)

[edit] Maple

The commands HermitianTranspose and IsUnitary are provided by the LinearAlgebra package.

M:=<<3|2+I>,<2-I|1>>:
 
with(LinearAlgebra):
IsNormal:=A->EqualEntries(A^%H.A,A.A^%H):
 
M^%H;
HermitianTranspose(M);
type(M,'Matrix'(hermitian));
IsNormal(M);
IsUnitary(M);

Output:

                               [  3    2 + I]
                               [            ]
                               [2 - I    1  ]

                               [  3    2 + I]
                               [            ]
                               [2 - I    1  ]

                                    true

                                    true

                                    false

[edit] Mathematica

NormalMatrixQ[a_List?MatrixQ] := Module[{b = Conjugate@Transpose@a},a.b === b.a]
UnitaryQ[m_List?MatrixQ] := (Conjugate@Transpose@m.m == IdentityMatrix@Length@m)
 
m = {{1, 2I, 3}, {3+4I, 5, I}};
m //MatrixForm
->
(1 2I 3
3+4I 5 I)
 
ConjugateTranspose[m] //MatrixForm
->
(1 3-4I
-2I 5
3 -I)
 
{HermitianMatrixQ@#, NormalMatrixQ@#, UnitaryQ@#}&@m
-> {False, False, False}

[edit] PARI/GP

conjtranspose(M)=conj(M~)
isHermitian(M)=M==conj(M~)
isnormal(M)=my(H=conj(M~));H*M==M*H
isunitary(M)=M*conj(M~)==1

[edit] PL/I

 
test: procedure options (main); /* 1 October 2012 */
declare n fixed binary;
 
put ('Conjugate a complex square matrix.');
put skip list ('What is the order of the matrix?:');
get (n);
begin;
declare (M, MH, MM, MM_MMH, MM_MHM, IDENTITY)(n,n) fixed complex;
declare i fixed binary;
 
IDENTITY = 0; do i = 1 to n; IDENTITY(I,I) = 1; end;
put skip list ('Please type the matrix:');
get list (M);
do i = 1 to n;
put skip list (M(i,*));
end;
do i = 1 to n;
MH(i,*) = conjg(M(*,i));
end;
put skip list ('The conjugate transpose is:');
do i = 1 to n;
put skip list (MH(i,*));
end;
if all(M=MH) then
put skip list ('Matrix is Hermitian');
call MMULT(M, MH, MM_MMH);
call MMULT(MH, M, MM_MHM);
 
if all(MM_MMH = MM_MHM) then
put skip list ('Matrix is Normal');
 
if all(ABS(MM_MMH - IDENTITY) < 0.0001) then
put skip list ('Matrix is unitary');
if all(ABS(MM_MHM - IDENTITY) < 0.0001) then
put skip list ('Matrix is unitary');
end;
 
MMULT: procedure (M, MH, MM);
declare (M, MH, MM)(*,*) fixed complex;
declare (i, j, n) fixed binary;
 
n = hbound(M,1);
do i = 1 to n;
do j = 1 to n;
MM(i,j) = sum(M(i,*) * MH(*,j) );
end;
end;
end MMULT;
end test;
 

Outputs from separate runs:

Please type the matrix: 

       1+0I                    1+0I                    1+0I       
       1+0I                    1+0I                    1+0I       
       1+0I                    1+0I                    1+0I       
The conjugate transpose is: 
       1-0I                    1-0I                    1-0I       
       1-0I                    1-0I                    1-0I       
       1-0I                    1-0I                    1-0I       
Matrix is Hermitian 
Matrix is Normal 

       1+0I                    1+0I                    0+0I
       0+0I                    1+0I                    1+0I       
       1+0I                    0+0I                    1+0I       
The conjugate transpose is: 
       1-0I                    0-0I                    1-0I       
       1-0I                    1-0I                    0-0I       
       0-0I                    1-0I                    1-0I       
Matrix is Normal 

Next test performed with declaration of matrixes changed to decimal precision (10,5).

Please type the matrix:

      0.70710+0.00000I        0.70710+0.00000I        0.00000+0.00000I
      0.00000+0.70710I        0.00000-0.70710I        0.00000+0.00000I
      0.00000+0.00000I        0.00000+0.00000I        0.00000+1.00000I
    
The conjugate transpose is: 
      0.70710-0.00000I        0.00000-0.70710I        0.00000-0.00000I
      0.70710-0.00000I        0.00000+0.70710I        0.00000-0.00000I
      0.00000-0.00000I        0.00000-0.00000I        0.00000-1.00000I

Matrix is Normal 
Matrix is unitary 
Matrix is unitary

[edit] Python

Internally, matrices must be represented as rectangular tuples of tuples of complex numbers.

def conjugate_transpose(m):
return tuple(tuple(n.conjugate() for n in row) for row in zip(*m))
 
def mmul( ma, mb):
return tuple(tuple(sum( ea*eb for ea,eb in zip(a,b)) for b in zip(*mb)) for a in ma)
 
def mi(size):
'Complex Identity matrix'
sz = range(size)
m = [[0 + 0j for i in sz] for j in sz]
for i in range(size):
m[i][i] = 1 + 0j
return tuple(tuple(row) for row in m)
 
def __allsame(vector):
first, rest = vector[0], vector[1:]
return all(i == first for i in rest)
 
def __allnearsame(vector, eps=1e-14):
first, rest = vector[0], vector[1:]
return all(abs(first.real - i.real) < eps and abs(first.imag - i.imag) < eps
for i in rest)
 
def isequal(matrices, eps=1e-14):
'Check any number of matrices for equality within eps'
x = [len(m) for m in matrices]
if not __allsame(x): return False
y = [len(m[0]) for m in matrices]
if not __allsame(y): return False
for s in range(x[0]):
for t in range(y[0]):
if not __allnearsame([m[s][t] for m in matrices], eps): return False
return True
 
 
def ishermitian(m, ct):
return isequal([m, ct])
 
def isnormal(m, ct):
return isequal([mmul(m, ct), mmul(ct, m)])
 
def isunitary(m, ct):
mct, ctm = mmul(m, ct), mmul(ct, m)
mctx, mcty, cmx, ctmy = len(mct), len(mct[0]), len(ctm), len(ctm[0])
ident = mi(mctx)
return isequal([mct, ctm, ident])
 
def printm(comment, m):
print(comment)
fields = [['%g%+gj' % (f.real, f.imag) for f in row] for row in m]
width = max(max(len(f) for f in row) for row in fields)
lines = (', '.join('%*s' % (width, f) for f in row) for row in fields)
print('\n'.join(lines))
 
if __name__ == '__main__':
for matrix in [
((( 3.000+0.000j), (+2.000+1.000j)),
(( 2.000-1.000j), (+1.000+0.000j))),
 
((( 1.000+0.000j), (+1.000+0.000j), (+0.000+0.000j)),
(( 0.000+0.000j), (+1.000+0.000j), (+1.000+0.000j)),
(( 1.000+0.000j), (+0.000+0.000j), (+1.000+0.000j))),
 
((( 2**0.5/2+0.000j), (+2**0.5/2+0.000j), (+0.000+0.000j)),
(( 0.000+2**0.5/2j), (+0.000-2**0.5/2j), (+0.000+0.000j)),
(( 0.000+0.000j), (+0.000+0.000j), (+0.000+1.000j)))]:
printm('\nMatrix:', matrix)
ct = conjugate_transpose(matrix)
printm('Its conjugate transpose:', ct)
print('Hermitian? %s.' % ishermitian(matrix, ct))
print('Normal?  %s.' % isnormal(matrix, ct))
print('Unitary?  %s.' % isunitary(matrix, ct))
Output:
Matrix:
3+0j, 2+1j
2-1j, 1+0j
Its conjugate transpose:
3-0j, 2+1j
2-1j, 1-0j
Hermitian? True.
Normal?    True.
Unitary?   False.

Matrix:
1+0j, 1+0j, 0+0j
0+0j, 1+0j, 1+0j
1+0j, 0+0j, 1+0j
Its conjugate transpose:
1-0j, 0-0j, 1-0j
1-0j, 1-0j, 0-0j
0-0j, 1-0j, 1-0j
Hermitian? False.
Normal?    True.
Unitary?   False.

Matrix:
0.707107+0j, 0.707107+0j,        0+0j
0-0.707107j, 0+0.707107j,        0+0j
       0+0j,        0+0j,        0+1j
Its conjugate transpose:
0.707107-0j, 0+0.707107j,        0-0j
0.707107-0j, 0-0.707107j,        0-0j
       0-0j,        0-0j,        0-1j
Hermitian? False.
Normal?    True.
Unitary?   True.

[edit] Racket

 
#lang racket
(require math)
(define H matrix-hermitian)
 
(define (normal? M)
(define MH (H M))
(equal? (matrix* MH M)
(matrix* M MH)))
 
(define (unitary? M)
(define MH (H M))
(and (matrix-identity? (matrix* MH M))
(matrix-identity? (matrix* M MH))))
 
(define (hermitian? M)
(equal? (H M) M))
 

Test:

 
(define M (matrix [[3.000+0.000i +2.000+1.000i]
[2.000-1.000i +1.000+0.000i]]))
(H M)
(normal? M)
(unitary? M)
(hermitian? M)
 

Output:

 
(array #[#[3.0-0.0i 2.0+1.0i] #[2.0-1.0i 1.0-0.0i]])
#t
#f
#f
 

[edit] REXX

/*REXX pgm performs a  conjugate transpose  on a  complex square matrix.*/
parse arg N elements; if N=='' then N=3
M.=0 /*Matrix has all elements = zero.*/
k=0; do r=1 for N
do c=1 for N; k=k+1; M.r.c=word(word(elements,k) 1,1); end /*c*/
end /*r*/
 
call showCmat 'M' ,N /*display a nice formatted matrix*/
identity.=0; do d=1 for N; identity.d.d=1; end /*d*/
call conjCmat 'MH', "M" ,N /*conjugate the M matrix ───► MH */
call showCmat 'MH' ,N /*display a nice formatted matrix*/
say 'M is Hermitian: ' word('no yes',isHermitian('M',"MH",N)+1)
call multCmat 'M', 'MH', 'MMH', N /*multiple two matrices together.*/
call multCmat 'MH', 'M', 'MHM', N /* " " " " */
say ' M is Normal: ' word('no yes',isHermitian('MMH',"MHM",N)+1)
say ' M is Unary: ' word('no yes',isUnary('M',N)+1)
say 'MMH is Unary: ' word('no yes',isUnary('MMH',N)+1)
say 'MHM is Unary: ' word('no yes',isUnary('MHM',N)+1)
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────CONJCMAT subroutine─────────────────*/
conjCmat: arg matX,matY,rows 1 cols; call normCmat matY,rows
do r=1 for rows; _=
do c=1 for cols; v=value(matY'.'r"."c)
rP=rP(v); cP=-cP(v); call value matX'.'c"."r, rP','cP
end /*c*/
end /*r*/
return
/*──────────────────────────────────ISHERMITIAN subroutine──────────────*/
isHermitian: arg matX,matY,rows 1 cols; call normCmat matX,rows
call normCmat matY,rows
do r=1 for rows; _=
do c=1 for cols
if value(matX'.'r"."c)\=value(matY'.'r"."c) then return 0
end /*c*/
end /*r*/
return 1
/*──────────────────────────────────ISUNARY subroutine──────────────────*/
isUnary: arg matX,rows 1 cols
do r=1 for rows; _=
do c=1 for cols; z=value(matX'.'r"."c); rP=rP(z); cP=cP(z)
if abs(sqrt(rP(z)**2+cP(z)**2)-(r==c))>=.0001 then return 0
end /*c*/
end /*r*/
return 1
/*──────────────────────────────────MULTCMAT subroutine─────────────────*/
multCmat: arg matA,matB,matT,rows 1 cols; call value matT'.',0
do r=1 for rows; _=
do c=1 for cols
do k=1 for cols; T=value(matT'.'r"."c); Tr=rP(T); Tc=cP(T)
A=value(matA'.'r"."k); Ar=rP(A); Ac=cP(A)
B=value(matB'.'k"."c); Br=rP(B); Bc=cP(B)
Pr=Ar*Br-Ac*Bc; Pc=Ac*Br+Ar*Bc; Tr=Tr+Pr; Tc=Tc+Pc
call value matT'.'r"."c,Tr','Tc
end /*k*/
end /*c*/
end /*r*/
return
/*──────────────────────────────────NORMCMAT subroutine─────────────────*/
normCmat: arg matN,rows 1 cols
do r=1 to rows; _=
do c=1 to cols; v=translate(value(matN'.'r"."c),,"IiJj")
parse upper var v real ',' cplx
if real\=='' then real=real/1
if cplx\=='' then cplx=cplx/1; if cplx=0 then cplx=
if cplx\=='' then cplx=cplx"j"
call value matN'.'r"."c,strip(real','cplx,"T",',')
end /*c*/
end /*r*/
return
/*──────────────────────────────────SHOWCMAT subroutine─────────────────*/
showCmat: arg matX,rows,cols; if cols=='' then cols=rows; pad=left('',6)
say; say center('matrix' matX,79,'─'); call normCmat matX,rows,cols
do r=1 to rows; _=
do c=1 to cols; _=_ pad left(value(matX'.'r"."c),9)
end /*c*/
say _
end /*r*/
say; return
/*──────────────────────────────────one─liner subroutines───────────────*/
cP: procedure; arg ',' p; return word(strip(translate(p,,'IJ')) 0,1)
rP: procedure; parse arg r ','; return word(r 0,1)
/*──────────────────────────────────SQRT subroutine─────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits()
numeric digits 11; g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g);end
numeric digits d; return g/1
.sqrtGuess: numeric form; m.=11; p=d+d%4+2
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2

output when using the default input

───────────────────────────────────matrix M────────────────────────────────────
        1                1                1
        1                1                1
        1                1                1


───────────────────────────────────matrix MH───────────────────────────────────
        1                1                1
        1                1                1
        1                1                1

M is Hermitian:   yes
  M is Normal:    yes
  M is Unary:     no
MMH is Unary:     no
MHM is Unary:     no

output when using the input of: 3 .7071 .7071 0 0,.7071 0-.7071 0 0 0 0,1

───────────────────────────────────matrix M────────────────────────────────────
        0.7071           0.7071           0
        0,0.7071j        0,-0.7071        0
        0                0                0,1j


───────────────────────────────────matrix MH───────────────────────────────────
        0.7071           0,-0.7071        0
        0.7071           0,0.7071j        0
        0                0                0,-1j

M is Hermitian:   no
  M is Normal:    yes
  M is Unary:     no
MMH is Unary:     yes
MHM is Unary:     yes

[edit] Ruby

Works with: Ruby version 2.0
require 'matrix'
 
# Start with some matrix.
i = Complex::I
matrix = Matrix[[i, 0, 0],
[0, i, 0],
[0, 0, i]]
 
# Find the conjugate transpose.
# Matrix#conjugate appeared in Ruby 1.9.2.
conjt = matrix.conj.t # aliases for matrix.conjugate.tranpose
print 'conjugate tranpose: '; puts conjt
 
if matrix.square?
# These predicates appeared in Ruby 1.9.3.
print 'Hermitian? '; puts matrix.hermitian?
print ' normal? '; puts matrix.normal?
print ' unitary? '; puts matrix.unitary?
else
# Matrix is not square. These predicates would
# raise ExceptionForMatrix::ErrDimensionMismatch.
print 'Hermitian? false'
print ' normal? false'
print ' unitary? false'
end

Note: Ruby 1.9 had a bug in the Matrix#hermitian? method. It's fixed in 2.0.

[edit] Sparkling

Sparkling has support for basic complex algebraic operations, but complex matrix operations are not in the standard library.

# Computes conjugate transpose of M
let conjTransp = function conjTransp(M) {
return map(range(sizeof M[0]), function(row) {
return map(range(sizeof M), function(col) {
return cplx_conj(M[col][row]);
});
});
};
 
# Helper for cplxMatMul
let cplxVecScalarMul = function cplxVecScalarMul(A, B, row, col) {
var M = { "re": 0.0, "im": 0.0 };
let N = sizeof A;
for (var i = 0; i < N; i++) {
let P = cplx_mul(A[row][i], B[i][col]);
M = cplx_add(M, P);
}
return M;
};
 
# Multiplies matrices A and B
# A and B are assumed to be square and of the same size,
# this condition is not checked.
let cplxMatMul = function cplxMatMul(A, B) {
var R = {};
let N = sizeof A;
for (var row = 0; row < N; row++) {
R[row] = {};
for (var col = 0; col < N; col++) {
R[row][col] = cplxVecScalarMul(A, B, row, col);
}
}
return R;
};
 
# Helper for creating an array representing a complex number
# given its textual representation
let _ = function makeComplex(str) {
let sep = indexof(str, "+", 1);
if sep < 0 {
sep = indexof(str, "-", 1);
}
let reStr = substrto(str, sep);
let imStr = substrfrom(str, sep);
return { "re": tofloat(reStr), "im": tofloat(imStr) };
};
 
# Formats a complex matrix
let printCplxMat = function printCplxMat(M) {
foreach(M, function(i, row) {
foreach(row, function(j, elem) {
printf("  %.2f%+.2fi", elem.re, elem.im);
});
print();
});
};
 
# A Hermitian matrix
let H = {
{ _("3+0i"), _("2+1i") },
{ _("2-1i"), _("0+0i") }
};
 
# A normal matrix
let N = {
{ _("1+0i"), _("1+0i"), _("0+0i") },
{ _("0+0i"), _("1+0i"), _("1+0i") },
{ _("1+0i"), _("0+0i"), _("1+0i") }
};
 
# A unitary matrix
let U = {
{ _("0.70710678118+0i"), _("0.70710678118+0i"), _("0+0i") },
{ _("0-0.70710678118i"), _("0+0.70710678118i"), _("0+0i") },
{ _("0+0i"), _("0+0i"), _("0+1i") }
};
 
 
print("Hermitian matrix:\nH = ");
printCplxMat(H);
print("H* = ");
printCplxMat(conjTransp(H));
print();
 
print("Normal matrix:\nN = ");
printCplxMat(N);
print("N* = ");
printCplxMat(conjTransp(N));
print("N* x N = ");
printCplxMat(cplxMatMul(conjTransp(N), N));
print("N x N* = ");
printCplxMat(cplxMatMul(N, conjTransp(N)));
print();
 
print("Unitary matrix:\nU = ");
printCplxMat(U);
print("U* = ");
printCplxMat(conjTransp(U));
print("U x U* = ");
printCplxMat(cplxMatMul(U, conjTransp(U)));
print();

[edit] Tcl

Tcl's matrixes (in Tcllib) do not assume that the contents are numeric at all. As such, they do not provide mathematical operations over them and this considerably increases the complexity of the code below. Note the use of lambda terms to simplify access to the complex number package.

Library: Tcllib (Package: math::complexnumbers)
Library: Tcllib (Package: struct::matrix)
package require struct::matrix
package require math::complexnumbers
 
proc complexMatrix.equal {m1 m2 {epsilon 1e-14}} {
if {[$m1 rows] != [$m2 rows] || [$m1 columns] != [$m2 columns]} {
return 0
}
# Compute the magnitude of the difference between two complex numbers
set ceq [list apply {{epsilon a b} {
expr {[mod [- $a $b]] < $epsilon}
} ::math::complexnumbers} $epsilon]
for {set i 0} {$i<[$m1 columns]} {incr i} {
for {set j 0} {$j<[$m1 rows]} {incr j} {
if {![{*}$ceq [$m1 get cell $i $j] [$m2 get cell $i $j]]} {
return 0
}
}
}
return 1
}
 
proc complexMatrix.multiply {a b} {
if {[$a columns] != [$b rows]} {
error "incompatible sizes"
}
# Simplest to use a lambda in the complex NS
set cpm {{sum a b} {
+ $sum [* $a $b]
} ::math::complexnumbers}
set c0 [math::complexnumbers::complex 0.0 0.0]; # Complex zero
set c [struct::matrix]
$c add columns [$b columns]
$c add rows [$a rows]
for {set i 0} {$i < [$a rows]} {incr i} {
for {set j 0} {$j < [$b columns]} {incr j} {
set sum $c0
foreach rv [$a get row $i] cv [$b get column $j] {
set sum [apply $cpm $sum $rv $cv]
}
$c set cell $j $i $sum
}
}
return $c
}
 
proc complexMatrix.conjugateTranspose {matrix} {
set mat [struct::matrix]
$mat = $matrix
$mat transpose
for {set c 0} {$c < [$mat columns]} {incr c} {
for {set r 0} {$r < [$mat rows]} {incr r} {
set val [$mat get cell $c $r]
$mat set cell $c $r [math::complexnumbers::conj $val]
}
}
return $mat
}

Using these tools to test for the properties described in the task:

proc isHermitian {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
return 0
}
set cc [complexMatrix.conjugateTranspose $matrix]
set result [complexMatrix.equal $matrix $cc $epsilon]
$cc destroy
return $result
}
 
proc isNormal {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
return 0
}
set mh [complexMatrix.conjugateTranspose $matrix]
set mhm [complexMatrix.multiply $mh $matrix]
set mmh [complexMatrix.multiply $matrix $mh]
$mh destroy
set result [complexMatrix.equal $mhm $mmh $epsilon]
$mhm destroy
$mmh destroy
return $result
}
 
proc isUnitary {matrix {epsilon 1e-14}} {
if {[$matrix rows] != [$matrix columns]} {
# Must be square!
return 0
}
set mh [complexMatrix.conjugateTranspose $matrix]
set mhm [complexMatrix.multiply $mh $matrix]
set mmh [complexMatrix.multiply $matrix $mh]
$mh destroy
set result [complexMatrix.equal $mhm $mmh $epsilon]
$mhm destroy
if {$result} {
set id [struct::matrix]
$id = $matrix; # Just for its dimensions
for {set c 0} {$c < [$id columns]} {incr c} {
for {set r 0} {$r < [$id rows]} {incr r} {
$id set cell $c $r \
[math::complexnumbers::complex [expr {$c==$r}] 0]
}
}
set result [complexMatrix.equal $mmh $id $epsilon]
$id destroy
}
$mmh destroy
return $result
}
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox