13
$\begingroup$

I have a Real packed array arr. It may have arbitrary depth but Last@Dimensions[arr] == 2. For example,

arr = N@Partition[Partition[Range[24], {2}], 3]

I want to convert it to a Complex packed array of dimensions Most@Dimensions[arr].

Thus {{1., 2.}, {3., 4.}} would be converted to {1. + 2. I, 3. + 4. I}.

What is the most efficient way of doing so? The key here is memory efficiency and speed.


One solution is

toComplex[
  arr_ /; 
    Developer`PackedArrayQ[arr, Real] && 
    Last@Dimensions[arr] == 2
 ] := #1 + I #2 & @@ Transpose[arr, Append[Range[2, Depth[arr] - 1], 1]]

This unpacks only to level 1, and does some arithmetic.

Can we avoid any unpacking and any arithmetic? LibraryLink's MTensors, which I suspect are very similar to the internal representation of packed arrays, would store both the Complex and the corresponding Real arrays identically in memory. Converting between them is possible without any unpacking or arithmetic. Only simple copying would be needed. Is this also possible in Mathematica?


Here's a small benchmark for the posted solutions. I am also comparing a LibraryLink-based solution.

Up to an array size limit, Carl Woll's solution is the fastest. For very large arrays, it performs the same as user21's solution.

arr = RandomReal[1, {13, 12, 11, 10, 9, 8, 2}];

r2 = toComplexU21[arr]; // RepeatedTiming
(* {0.026, Null} *)

r3 = toComplexLLIAM[arr]; // RepeatedTiming
(* {0.17, Null} *)

r4 = toComplexCW[arr]; // RepeatedTiming
(* {0.017, Null} *)

r5 = toComplexLibLink[arr]; // RepeatedTiming
(* {0.0050, Null} *)

r2 == r3 == r4 == r5
(* True *)

Larger arrays:

arr = RandomReal[1, {16, 15, 14, 13, 12, 11, 2}];

r2 = toComplexU21[arr]; // RepeatedTiming
(* {0.141, Null} *)

r3 = toComplexLLIAM[arr]; // RepeatedTiming
(* {0.78, Null} *)

r4 = toComplexCW[arr]; // RepeatedTiming
(* {0.141, Null} *)

r5 = toComplexLibLink[arr]; // RepeatedTiming
(* {0.032, Null} *)

The code:

c = Developer`ToPackedArray[{1, I}, Complex];
toComplexCW[arr_] := arr.c;

toComplexU21[arr_] := #[[1]] + #[[2]]*I &[
  Transpose[arr, Append[Range[2, Depth[arr] - 1], 1]]]

toComplexSz[arr_] := #1 + I #2 & @@ 
  Transpose[arr, Append[Range[2, Depth[arr] - 1], 1]]

toComplexLLIAM = 
  Compile[{{pair, _Real, 1}}, First@pair + I Last@pair, 
   RuntimeAttributes -> {Listable}];

(* This function has no safety checks.  Be careful about what arrays 
   you are passing to it, as it WILL crash the kernel if
   given invalid input! *)
src = "
  #include \"WolframLibrary.h\"
  #include \"string.h\"

  DLLEXPORT mint WolframLibrary_getVersion(){ return WolframLibraryVersion; }

  DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) { return 0; }

  DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) { return; }

  DLLEXPORT int toComplex(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res){
     MTensor rt = MArgument_getMTensor(Args[0]);
     MTensor ct;
     libData->MTensor_new(MType_Complex, libData->MTensor_getRank(rt) - 1, libData->MTensor_getDimensions(rt), &ct);
     memcpy((void *) libData->MTensor_getComplexData(ct), (void *) libData->MTensor_getRealData(rt), libData->MTensor_getFlattenedLength(rt) * sizeof(double));
     MArgument_setMTensor(Res, ct);
     return LIBRARY_NO_ERROR;
  }
  ";

Needs["CCompilerDriver`"]
lib = CreateLibrary[src, "toComplex"];

toComplexLibLink = 
  LibraryFunctionLoad[lib, 
   "toComplex", {{Real, _, "Constant"}}, {Complex, _}];
$\endgroup$
  • $\begingroup$ Very nice. Favorited as a reference for LibraryLink. $\endgroup$ – LLlAMnYP Oct 18 '17 at 14:10
10
$\begingroup$

Using Dot is pretty fast:

a = RandomReal[1, {8, 7, 6, 5, 4, 3, 2}];
c = N[{1, I}];

r1 = a.c;// RepeatedTiming

{0.000100, Null}

$\endgroup$
  • $\begingroup$ Looks like the way to go for tensors. Mine works for non-tensor objects, but then starts unpacking seemingly deeper than necessary. $\endgroup$ – LLlAMnYP Oct 18 '17 at 13:45
6
$\begingroup$

Get rid of the Apply:

#[[1]] + #[[2]]*I &[
 Transpose[arr, Append[Range[2, Depth[arr] - 1], 1]]]
$\endgroup$
4
$\begingroup$
toCmplx = 
 Compile[{{pair, _Real, 1}}, First@pair + I Last@pair, 
  RuntimeAttributes -> {Listable}]

I believe this works with arbitrary depth, including ragged arrays, so long as the bottom dimension is 2 for all elements. Primitive test with On["Packing"] returned no messages.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.