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
>
>