MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: multiple outputs from compiled function

  • To: mathgroup at smc.vnet.net
  • Subject: [mg114387] Re: multiple outputs from compiled function
  • From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
  • Date: Thu, 2 Dec 2010 05:43:51 -0500 (EST)

Hi Eric,

On Tue, 30 Nov 2010, Eric Michielssen wrote:

> Thanks for the suggestion.

most welcome.

>
> I was thinking more about a routine that populates a huge sparse matrix, say
> for finite element analysis. SparseArray does not work within Compile, and
> that's understandable.
>
> So the next best thing would be to have a compiled function generate a long
> list of positions ( pos = {pos1,pos2,pos3...} with posi = {n1,n2} ) and a
> list of values ( val = {val1,val2,val3,...} ), and then construct the
> SparseArray externally to a compiled function.
>
> But how do I get my compiled function to return both pos and val? Upon

You don not. Let me show you why.

> reaching Return[{pos,val}], Mathematica reverts to the noncompiled function,
> in essence starting over from scratch. (I can only think of a few not so
> elegant solutions, none of which would apply in situations more complex than
> the above).
>

First you might find this link interesting:

http://library.wolfram.com/infocenter/Conferences/7549/

(*Some setup:*)

$HistoryLength = 0;
On["Packing"]
pack = Developer`ToPackedArray;
packedQ = Developer`PackedArrayQ;

(*OK, lets create an example mesh. For this we use the new TetGenLink:*)

Needs["TetGenLink`"]
pts = {{0., 0., 0.}, {10., 0., 0.}, {10., 2., 0.}, {0., 2., 0.}, {0.,
     0., 1.}, {10., 0., 1.}, {10., 2., 1.}, {0., 2., 1.}};
facets = {{{1, 2, 3, 4}}, {{5, 6, 7, 8}}, {{1, 2, 6, 5}}, {{2, 3, 7,
      6}}, {{3, 4, 8, 7}}, {{4, 1, 5, 8}}};
inInst = TetGenCreate[];
TetGenSetPoints[inInst, pts];
TetGenSetFacets[inInst, facets];
(* to get more elements use: "pqa0.000075" *)
outInst = TetGenTetrahedralize[inInst, "pqa0.075"];
coords = TetGenGetPoints[outInst];
surface = TetGenGetFaces[outInst];
meshElements = TetGenGetElements[outInst];
dofs = Length[coords]
Length[meshElements]

(* visualize *)
TetrahedraWireframe[i_] :=
  Line[ Join @@ (i[[All, #]] & /@ {{1, 2}, {2, 3}, {3, 1}, {1, 4}, {2,
        4}, {3, 4}})]
Graphics3D[GraphicsComplex[coords, TetrahedraWireframe[meshElements]],
   Boxed -> False]


(* find position of Dirichlet values *)
diriPos1 =
   Flatten[ pack[
     Position[
      coords[[All, 1]], _?(Function[xCoord, xCoord <= 0.]), {1},
      Heads -> False ]]];
diriPos2 =
   Flatten[ pack[
     Position[
      coords[[All, 1]], _?(Function[xCoord, xCoord >= 10.]), {1},
      Heads -> False ]]];
diriPos = Join[ diriPos1, diriPos2];
diriVals =
   Join[ ConstantArray[ 1., {Length[ diriPos1]}],
    ConstantArray[ 0., {Length[ diriPos2]}]];

(* visualize *)
Show[ {
   Graphics3D[GraphicsComplex[coords, Polygon[surface]]],
   Graphics3D[ {
     {Red, PointSize[0.02], Point[ coords[[diriPos1]]]},
     {Blue, PointSize[0.02], Point[ coords[[diriPos2]]]}}]
   }, Boxed -> False]


(* symbolically precompute a FEM diffusion kernel and inject in compile *)
tmp = Join[ {{1, 1, 1, 1}},
    Transpose[ Quiet[Array[Part[var, ##] &, {4, 3}]]]];
me = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
p = Inverse[ tmp].me;
help = Transpose[ (p.Transpose[p])*Abs[Det[ tmp]]/6 ];
femk = With[{code = help},
    Compile[{{coords, _Real, 2}, {incidents, _Integer, 1}}, Block[{var},
      var = coords[[incidents]];
      code
      ], RuntimeAttributes -> Listable, CompilationTarget -> "C"]];


(* compute the local elements. - note just the elements. we do the 
positions later. For fun, you may want to scale the number of elements 
later *)
AbsoluteTiming[ localElements = femk[coords, meshElements]; ]
packedQ[localElements]

(* this is the function to assemble the local elements into a sparse array 
*)
matrixAssembly[ values_, pos_, dim_] := Block[{matrix},
   System`SetSystemOptions[
    "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
   matrix = SparseArray[ pos -> Flatten[ values], dim];
   System`SetSystemOptions[
    "SparseArrayOptions" -> {"TreatRepeatedEntries" -> 0}];
   Return[ matrix]]


(*
"TreatRepeatedEntries" does the following. If during the SparseArray 
creation duplicate positions are found such {1,1}->1 and {1,1}->2, the 
values are added and you'll find the content to be {1,1}->3 *)

(* now the positions:
The reason to compute the position independent from the elements is to 
save them for later. Consider this, if you want to assemble a mass matrix 
the entries will be at exactly the same positions as in the stiffness 
matrix, just different values. This is quite efficient.
  *)
AbsoluteTiming[
  pos = Flatten[ Map[ Outer[ List, #, #] &, meshElements], 2]; ]


(*assemble the system*)
AbsoluteTiming[
  stiffness = matrixAssembly[ localElements, pos, dofs ] ]

(*now for the Dirichlet conditions, we need something to pass the 
stiffness by reference and modify that matrix directly, you can do so*)

SetAttributes[dirichletBoundary, HoldFirst]
dirichletBoundary[ {load_, matrix_}, fPos_List, fValue_List] :=
  Block[{},
   load -= matrix[[All, fPos]].fValue;
   load[[fPos]] = fValue;
   matrix[[All, fPos]] = matrix[[fPos, All]] = 0.;
   Transpose[ {fPos, fPos}] -> ConstantArray[ 1., {Length[ fPos]}];
   matrix += SparseArray[
     Transpose[ {fPos, fPos}] -> ConstantArray[ 1., {Length[ fPos]}],
     Dimensions[ matrix], 0];
   ]

load = Table[ 0., {dofs}];
dirichletBoundary[{load, stiffness}, diriPos, diriVals]

(* dirichletBoundary returns Null and operates on stiffness and load 
directly *)

(* solve with a direct solver *)
AbsoluteTiming[
  solution = LinearSolve[ stiffness, load]; ]

(* visualize *)
Graphics3D[{EdgeForm[],
   GraphicsComplex[ coords, Polygon[ surface ],
    VertexColors -> (ColorData["TemperatureMap"][#] & /@ solution) ]},
  Boxed -> False, Background -> Gray ]



I hope this helps to get you FEM code going.

Oliver



> Eric
>
> -----Original Message-----
> From: Oliver Ruebenkoenig [mailto:ruebenko at wolfram.com]
> Sent: Tuesday, November 30, 2010 6:24 AM
> To: mathgroup at smc.vnet.net
> Subject: [mg114315] Re: multiple outputs from compiled function
>
> On Tue, 30 Nov 2010, Eric Michielssen wrote:
>
>> I have a compiled function that internally generates 10 large lists (list1
>> though list10), all of different ranks and sizes. Problem is that a
> compiled
>> function can only return one packed array, so returning them to the main
>> code is an issue. Of course, I could combine/pack all these lists into one
>> list by Flatten[list1,...,list10], return the new structure, and then
> unpack
>> it in the main code. Is there a less cumbersome way of doing this?
>>
>> Eric Michielssen
>>
>>
>>
>>
> Eric,
>
> Here are some thoughts:
>
> It is correct that compile can only work with rectangular tensors. See the
> following:
>
> f = Compile[{{t1, _Real, 2}, {t2, _Real, 2}}, Join[t1, t2]]
>
> f[RandomReal[{0, 1}, {2, 2}], RandomReal[{0, 1}, {2, 1}]]
>
> compared to:
>
> f[RandomReal[{0, 1}, {2, 2}], RandomReal[{0, 1}, {2, 2}]]
>
> But if the result of the base CompiledFunction is a list, and the list
> results come out different lengths, it returns a list of PackedArrays,
>
> cf = Compile[{{x, _Real}, {n, _Integer}},
>    Module[{y = 1.}, Table[y = (y + x/y)/2, {n}]],
>    RuntimeAttributes -> Listable];
>
> res = cf[{2., 3.}, {3, 2}]
>
> {Developer`PackedArrayQ[res], Map[Developer`PackedArrayQ, res]}
>
> resp = cf[{2., 3.}, 3]
>
> Developer`PackedArrayQ[resp]
>
> What do exactly in your case depends on the code.
>
> Hope this helps,
> Oliver
>
>


  • Prev by Date: Re: How to short-circuit match failure?
  • Next by Date: Re: Mathematica 8 ImageKeypoints error
  • Previous by thread: Re: multiple outputs from compiled function
  • Next by thread: Re: multiple outputs from compiled function