3
$\begingroup$

I'd like to do this: solve some linear systems, compiling the task of getting the building up the matrices from parameter values, and then I have a compiled implementation of LinearSolve. The problem is that there are a huge number of different LHS matrices, and it take time to build them up. I only want to build one per function call. So this is the idea:

RHSarray = {{Indexed[arg1, 1], 5.}, {Indexed[arg1, 1], 
   Indexed[arg1, 2]}, {7., 9.}}
LHSarray = {{{Indexed[arg1, 1], Indexed[arg1, 2]}, {3., 
    4.}}, {{Indexed[arg1, 1], 0.}, {0, 1}}, {{1., 0.}, {0, 
    Indexed[arg1, 1]}}}
cFunStraightforward = Compile[{{i, _Integer}, {arg1, _Real, 1}},
   Module[{LHStosolve = {{}}, RHStosolve = {}},
    LHStosolve = LHSarray[[i]]; RHStosolve = RHSarray[[i]];
    LinearSolve[LHStosolve, RHStosolve]],
   CompilationOptions -> {"InlineExternalDefinitions" -> True}];

But it creates the entire rank 3 tensor LHSarray, and then selects the ith matrix.

Do[ToExpression[
  "LHSarray" <> ToString[i] <> "= LHSarray[[" <> ToString[i] <> 
   "]];"]; ToExpression[
  "RHSarray" <> ToString[i] <> "= RHSarray[[" <> ToString[i] <> 
   "]];"];, {i, 3}]
cFun = Compile[{{i, _Integer}, {arg1, _Real, 1}},
   Module[{LHStosolve = {{}}, RHStosolve = {}},
    If[i == 1, LHStosolve = LHSarray1; RHStosolve = RHSarray1];
    If[i == 2, LHStosolve = LHSarray2; RHStosolve = RHSarray2];
    If[i == 3, LHStosolve = LHSarray3; RHStosolve = RHSarray3];
    LinearSolve[LHStosolve, RHStosolve]],
   CompilationOptions -> {"InlineExternalDefinitions" -> True}];

This does what I want. Can you help me do this programmatically? And I would love to get rid of the call to ToExpression too, I just don't know how best to do that.

Edit: I found this, which works. I know there has to be a better way, but it does work.

nestIfs[expr_, n_] := 
 "If[i\[Equal]" <> ToString[n] <> ",LHStosolve=LHSarray" <> 
  ToString[n] <> ";RHStosolve=RHSarray" <> ToString[n] <> "," <> 
  expr <> "]"
ToExpression[
 "cFun = Compile[{{i,_Integer},{arg1,_Real,1}},\[IndentingNewLine]\
Module[{LHStosolve = {{}},RHStosolve = {}}," <> 
  Fold[nestIfs, "0", Reverse[Range[3]]] <> 
  ";LinearSolve[LHStosolve,RHStosolve]],\[IndentingNewLine]\
CompilationOptions\[Rule]{\"InlineExternalDefinitions\"\[Rule]True}]"]
cFun[1, {1., 2.}]
$\endgroup$
  • $\begingroup$ Well, are you aware that LinearSolve is not compilable so you'll gain little speed-up by compiling your function? $\endgroup$ – xzczd Jul 14 '20 at 3:28
  • $\begingroup$ @xzczd I have a compiled version of LinearSolve borrowed from here, which is somewhat faster. And I get a good speed up from all of the complex building up of each matrix (i.e. all the manipulations on what I call arg1.) I could really use some help writing code that isn't ToExpression[string] though! $\endgroup$ – NathanRL Jul 14 '20 at 3:44
  • $\begingroup$ 1. I think using ToExpression in poster's definitions is just fine. 2. Maybe a better fitting question title is "Code generation for Compile"? (I am not criticizing, by the way, I am asking in order to clarify poster's intent.) $\endgroup$ – Anton Antonov Jul 14 '20 at 12:03
2
$\begingroup$

How about having n separate compiled functions?:

cfunclst = MapThread[Compile[{{arg1, _Real, 1}}, LinearSolve@##] &, {LHSarray, RHSarray}]

Then just use e.g. cFunclst[[1]][{1., 2.}].

If you insist on using a single compiled function, then a possible solution is

toseq = Flatten[#, 1] &@Transpose@{Range@Length@#, #} &;

cfunc = Hold@
     Compile[{{i, _Integer}, {arg1, _Real, 1}}, 
      LinearSolve[sw1, sw2]] /. {sw1 -> switch[i, ##, _, {{0.}}] & @@ toseq@LHSarray, 
     sw2 -> switch[i, ##, _, {0.}] & @@ toseq@RHSarray} /. switch -> Switch // 
  ReleaseHold
$\endgroup$
  • $\begingroup$ Yeah, that's how I have it now. It seemed better to have it all in one function, for distributing it to the parallel workers and loading and unloading, but what do I know. $\endgroup$ – NathanRL Jul 14 '20 at 5:17
  • $\begingroup$ @NathanRL Check my update. $\endgroup$ – xzczd Jul 14 '20 at 6:35

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.