19
$\begingroup$

For odd n I'm looking for a short and swift way to construct with (f.e.)

n = 3

enter image description here

n = 11

enter image description here

n=1111

enter image description here

$\endgroup$
4
  • 4
    $\begingroup$ How come the question got fewer upvotes than it received answers? It should have at least that number (plus one). $\endgroup$
    – Jens
    Jun 30, 2014 at 0:01
  • 1
    $\begingroup$ @Jens ......... Why? $\endgroup$
    – eldo
    Jun 30, 2014 at 0:15
  • 2
    $\begingroup$ @eldo Well, if one deems a question worthy of answering, then it should probably be worthy of an upvote too. I upvote questions I answer almost always, though I often forget to do it right away, so the upvote may come only an hour after the answer. This is a question that clearly interests people and they're having fun with it. $\endgroup$
    – Szabolcs
    Jun 30, 2014 at 0:21
  • 1
    $\begingroup$ Related link: mathematica.stackexchange.com/questions/32307/… $\endgroup$
    – chyanog
    Jul 1, 2014 at 6:09

11 Answers 11

18
$\begingroup$

Here's my take using NestList

cm[n_] := NestList[# + 1 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], n - 1]

Then

cm[11]

Mathematica graphics

Here's a FoldList version (just as fast):

cmf[n_] := FoldList[#1 + #2 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], 
            ConstantArray[1, n - 1]]

The above methods according to the benchmarks posted are already as fast as the fastest methods. Here's a method that's just as fast as Belisarius's (I took a page from his rule based solution):

a4[n_] := With[{t = Join[Range[n/2 + 1], Reverse@Range[n/2]]}, Array[s + # &, n, 0] /. s -> t]

Timings:

I have taken the fastest methods from the five answers with the fastest times and for each function I averaged 7 runs per input for 9 different sizes (1001 to 9001), here are the results:

Mathematica graphics

Edit

After working with JacobAkerbbom to get his LibraryLink method working on Windows see here, I've now included his method which obviously is now the king of the hill.

$\endgroup$
8
  • $\begingroup$ @Give my frx update a whirl on large n if you have the time... $\endgroup$
    – ciao
    Jun 30, 2014 at 21:28
  • $\begingroup$ @rasher, that's the one I included in this latest update. The f6 was slightly slower than Szabolcs`. $\endgroup$
    – RunnyKine
    Jun 30, 2014 at 21:30
  • $\begingroup$ No, I mean the one called FRX - just added it ;-) $\endgroup$
    – ciao
    Jun 30, 2014 at 21:31
  • 1
    $\begingroup$ @rasher, I knew when you added it and I immediately tested and updated my plot. $\endgroup$
    – RunnyKine
    Jun 30, 2014 at 21:32
  • $\begingroup$ @rasher Like I said frx is slightly faster than Szabolcs` wheras f6 was slightly slower. $\endgroup$
    – RunnyKine
    Jun 30, 2014 at 21:34
12
$\begingroup$

Edit: See end of post for latest performance enhancement.

f=With[{c = Ceiling[#/2]}, c - 1 + Array[#1 - Abs[c - #2] &, {#, #}]] &;

f[5]

(* {{1, 2, 3, 2, 1}, {2, 3, 4, 3, 2}, {3, 4, 5, 4, 3}, {4, 5, 6, 5, 4}, {5, 6, 7, 6, 5}} *)

Short, sweet, fast.

For more speed,

f5 = With[{c = Ceiling[#/2]},
   Subtract[
    ArrayPad[ConstantArray[Range[#, # + c - 1], c], {{c - 1, 0}, {0, c - 1}}, "Reflected"], 
    Range[# - 1, 0, -1]]] &

And even faster...

f6[n_] := ConstantArray[Join[Range[n/2 + 1], Reverse@Range[n/2]], n] + Range[0, n - 1]

Update: For large n, this seems to do quite well (takes a hit for smallish n before auto-compile kicks in, could just compile overall):

frx[n_] := With[{rng = Range@(n + Floor[n/2]),fl = Floor[n/2], r = Range@n},
   Table[Join[rng[[x ;; x + fl]], rng[[Subtract[x + fl, 1] ;; x ;; -1]]], {x, r}]];

Finally, if this is something you're going to call repeatedly with varying n, particulary larger n, consider injecting a precomputed base of maximum needed size into a lookup function (see later fastest only benchmark for performance illustration):

x2 = With[{base = fr[5001]}, Drop[base, # + 1 ;;, {Ceiling[#/2] + 1, -Ceiling[#/2]}] &];

As requested, a quick bench chart. Usual loungebook caveats apply, and I did not include nasser's or aky's, both get slow on large problems. Looks like f6 and Runnykine`s solutions trade blows for fastest. Update: I'd spaced on Szabolc's answer, added it - seems it and F6 vie for fastest (I'd write off differences to noise), both having a slight but consistent edge over Runnykine's solution (and that itself is small).

Update 2: Added belisarius' second method. Quite quick!

enter image description here

Here's a new set of benchmarks of just the fastest few. Note I've used Timing vs AbsoluteTiming - I think this a much better indicator of efficiency for this kind of test. That said, the earlier ones seemed pretty consistent with Szabolcs' benchmarks, so beyond platform / version / AbsoluteTiming clouding things, hypothesis non fingo why Runnykine's benchmarks seem to differ significantly. Times are average of ten runs per size per tested.

Code I used for this Q-N-D:

tested = {x2,f6, frx, szabolcs, belisaurius2, runnykinea4};
times = 10;
cur = {}
Monitor[
 fastestBmark = 
  Table[cur = {size, (Mean[
         Table[(ClearSystemCache[]; 
           First@Timing@(#@size)), {times}]]) & /@ tested}, {size, 
    301, 5101, 200}], cur]

enter image description here

Finally (really...), one can extend the latter idea if the need is repeated calls over a wide range of n by caching "touchstone" base cases and using the appropriate case for a given n (since dropping has overhead inverse to amount dropped - one could also switch between Drop and Part depending on n vs touchstone size):

x3 = With[{base = 
     frx[#] & /@ {101, 301, 501, 701, 1001, 2001, 3001, 4001, 5001}}, 
   With[{p = Position[{101, 301, 501, 701, 1001, 2001, 3001, 4001, 5001}, 
                       x_ /; x >= #, 1, 1] &@#}, 
     If[p == {}, {}, 
          Drop[base[[p[[1, 1]]]], # + 1 ;;, {Ceiling[#/2] + 1, -Ceiling[#/2]}]]] &];

Here's a quick test against the fastest answer excluding my others:

enter image description here

$\endgroup$
17
  • $\begingroup$ We posted nearly the same thing at almost the same time. $\endgroup$
    – Szabolcs
    Jun 29, 2014 at 22:05
  • $\begingroup$ @Szabolcs: LOL, like minds. I recently posted and answer that was exactly the same as one by Mr. Wizard. We got a chuckle, he posted a hilarious GIF reply. +1 on yours, of course ;-) $\endgroup$
    – ciao
    Jun 29, 2014 at 22:07
  • $\begingroup$ @rasher All the answers so far are much nicer than the ugly code with which I produced the question. Yesterday, you made this logarithmic ("Ferrari") speed test. It would certainly facilitate my acceptance of one of the answers if you would ran "Ferrari" again against this case. $\endgroup$
    – eldo
    Jun 29, 2014 at 22:11
  • $\begingroup$ I feel left out from the chart :( On my machine my solution is faster than either f6 or RunnyKine's solution, by a tiny margin. Algorithmically it's the same as your f, but for some reason it's the fastest in my benchmarking. Probably because Table auto-compiles, while Array (probably) doesn't. $\endgroup$
    – Szabolcs
    Jun 30, 2014 at 0:17
  • 2
    $\begingroup$ @JacobAkkerboom copying elements helps in the sense that the branch is between copying another (continuing the loop) and stopping (exiting the loop). Doing it another way will introduce extra conditionals and special cases that result in more mispredicts. Of course we want to keep the source range in L1 cache, while writing an entire cache line at a time prevents cache thrashing. And the memory access predictor only deals with quite simple access patterns, so it may not work if we do non-consecutive accesses. $\endgroup$ Jul 11, 2014 at 4:06
11
$\begingroup$

This is pretty straightforward and very easy to follow even for someone who just started learning Mathematica. This has its value when you need to read your code a year later, even if you're an experienced user.

n = 11;
k = (n + 1)/2;
row = k - Abs[k - Range[n]];
Table[row + i, {i, 0, n - 1}]

Should be fast enough for most application.


Benchmarks

The timings from the plot below are outdated. Several people have posted faster versions since, but I haven't had time to update the plot. I'll keep the benchmark code below in case someone wants to experiment with it.

nasser[n_] := Module[{mid},
  mid = Ceiling[n/2];
  SparseArray[{{i_, j_} /; j > mid, {i_, j_} /; j <= mid} :> {n - j + 
      i, i + j - 1}, {n, n}]
  ]

kuba1[n_] :=     
 ArrayPad[HankelMatrix[Range@n, 
    Range[n, n + Floor[n/2]]][[;; , ;; Ceiling[n/2]]], {{0, 0}, {0, 
    Floor[n/2]}}, "Reflected"]

kuba2[n_] := 
 ArrayPad[Array[Range[#, # + Floor[n/2]] &, n], {{0, 0}, {0, Floor[n/2]}}, "Reflected"]


kuba3[n_] := 
 Transpose[Range[n] + # & /@ Join[#, #[[-2 ;; 1 ;; -1]]] &@Range[0, Floor[n/2]]]

rasher1 = 
 With[{c = Ceiling[#/2]}, c - 1 + Array[#1 - Abs[c - #2] &, {#, #}]] &


rasher5 = 
  With[{c = Ceiling[#/2]}, 
    Subtract[
     ArrayPad[
      ConstantArray[Range[#, # + c - 1], c], {{c - 1, 0}, {0, c - 1}},
       "Reflected"], Range[# - 1, 0, -1]]] &;

rasher6[n_] := 
 ConstantArray[Join[Range[n/2 + 1], Reverse@Range[n/2]], n] + 
  Range[0, n - 1]

belisarius[n_] := (Join[#, Rest@Reverse@#] &@Range[n/2 + 1]) + # & /@ Range[0, n - 1]

belisarius2[n_] := With[{k = (Join[#, Rest@Reverse@#] &@Range[n/2 + 1])}, s + Range[0, n - 1] /. s -> k]

algohi[n_] := 
 Module[{k},
  k = Table[i, {i, 1 + #, n + #}] & /@ Range[0, n/2]; 
  Transpose@Join[Most@k, Reverse@k]
 ]

runnykine[n_] := 
 NestList[# + 1 &, Join[Range[n/2 + 1], Reverse@Range[n/2]], n - 1]

szabolcs[n_] := 
 Module[{k, row},
  k = (n + 1)/2;
  row = k - Abs[k - Range[n]];
  Table[row + i, {i, 0, n - 1}]
 ]

Verify they all do what they should:

functions = {"nasser", "kuba1", "kuba2", "kuba3", "rasher1", "rasher5", "rasher6", "belisarius", "belisarius2", "algohi", "runnykine", "szabolcs"};

Equal @@ Through[(ToExpression /@ functions)[11]]
(* True *)

Benchmark, be careful to choose odd n values only:

bench[fun_String] := 
 Module[{f = ToExpression[fun]}, 
  Table[{n, Min@Table[First@AbsoluteTiming[f[n]], {3}]}, {n, 2 Round[2^Range[4, 11, 1/2]] + 1}]
]

results = bench /@ functions;

ListLogLogPlot[results, Joined -> True, 
 PlotRange -> {10^-5, 1}, PlotLegends -> functions, 
 PlotMarkers -> Automatic, PlotStyle -> colours]

enter image description here

Observe how both my and Algohi's solution gets a sudden boost at a certain size threshold. This is due to Table automatically compiling its argument above SystemOptions["CompileOptions" -> "TableCompileLength"].

Further speedups are certainly possible by manually compiling some of the other solutions.


Numerical results from my machine, so you don't have to re-run it if you just want to re-plot it.

Caveat: This was run on a laptop. Laptop CPUs have the habit of throttling their frequency when they heat up, which may affect results.

{{{33, 0.006806}, {47, 0.008816}, {65, 0.013106}, {91, 
   0.022375}, {129, 0.044451}, {183, 0.096939}, {257, 0.196829}, {363,
    0.413977}, {513, 0.885362}, {725, 1.889206}, {1025, 
   4.141960}, {1449, 9.229286}, {2049, 20.484070}, {2897, 
   47.915846}, {4097, 123.148882}}, {{33, 0.000029}, {47, 
   0.000046}, {65, 0.000079}, {91, 0.000140}, {129, 0.000270}, {183, 
   0.000531}, {257, 0.001022}, {363, 0.002025}, {513, 0.004224}, {725,
    0.008757}, {1025, 0.019102}, {1449, 0.039411}, {2049, 
   0.078217}, {2897, 0.153589}, {4097, 0.338157}}, {{33, 
   0.000081}, {47, 0.000115}, {65, 0.000159}, {91, 0.000229}, {129, 
   0.000351}, {183, 0.000594}, {257, 0.002159}, {363, 0.004096}, {513,
    0.008323}, {725, 0.016285}, {1025, 0.035317}, {1449, 
   0.070283}, {2049, 0.149849}, {2897, 0.288016}, {4097, 
   0.607098}}, {{33, 0.000096}, {47, 0.000136}, {65, 0.000194}, {91, 
   0.000294}, {129, 0.000477}, {183, 0.000937}, {257, 0.001827}, {363,
    0.003632}, {513, 0.007367}, {725, 0.014935}, {1025, 
   0.034238}, {1449, 0.066544}, {2049, 0.135146}, {2897, 
   0.267405}, {4097, 0.665833}}, {{33, 0.000211}, {47, 0.000402}, {65,
    0.000751}, {91, 0.001488}, {129, 0.002925}, {183, 0.005780}, {257,
    0.011469}, {363, 0.022900}, {513, 0.045602}, {725, 
   0.091659}, {1025, 0.185086}, {1449, 0.370420}, {2049, 
   0.745876}, {2897, 1.493142}, {4097, 3.090481}}, {{33, 
   0.000041}, {47, 0.000070}, {65, 0.000123}, {91, 0.000231}, {129, 
   0.000453}, {183, 0.000896}, {257, 0.001742}, {363, 0.003533}, {513,
    0.006994}, {725, 0.013959}, {1025, 0.029870}, {1449, 
   0.059805}, {2049, 0.120957}, {2897, 0.242442}, {4097, 
   0.567620}}, {{33, 0.000023}, {47, 0.000033}, {65, 0.000054}, {91, 
   0.000090}, {129, 0.000155}, {183, 0.000297}, {257, 0.000567}, {363,
    0.001158}, {513, 0.002287}, {725, 0.005820}, {1025, 
   0.013108}, {1449, 0.025564}, {2049, 0.051499}, {2897, 
   0.100180}, {4097, 0.319559}}, {{33, 0.000198}, {47, 0.000282}, {65,
    0.000398}, {91, 0.000582}, {129, 0.000879}, {183, 0.001349}, {257,
    0.002032}, {363, 0.003323}, {513, 0.005294}, {725, 
   0.008768}, {1025, 0.016870}, {1449, 0.030577}, {2049, 
   0.082289}, {2897, 0.154232}, {4097, 0.285149}}, {{33, 
   0.000086}, {47, 0.000120}, {65, 0.000168}, {91, 0.000249}, {129, 
   0.000392}, {183, 0.000618}, {257, 0.001007}, {363, 0.001706}, {513,
    0.002993}, {725, 0.005540}, {1025, 0.011265}, {1449, 
   0.021607}, {2049, 0.068486}, {2897, 0.128332}, {4097, 
   0.242391}}, {{33, 0.000112}, {47, 0.000175}, {65, 0.000286}, {91, 
   0.000493}, {129, 0.000909}, {183, 0.001677}, {257, 0.001292}, {363,
    0.002553}, {513, 0.005325}, {725, 0.011771}, {1025, 
   0.028106}, {1449, 0.055791}, {2049, 0.110606}, {2897, 
   0.217962}, {4097, 0.534833}}, {{33, 0.000074}, {47, 0.000108}, {65,
    0.000138}, {91, 0.000206}, {129, 0.000188}, {183, 0.000375}, {257,
    0.000625}, {363, 0.001197}, {513, 0.002524}, {725, 
   0.006462}, {1025, 0.014223}, {1449, 0.027980}, {2049, 
   0.056060}, {2897, 0.111271}, {4097, 0.342107}}, {{33, 
   0.000087}, {47, 0.000109}, {65, 0.000152}, {91, 0.000204}, {129, 
   0.000330}, {183, 0.000533}, {257, 0.000695}, {363, 0.001227}, {513,
    0.002388}, {725, 0.006340}, {1025, 0.013557}, {1449, 
   0.025133}, {2049, 0.047835}, {2897, 0.094955}, {4097, 0.262013}}}
$\endgroup$
5
  • $\begingroup$ Do you mind to use my second try on the comparison? $\endgroup$ Jun 30, 2014 at 1:23
  • $\begingroup$ @belisarius Done. $\endgroup$
    – Szabolcs
    Jun 30, 2014 at 3:11
  • $\begingroup$ Thanks a lot. For some reason my machine hangs running your test suite. $\endgroup$ Jun 30, 2014 at 3:16
  • $\begingroup$ @belisarius It seems to use a fair bit of memory ... the monitor goes to 75% and I have 16 GB. Updated now with more reliable results (minimum of 3 runs). That's it for now, even though there could be plenty of improvements ... it's not really worth it. It takes a lot of time to run. $\endgroup$
    – Szabolcs
    Jun 30, 2014 at 3:29
  • $\begingroup$ Nicely done ....+1 $\endgroup$ Jun 30, 2014 at 3:49
10
$\begingroup$
n = 11;
mid = Ceiling[n/2];
mat = SparseArray[{{i_, j_} /; j > mid, {i_, j_} /; j <= mid}:>{n-j+i,i+j-1}, {n, n}];

MatrixForm@mat

Mathematica graphics

$\endgroup$
9
$\begingroup$
b[n_] := (Join[#, Rest@Reverse@#] &@Range[n/2 + 1]) + # & /@ Range[0, n-1]
b[11] // MatrixForm

Mathematica graphics

Edit

Enhanced for some speedup, and curiously enough, it competes well with the fastest answers so far (see @rasher's benchamrk):

bs[n_] := With[{k = (Join[#, Rest@Reverse@#] &@Range[n/2 + 1])},  s+Range[0, n-1] /.s-> k]
$\endgroup$
3
  • $\begingroup$ Please see my update. It's now just neck and neck with yours and I took a page from your rule based approach. $\endgroup$
    – RunnyKine
    Jun 30, 2014 at 21:29
  • $\begingroup$ @RunnyKine Ha! Replacements are usually quite slow artifacts, I still don't understand why they beat "better" approaches here :) $\endgroup$ Jun 30, 2014 at 23:04
  • 1
    $\begingroup$ My guess is that since we use a symbol the addition is not done and all the lists are generated with no computation done. Once the replacement is made a fast vector operation happens almost concertedly giving the speed-up. $\endgroup$
    – RunnyKine
    Jun 30, 2014 at 23:07
8
$\begingroup$
n = 11;
k = Table[i, {i, 1 + #, n + #}] & /@ Range[0, n/2];
(Transpose@Join[Most@k, Reverse@k]) // TableForm

or

n = 11;
Table[j, {i, 1, 
   n}, {j, (n - 1)/2 - Abs[Range[-(n - 1)/2, (n - 1)/2]] + 
    i}] // TableForm

enter image description here

$\endgroup$
7
$\begingroup$
n = 11;
ArrayPad[
  HankelMatrix[Range@n, Range[n, n + Floor[n/2]]][[;; , ;; Ceiling[n/2]]],
  {{0, 0}, {0, Floor[n/2]}}, 
  "Reflected"]

or

ArrayPad[
  Array[Range[#, # + Floor[n/2]] &, n],
  {{0, 0}, {0, Floor[n/2]}}, 
  "Reflected"]

or

Transpose[
 Range[n] + # & /@ Join[#, #[[-2 ;; 1 ;; -1]]] &@Range[0, Floor[n/2]] 
 ]
$\endgroup$
2
  • $\begingroup$ +1 on use of Hankel. Started with similar idea, saw your post so would be redundant. $\endgroup$
    – ciao
    Jun 29, 2014 at 22:04
  • $\begingroup$ @rasher thanks, but it is almost the longest one :P $\endgroup$
    – Kuba
    Jun 30, 2014 at 4:49
7
$\begingroup$

Nothing special here, but as there wasn't any solution using Outer, I thought I'd post this:

With[{n = 11}, (* adjust n *)
Outer[#1 + #2 - 1 &, Range[n], Range[n/2 + 1]~Join~Reverse@Range[n/2]]]
$\endgroup$
6
$\begingroup$

Since no answer used the "Extrapolated" option for ArrayPad, which in my opinion is very straightforward (although not as efficient as the top ones), here is my trying:

Clear[extrpol]
extrpol[n_?OddQ] := 
    With[{k = (n + 1)/2}, 
         ArrayPad[2 k - {# + 1, #} & @ {2, 1, 2}, k - {{2, 1}, 2}, "Extrapolated"]
        ]
$\endgroup$
2
  • $\begingroup$ Mine looked like ArrayPad[# - {{2, 1, 2}, {1, 0, 1}, {0, -1, 0}}, Floor[#/2] - 1, "Extrapolated"] &, dropped idea - too slow. +1 in any case! $\endgroup$
    – ciao
    Jul 1, 2014 at 22:07
  • $\begingroup$ @rasher Thanks. On my machine, it's as fast as Kuba's first approach. $\endgroup$
    – Silvia
    Jul 2, 2014 at 5:11
5
$\begingroup$
g[n_] := With[{m=(n+1)/2},Table[n+j-m-Abs[m-i],{i,n},{j,n}]]
g[11]//Transpose//MatrixForm

enter image description here

$\endgroup$
5
$\begingroup$
<< SymbolicC`
<< Developer`
<< CCompilerDriver`
<< CCodeGenerator`

Please don't mind these unnecessary abstractions.

type = "mint";

abstractFunctionName = "makeMatr";

mainFunctionName = abstractFunctionName <> "I_T";

argumentSingletonGetterFunctionName[type_String] := 
  StringJoin["MArgument_get", type];

getter = argumentSingletonGetterFunctionName["Integer"];

typeSpecWL = "MType_Integer";

dataGetterAbstractor[type_String] :=
 "MTensor_get" <> type <> "Data"

dataGetter = dataGetterAbstractor["Integer"];

Generate some SymbolicC

makeMatrSC =
 CFunction[
  "int",
  mainFunctionName,
  {{"WolframLibraryData", "libData"}, {"mint", 
    "Argc"}, {CPointerType["MArgument"], "Args"}, {"MArgument", "Res"}}
  ,
  CBlock[
   {
    CDeclare["int", CAssign["err", "LIBRARY_NO_ERROR"]],
    CDeclare[type, "input"],
    CDeclare["MTensor", "result"], 
    CDeclare[type, CArray["resultDimensions", 2]],

    CAssign[
     "input", 
     CCall[getter, CArray["Args", 0]]
     ],

    CAssign[CArray["resultDimensions", 0], "input"], 
    CAssign[CArray["resultDimensions", 1], "input"], 
    CAssign["err", 
     CCall[CPointerMember["libData", "MTensor_new"], {typeSpecWL, 2, 
       "resultDimensions", CAddress["result"]}]],

    CDeclare[CPointerType[type], "resultDataPtr"],
    CDeclare[type, CAssign["value", 1]],
    CDeclare[type, 
     CAssign["square", COperator[Times, {"input", "input"}]]],
    CDeclare[type, CAssign["half", COperator[Divide, {"input", 2}]]],

    CAssign["resultDataPtr", 
     CCall[CPointerMember["libData", dataGetter], {"result"}]], 
    CDeclare[type, CAssign["iter", 1]],
    CDeclare[type, CAssign["iterRest", 1]],
    CDeclare[type, 
     CAssign["inputMOne" , COperator[Minus, {"input", 1}]]],

    CWhile[
     COperator[LessEqual, {"iter", "square"}],
     CBlock[
      {
       CAssign[CDereference["resultDataPtr"], "value"],
       COperator[Increment, "resultDataPtr"],
       CIf[
        COperator[LessEqual, {"iterRest", "half"}],
        COperator[Increment, "value"],
        COperator[Decrement, "value"]
        ],
       COperator[Increment, "iter"],
       CIf[
        COperator[Equal, {"iterRest", "inputMOne" }],
        CAssign["iterRest", 0],
        COperator[Increment, "iterRest"]
        ]
       }
      ]
     ],
    CCall["MArgument_setMTensor", {"Res", "result"}],
    CReturn["err"]
    }
   ]
  ];

Make it into a string

cCodeString = "DLLEXPORT"<> " " <> ToCCodeString[makeMatrSC];

boilerPlate = "
#include \"WolframLibrary.h\"

/* Return the version of Library Link */
DLLEXPORT mint WolframLibrary_getVersion( ) {
\treturn WolframLibraryVersion;
}

/* Initialize Library */
DLLEXPORT int WolframLibrary_initialize( WolframLibraryData \
libData) {
\treturn LIBRARY_NO_ERROR;
}

/* Uninitialize Library */
DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData \
libData) {
\treturn;
}

";

totalCString = boilerPlate <> cCodeString;

Create a library (code can be reused, as it makes different versions of the library. Todo: make code that cleans up libraries)

If[! ValueQ[counter], counter = 1;];
counter++;
counterString = ToString[counter];
libraryName = abstractFunctionName <> "Lib" <> counterString;
lib = CreateLibrary[totalCString, libraryName];

Load the library

LibraryLoad[libraryName]

makeMatrLL = 
 LibraryFunctionLoad[libraryName, 
  mainFunctionName, {{Integer}}, {Integer, 2}]

and profit

makeMatrLL[5]

enter image description here

Timing

My previous timings were kind of bad. This is because rI = (rI+1)%input is much slower than

if(rI == input -1){rI = 0}else{rI++}`

So I guess I learned something :).

My code seems to be the fastest code in this Q&A. Below I compare with RunnyKine, but really you should see his answer for a nice graph.

Timing data

timeTable[func_] := 
 Table[Mean@Table[First@Timing@func[in], {5}], {in, 1001, 9001, 1000}]

timings = timeTable[makeMatrLL];
a4Timings = timeTable[a4];
{timings, a4Timings}
 {{0.002772, 0.009899, 0.021801, 0.038221, 0.109890, 0.159009, 
   0.213708, 0.284732, 0.377101}, {0.005047, 0.028429, 0.057490, 
   0.101514, 0.152389, 0.212622, 0.293324, 0.378144, 0.418646}}

Again, see RunnyKine's answer for a nice graph :)

Further improvements

It seems the main thing that is unnecessary in this code is that we check the value of restIter at every iteration. Maybe we can assume that the input size will always be larger than a certain number. In this case, we can inline some code that many times in order to avoid "branches", i.e. using CIf too much. Using a bunch of Ranges in Mathematica also has to deal with conditions when copying the Ranges into the big matrix. My LibraryLink code does not have to copy data, which is something normal Mathematica code cannot avoid. But it seems these conditionals take much more time than copying of data.

Finally Oleksandr made a good point that copying ranges from a single range array would probably be faster.

$\endgroup$
8
  • $\begingroup$ I was out of town, I just saw this now. will give it a spin and include the timing in my answer. +1, I was hoping someone would use LibraryLink since I'm not very good at it. $\endgroup$
    – RunnyKine
    Jul 6, 2014 at 20:29
  • $\begingroup$ I'm getting an error: "LibraryFunction::libload: The function makeMatrI_T was not loaded from the file C:\Users\... when I run the LibraryFunctionLoad $\endgroup$
    – RunnyKine
    Jul 6, 2014 at 20:47
  • $\begingroup$ @RunnyKine ah that is unfortunate. I suppose using CompilationTarget -> "C" in Compile does work for you? I hear it can be a bit tricky under Windows. $\endgroup$ Jul 6, 2014 at 22:40
  • $\begingroup$ Yes that works for me in Compile $\endgroup$
    – RunnyKine
    Jul 6, 2014 at 23:00
  • $\begingroup$ @RunnyKine hmm, maybe WolframLibrary.h has to be put in some include folder that the compiler can find. You can find WolframLibrary.h here: FileNameJoin[{$InstallationDirectory, "SystemFiles", "IncludeFiles", "C", "WolframLibrary.h"}]. Alternatively, you can convert the cCodeString into a file and put WolframLibrary.h in the same folder. $\endgroup$ Jul 7, 2014 at 5:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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