У кого-нибудь есть идея, как использовать все ядра для расчета интеграции? Мне нужно использовать распараллеливание или параллельную таблицу, но как?
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
Если кто-то немного поможет интегрировать, сначала расширив элементы матрицы, то все можно сделать, приложив немного усилий.
На четырехъядерном ноутбуке с Windows и Mathematica 8.0.4 приведенный ниже код выполняется для запрошенного DIM=200 примерно за 13 минут, для DIM=50 код выполняется за 6 секунд.
$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]];
Я изменил интеграцию списка в список интеграций, чтобы можно было использовать 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 *)
В моих таймингах с {j, 5, 30, 1} вместо {j, 5, 200, 1} для ограничения используемого времени это примерно в 3,4 раза быстрее на моем четырехъядерном процессоре. Но это можно сделать еще быстрее, если:
X2par = Parallelize[Integrate[#, {x, -1, 1}] & /@ (Nw.Transpose[Nw])]
X2par == X1par == X1
(* ===> True *)
Это примерно в 6,8 раза быстрее, в 2,3 раза больше из-за Parallelize.
Timing и AbsoluteTiming не очень заслуживают доверия, когда речь идет о параллельном выполнении. Я использовал AbsoluteTime до и после каждой строки и взял разницу.
ИЗМЕНИТЬ
Мы не должны забывать ParallelMap:
На самом грубом уровне списка (1):
ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {1}]
На самом глубоком уровне списка (наиболее тонкое распараллеливание):
ParallelMap[Integrate[#, {x, -1, 1}] &, Nw.Transpose[Nw], {2}]
