Anonymous user
Conjugate transpose: Difference between revisions
+ D entry
m (→{{header|FORTRAN}}: Corrected header.) |
(+ D entry) |
||
Line 302:
Complex Matrix A is not normal</pre>
=={{header|D}}==
{{trans|Python}}
<lang d>import std.stdio, std.complex, std.math, std.range, std.algorithm,
std.numeric;
// alias CM(T) = Complex!T[][]; // Not yet useful.
Complex!T[][] conjugateTranspose(T)(in Complex!T[][] m)
/*pure nothrow*/ {
return iota(m[0].length)
.map!(i => m.transversal(i).map!conj.array)
.array;
}
T[][] matMul(T)(in T[][] A, in T[][] B) /*pure nothrow*/ {
const Bt = B[0].length.iota.map!(i=> B.transversal(i).array).array;
return A.map!(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 {
return v[1 .. $].all!(c => c == v[0]);
}
bool allNearSame(in Complex!T[] v) pure nothrow {
Complex!T v0 = v[0]; // To avoid another cast.
return v[1 .. $].all!(c=> feqrel(v0.re, cast()c.re) >= nBits &&
feqrel(v0.im, cast()c.im) >= nBits);
}
/*immutable*/ const x = matrices.map!(m => m.length).array;
if (!allSame(x))
return false;
/*immutable*/ const 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 iota(side)
.map!(r => side.iota.map!(c => complex(r == c)).array)
.array;
}
bool isUnitary(T)(in Complex!T[][] m, in Complex!T[][] ct)
/*pure nothrow*/ {
const mct = matMul(m, ct);
const 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)]]
]) {
writefln("Matrix:\n[%([%(%s, %)],\n %)]]", matrix);
/*immutable*/ const ct = conjugateTranspose(matrix);
"Its conjugate transpose:".writeln;
writefln("[%([%(%s, %)],\n %)]]", ct);
writefln("Hermitian? %s.", isHermitian(matrix, ct));
writefln("Normal? %s.", isNormal(matrix, ct));
writefln("Unitary? %s.\n", isUnitary(matrix, ct));
}
}</lang>
{{out}}
<pre>Matrix:
[[3+0i, 2+1i],
[2-1i, 1+0i]]
Its conjugate transpose:
[[3-0i, 2+1i],
[2-1i, 1-0i]]
Hermitian? true.
Normal? true.
Unitary? false.
Matrix:
[[1+0i, 1+0i, 0+0i],
[0+0i, 1+0i, 1+0i],
[1+0i, 0+0i, 1+0i]]
Its conjugate transpose:
[[1-0i, 0-0i, 1-0i],
[1-0i, 1-0i, 0-0i],
[0-0i, 1-0i, 1-0i]]
Hermitian? false.
Normal? true.
Unitary? false.
Matrix:
[[0.707107+0i, 0.707107+0i, 0+0i],
[0-0.707107i, 0+0.707107i, 0+0i],
[0+0i, 0+0i, 0+1i]]
Its conjugate transpose:
[[0.707107-0i, 0+0.707107i, 0-0i],
[0.707107-0i, 0-0.707107i, 0-0i],
[0-0i, 0-0i, 0-1i]]
Hermitian? false.
Normal? true.
Unitary? true.
</pre>
=={{header|Factor}}==
| |||