Conjugate transpose: Difference between revisions
(PL/I version created */) |
(→{{header|REXX}}: added the REXX language. -- ~~~~) |
||
| Line 624: | Line 624: | ||
Matrix is unitary |
Matrix is unitary |
||
Matrix is unitary |
Matrix is unitary |
||
</pre> |
|||
=={{header|REXX}}== |
|||
<lang 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</lang> |
|||
'''output''' when using the default input |
|||
<pre style="overflow:scroll"> |
|||
───────────────────────────────────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 |
|||
</pre> |
|||
'''output''' when using the input of: <tt> 3 .7071 .7071 0 0,.7071 0-.7071 0 0 0 0,1 </tt> |
|||
<pre style="overflow:scroll"> |
|||
───────────────────────────────────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 |
|||
</pre> |
</pre> |
||
Revision as of 03:33, 27 December 2012
Suppose that a matrix contains complex numbers. Then the conjugate transpose of is a matrix containing the complex conjugates of the matrix transposition of .
This means that row , column of the conjugate transpose equals the complex conjugate of row , column of the original matrix.
In the next list, must also be a square matrix.
- A Hermitian matrix equals its own conjugate transpose: .
- A normal matrix is commutative in multiplication with its conjugate transpose: .
- A unitary matrix has its inverse equal to its conjugate transpose: . This is true iff and iff , where is the identity 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.
- MathWorld: conjugate transpose, Hermitian matrix, normal matrix, unitary matrix
Ada
<lang 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;</lang>
- 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
C
<lang 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 = (complex **) malloc (b.rows * sizeof (complex *));
for (i = 0; i < b.rows; i++)
{
b.z[i] = (complex *) 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 = (complex **) malloc (c.rows * (sizeof (complex *)));
for (i = 0; i < c.rows; i++)
{
c.z[i] = (complex *) 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 = (complex **) malloc (a.rows * sizeof (complex *));
printf ("Randomly Generated Complex Matrix A is : ");
for (i = 0; i < a.rows; i++)
{
printf ("\n");
a.z[i] = (complex *) 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] = (complex *) 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;
} </lang>
- 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
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.
<lang factor>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 = ;</lang>
Usage:
USE: rosetta.hermitian
IN: scratchpad { { C{ 1 2 } 0 }
{ 0 C{ 3 4 } } }
[ hermitian-matrix? . ]
[ normal-matrix? . ]
[ unitary-matrix? . ] tri
f
t
f
Go
<lang 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
}</lang> 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
J
Solution: <lang j> ct =: +@|: NB. Conjugate transpose (ct A is A_ct)</lang> Examples: <lang j> 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</lang>
Reference (example matrices for other langs to use):<lang j> 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)</lang>
Mathematica
<lang 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}</lang>
PL/I
<lang 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; </lang> 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
REXX
<lang 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</lang>
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
Ruby
This example is incorrect. Please fix the code and remove this message.
|
<lang ruby>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</lang>
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.
<lang tcl>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
}</lang> Using these tools to test for the properties described in the task: <lang tcl>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
}</lang>