Comment paralléliser l'intégration dans Mathematica 8

Question

Quelqu'un a une idée de l'utilisation de tous les cœurs pour le calcul de l'intégration? J'ai besoin d'utiliser la paralléliseur ou la table parallèle mais comment?

f[r_] := Sum[(((-1)^n*(2*r - 2*n - 7)!!)/(2^n*n!*(r - 2*n - 1)!))* 
x^(r - 2*n - 1), {n, 0, r/2}]; 


Nw := Transpose[Table[f[j], {i, 1}, {j, 5, 200, 1}]]; 

X1 = Integrate[Nw . Transpose[Nw], {x, -1, 1}]; 

Y1 = Integrate[D[Nw, {x, 2}] . Transpose[D[Nw, {x, 2}]], {x, -1, 1}]; 

X1//MatrixForm 
Y1//MatrixForm 

Answer 1

Si on aide à intégrer un peu en développant les éléments de la matrice d'abord, choses sont réalisables avec un peu d'effort. Sur un ordinateur portable quad-core avec Windows et Mathematica 8.0.4, le code ci-dessous s'exécute pour le demandé DIM = 200 en environ 13 minutes, pour DIM = 50 le code s'exécute en 6 secondes.

---

$starttime = AbsoluteTime[]; Quiet[LaunchKernels[]]; 
DIM = 200; 
Print["$Version = ", $Version, " ||| ", "Number of Kernels : ", Length[Kernels[]]]; 
f[r_] := f[r] = Sum[(((-1)^n*(-(2*n) + 2*r - 7)!!)*x^(-(2*n) + r - 1))/(2^n*n!*(-(2*n) + r - 1)!), {n, 0, r/2}]; 
Nw = Transpose[Table[f[j], {i, 1}, {j, 5, DIM, 1}]]; 
nw2 = Nw . Transpose[Nw]; 
Print["Seconds for expanding Nw.Transpose[Nm] ", Round[First[AbsoluteTiming[nw3 = ParallelMap[Expand, nw2]; ]]]]; 
Print["do the integral once: ", Integrate[x^n, {x, -1, 1}, Assumptions -> n > -1]]; 
Print["the integration can be written as a simple rule: ", intrule = (pol_Plus)?(PolynomialQ[#1, x] &) :> 
    (Select[pol, !FreeQ[#1, x] & ] /. x^(n_.) /; n > -1 :> ((-1)^n + 1)/(n + 1)) + 2*(pol /. x -> 0)]; 
Print["Seconds for integrating Nw.Transpose[Nw] : ", Round[First[AbsoluteTiming[X1 = ParallelTable[row /. intrule, {row, nw3}]; ]]]]; 
Print["expanding: ", Round[First[AbsoluteTiming[preY1 = ParallelMap[Expand, D[Nw, {x, 2}] . Transpose[D[Nw, {x, 2}]]]; ]]]]; 
Print["Seconds for integrating : ", Round[First[AbsoluteTiming[Y1 = ParallelTable[py /. intrule, {py, preY1}]; ]]]]; 
Print["X1 = ", (Shallow[#1, {4, 4}] &)[X1]]; 
Print["Y1 = ", (Shallow[#1, {4, 4}] &)[Y1]]; 
Print["seq Y1 : ", Simplify[FindSequenceFunction[Diagonal[Y1], n]]]; 
Print["seq X1 0 : ",Simplify[FindSequenceFunction[Diagonal[X1, 0], n]]]; 
Print["seq X1 2: ",Simplify[FindSequenceFunction[Diagonal[X1, 2], n]]]; 
Print["seq X1 4: ",Simplify[FindSequenceFunction[Diagonal[X1, 4], n]]]; 
Print["overall time needed in seconds: ", Round[AbsoluteTime[] - $starttime]]; 

Answer 2

j'ai changé l'intégration d'une liste dans une liste des intégrations afin que je puisse utiliser ParallelTable:

X1par=ParallelTable[Integrate[i, {x, -1, 1}], {i, Nw.Transpose[Nw]}]; 

X1par==X1 

(* ===> True *) 

Y1par = ParallelTable[Integrate[i,{x,-1,1}],{i,D[Nw,{x,2}].Transpose[D[Nw,{x,2}]]}] 

Y1 == Y1par 

(* ===> True *) 

Dans mes horaires, avec {j, 5, 30, 1} au lieu de {j, 5, 200, 1} pour limiter le temps utilisé un peu, ce est environ 3,4 fois plus rapide sur mon quod-core. Mais il peut faire encore plus vite avec:

X2par = Parallelize[Integrate[#, {x, -1, 1}] & /@ (Nw.Transpose[Nw])] 

X2par == X1par == X1 

(* ===> True *) 

Ceci est environ 6,8 fois plus rapide, un facteur de 2,3 dont est due à Parallelize.

Timing et AbsoluteTiming ne sont pas très fiables en cas d'exécution parallèle. J'ai utilisé AbsoluteTime avant et après chaque ligne et a pris la différence.

---

EDIT

Nous ne devons pas oublier ParallelMap:

Au niveau de la liste grossière (1):

ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {1}] 

Au plus profond de la liste (la plus fine parallélisation grainée):

ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {2}]