Conjugate transpose: Difference between revisions

{{trans|Python}} A well typed and mostly imperative version:

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 {

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;

pure nothrow {

pure nothrow {

immutable mct = matMul(m, ct);

immutable ident = mct.length.complexIdentitymatrix;

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) {

writefln("Matrix:\n" ~ mFormat, mat);

immutable ct = conjugateTranspose(mat);

writefln("Hermitian? %s.", isHermitian(mat, ct));

writefln("Normal? %s.", isNormal(mat, ct));

writefln("Unitary? %s.\n", isUnitary(mat, ct));

===Alternative Version===

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

<lang d>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 {

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 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));

}

}</lang>