       Re: DirectSum (feature request)

• To: mathgroup at smc.vnet.net
• Subject: [mg97818] Re: DirectSum (feature request)
• From: Raffy <raffy at mac.com>
• Date: Sun, 22 Mar 2009 05:50:17 -0500 (EST)
• References: <gq2eu9\$ech\$1@smc.vnet.net>

```On Mar 21, 3:18 am, Maris Ozols <maroz... at gmail.com> wrote:
> Taking a direct sum of a given list of matrices is a very common task
> (unless you are a quantum physicist and use only KroneckerProduct).
> Unfortunately there is no built-in function (that I know of) for doing
> this in Mathematica. The closest thing we have is ArrayFlatten. So I
> usually do something like this to compute a direct sum:
>
> DirectSum[Ms_] := Module[{n = Length[Ms], z, i},
>   z = ConstantArray[0, n];
>   ArrayFlatten@Table[ReplacePart[z, i -> Ms[[i]]], {i, 1, n}]
> ];
>
> Is there a better way of doing this?
>
> Note: A nice way to implement it would be
>
> DirectSum[Ms_] := ArrayFlatten@DiagonalMatrix[Ms];
>
> Unfortunately this gives "DiagonalMatrix::vector" error, since
> DiagonalMatrix is not flexible enough to accept a list of matrices.
> The way DiagonalMatrix is used in the above code might cause some
> confusion for beginners, but in general I don't see why DiagonalMatrix
> should be limited in this way.
>
> ~Maris Ozols~

directSum1[vComp_?(VectorQ[#, MatrixQ] &)] := Block[
{mForm = ConstantArray[0, {Length[vComp], Length[vComp]}]},
Do[mForm[[i, i]] = vComp[[i]], {i, Length[vComp]}];
ArrayFlatten[mForm]];

directSum2[vComp_?(VectorQ[#, MatrixQ] &)] :=
ArrayFlatten@
Array[If[Equal[##], vComp[[#1]], 0] &, {Length[vComp],
Length[vComp]}];

directSum3[vComp_?(VectorQ[#, MatrixQ] &)] :=
ArrayFlatten[
DiagonalMatrix[Range@Length[vComp]] /. i_ /; i > 0 :> vComp[[i]]];

directSum4[vComp_?(VectorQ[#, MatrixQ] &)] :=
ArrayFlatten@ReleaseHold@DiagonalMatrix[Hold /@ vSum];

directSum5[vComp_?(VectorQ[#, MatrixQ] &)] :=
ArrayFlatten[DiagonalMatrix[\$ /@ vSum] /. \$ -> Sequence];

vSum= {{{a, b}, {c, d}}, {{n}}, {{r, s, t}, {u, v, w}, {x, y, z}}};

Out= {0.025205, Null}
In:= Do[r1 = directSum1[vSum], {10000}] // Timing
Out= {0.751228, Null}
In:= Do[r2 = directSum2[vSum], {10000}] // Timing
Out= {0.79516, Null}
In:= Do[r3 = directSum3[vSum], {10000}] // Timing
Out= {0.857756, Null}
In:= Do[r4 = directSum4[vSum], {10000}] // Timing
Out= {0.65793, Null}
In:= Do[r5 = directSum5[vSum], {10000}] // Timing
Out= {0.71162, Null}

directSum4 seems to be the fastest, it is also the simplest and most
straightforward.
The Hold/ReleaseHold is used to cheat the DiagonalMatrix VectorQ check.

```

• Prev by Date: Re: Using findroot and NIntegrate with a double integral
• Next by Date: Re: ParametricPlot problem
• Previous by thread: Re: DirectSum (feature request)
• Next by thread: Re: DirectSum (feature request)