В настоящее время я работаю с некоторым кодом Mathematica для выполнения итерации Пикара.Сам код работает нормально, но я пытаюсь сделать его более эффективным.У меня был некоторый успех, но я ищу предложения.Возможно, уже не удастся ускорить его, но у меня кончились идеи, и я надеюсь, что люди с большим опытом программирования / Mathematica, чем я, могли бы сделать некоторые предложения.Я публикую только саму итерацию, но могу предоставить дополнительную информацию по мере необходимости.
Приведенный ниже код был отредактирован как полностью исполняемый в соответствии с запросом
Также я изменил его сЦикл while to Do для упрощения тестирования, поскольку конвергенция не требуется.
Clear["Global`*"]
ngrid = 2048;
delr = 4/100;
delk = \[Pi]/delr/ngrid;
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}];
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}];
wa[x_] := (19 + .5 x) Exp[-.7 x] + 1
wb[x_] := (19 + .1 x) Exp[-.2 x] + 1
wd = SetPrecision[
Table[{{wa[(i - 1/2) delk], 0}, {0, wb[(i - 1/2) delk]}}, {i, 1,
ngrid}], 26];
sigmaAA = 1;
hcloseAA = {};
i = 1;
While[(i - 1/2)*delr < sigmaAA, hcloseAA = Append[hcloseAA, -1]; i++]
hcloselenAA = Length[hcloseAA];
hcloseAB = hcloseAA;
hcloselenAB = hcloselenAA;
hcloseBB = hcloseAA;
hcloselenBB = hcloselenAA;
ccloseAA = {};
i = ngrid;
While[(i - 1/2)*delr >= sigmaAA, ccloseAA = Append[ccloseAA, 0]; i--]
ccloselenAA = Length[ccloseAA];
ccloselenAA = Length[ccloseAA];
ccloseAB = ccloseAA;
ccloselenAB = ccloselenAA;
ccloseBB = ccloseAA;
ccloselenBB = ccloselenAA;
na = 20;
nb = 20;
pa = 27/(1000 \[Pi]);
pb = 27/(1000 \[Pi]);
p = {{na pa, 0}, {0, nb pb}};
id = {{1, 0}, {0, 1}};
AFD = 1;
AFDList = {};
timelist = {};
gammainitial = Table[{{0, 0}, {0, 0}}, {ngrid}];
gammafirst = gammainitial;
step = 1;
tol = 10^-7;
old = 95/100;
new = 1 - old;
Do[
t = AbsoluteTime[];
extractgAA = Table[Extract[gammafirst, {i, 1, 1}], {i, hcloselenAA}];
extractgBB = Table[Extract[gammafirst, {i, 2, 2}], {i, hcloselenBB}];
extractgAB = Table[Extract[gammafirst, {i, 1, 2}], {i, hcloselenAB}];
csolutionAA = (Join[hcloseAA - extractgAA, ccloseAA]) rvalues;
csolutionBB = (Join[hcloseBB - extractgBB, ccloseBB]) rvalues;
csolutionAB = (Join[hcloseAB - extractgAB, ccloseAB]) rvalues;
chatAA = FourierDST[SetPrecision[csolutionAA, 32], 4];
chatBB = FourierDST[SetPrecision[csolutionBB, 32], 4];
chatAB = FourierDST[SetPrecision[csolutionAB, 32], 4];
chatmatrix =
2 \[Pi] delr Sqrt[2*ngrid]*
Transpose[{Transpose[{chatAA, chatAB}],
Transpose[{chatAB, chatBB}]}]/kvalues;
gammahat =
Table[(wd[[i]].chatmatrix[[i]].(Inverse[
id - p.wd[[i]].chatmatrix[[i]]]).wd[[i]] -
chatmatrix[[i]]) kvalues[[i]], {i, ngrid}];
gammaAA =
FourierDST[SetPrecision[Table[gammahat[[i, 1, 1]], {i, ngrid}], 32],
4];
gammaBB =
FourierDST[SetPrecision[Table[gammahat[[i, 2, 2]], {i, ngrid}], 32],
4];
gammaAB =
FourierDST[SetPrecision[Table[gammahat[[i, 1, 2]], {i, ngrid}], 32],
4];
gammasecond =
Transpose[{Transpose[{gammaAA, gammaAB}],
Transpose[{gammaAB, gammaBB}]}]/(rvalues 2 \[Pi] delr Sqrt[
2*ngrid]);
AFD = Sqrt[
1/ngrid Sum[((gammafirst[[i, 1, 1]] -
gammasecond[[i, 1, 1]])/(gammafirst[[i, 1, 1]] +
gammasecond[[i, 1, 1]]))^2 + ((gammafirst[[i, 2, 2]] -
gammasecond[[i, 2, 2]])/(gammafirst[[i, 2, 2]] +
gammasecond[[i, 2, 2]]))^2 + ((gammafirst[[i, 1, 2]] -
gammasecond[[i, 1, 2]])/(gammafirst[[i, 1, 2]] +
gammasecond[[i, 1, 2]]))^2 + ((gammafirst[[i, 2, 1]] -
gammasecond[[i, 2, 1]])/(gammafirst[[i, 2, 1]] +
gammasecond[[i, 2, 1]]))^2, {i, 1, ngrid}]];
gammafirst = old gammafirst + new gammasecond;
time2 = AbsoluteTime[] - t;
timelist = Append[timelist, time2], {1}]
Print["Mean time per calculation = ", Mean[timelist]]
Print["STD time per calculation = ", StandardDeviation[timelist]]
Несколько замечаний о вещах
ngrid, delr, delk, rvalues, kvalues - это только те значения, которые используются для решения проблемы.Обычно это
ngrid = 2048;
delr = 4/100;
delk = \[Pi]/delr/ngrid;
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}];
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}];
Все используемые матрицы имеют размер 2 x 2 с одинаковыми недиагоналями
Матрица тождества и матрица P (это фактически для плотности)
p = {{na pa, 0}, {0, nb pb}};
id = {{1, 0}, {0, 1}};
Основными медленными точками в вычислении, которое я идентифицировал, являются вычисления FourierDST (прямое и обратное преобразования составляют около 40% времени вычисления). Расчет гаммахата составляет 40% времени прив оставшееся время преобладает вычисление AFD.) На моем процессоре i7 среднее время расчета за цикл составляет 1,52 секунды.Я надеюсь получить его менее чем за секунду, но это может оказаться невозможным.Я надеялся ввести несколько параллельных вычислений, которые были опробованы как с командами ParallelTable, так и с использованием ParallelSubmit WaitAll.Тем не менее, я обнаружил, что любое ускорение от параллельных вычислений было компенсировано временем связи от главного ядра к другим ядрам (по крайней мере, я так предполагаю, поскольку вычисления для новых данных занимают вдвое больше времени, чем просто пересчет существующих данных).Я предположил, что это означало, что замедление было в распространении новых списков) Я играл с DistributDefinitions, а также SetSharedVariable, однако не смог заставить это что-либо сделать.
Одна вещь, которая меня интересуетесли использование Table для моих дискретных вычислений - лучший способ сделать это?
Я также подумал, что мог бы переписать это таким образом, чтобы иметь возможность скомпилировать его, но я понимаю, что это сработает, только если вы работаете с точностью станка, где мне нужно работать с большей точностьючтобы получить сближение.
Заранее спасибо за любые предложения.