# RootSearch problem when called from within ParallelDo

I have been using Ted Ersek's RootSearch package for Wolfram Mathematica for many years. The package worked great, including inside the Do loop. However, inside the ParallelDo loop, it issues a diagnostic message, for example:

Thread :: tdlen: Objects of unequal length in {0,0} + {0,0} + <<8>> + <<9>> cannot be combined.

At the same time, the calculation result seems to remain correct.

At the moment I am not ready to post minimal working example here as my code is quite huge (Currently its size is 1k). But I noted that Ted Ersek is active at this forum so I hope that he let me know how I could send him my code.

Here is an MWE derived from my code. I am guessing that the problem has something to do with the fact that the equation for RootSearch uses the ParametricNDSolve solution.

Needs["ErsekRootSearch"];

debug = 0; (* controls how extensive is debug printout *)

Λs = {1.1, 2, 5, 10, 50, 100, 500}; (* list of values for ParallelDo cycle *)

results = {{"k", "μ", "ν", "K", "Λ", "βcrit1", "βcrit2", "$KernelID"}}; SetSharedVariable[results]; DistributeDefinitions[RootSearch]; Km = 8; k = ∞; μ = 1; ν = 2; (* particular parameters from outer Do cycle not shown here *) "Boundary condition at z=1: BC1[β_,K_,Λ_]" (* used in RootSearch *) BC1[β_, K_, Λ_] = -(((-1 + K) (3 K^2 - 2 β) β ϕ[1])/( 2 (K^2 - β) (K^2 - β + K^2 Λ))) "Equation to be solved: eqn=" eqn = {1/256 (49 β (1 - ( 64 (8 - 7 z)^2 β (128 - (8 - 7 z)^2 β))/(512 - 8 (8 - 7 z)^2 β)^2) ϕ[z] + 1/(512 - 8 (8 - 7 z)^2 β)^2 6272 β (32 (-1 + z) z^3 β^2 - z^4 β^2 + 64 z^2 β (1 - 6 (-1 + z)^2 β) + 1024 (-1 + z) z β (-1 + 2 (-1 + z)^2 β) + 4096 (-2 + (-1 + z)^2 β - (-1 + z)^4 β^2)) ϕ[z] + 56 (8 - 7 z) β Derivative[1][ϕ][z] + 4 (64 - (8 - 7 z)^2 β + 64 Λ) (ϕ′′[z]) == 0, ϕ[0] == 1, Derivative[1][ϕ][0] == (7 β (-3 + 2 β) ϕ[0])/( 16 (-1 + β) (-1 + β - Λ))} "Solution of ParametricNDSolve:" sol1p = ParametricNDSolve[ eqn, ϕ, {z, 0, 1}, {{β, 0, 1}, {Λ, 1, 10000}}] ϕ1p[z_, β_, Λ_] = ϕ[β, Λ][z] /. sol1p ϕ1p1[z_, β_, Λ_] = D[ϕ1p[z, β, Λ], z] (*===================================================================*) DistributeDefinitions[ϕ1p, ϕ1p1, BC1]; ParallelDo[ ClearAll[βcrit1, βcrit2, Err1, Err2]; If[debug > 2 , Print@ Plot[{ϕ1p1[1, β, Λm], BC1[β, Km, Λm] /. {ϕ[1] -> ϕ1p[ 1, β, Λm]}}, {β, 0, 0.9999} , AxesLabel -> {"β", {ϕ', BC1}} , PlotLegends -> {ϕ', BC1} , PlotLabel -> StringForm[ "k=, μ=, ν=, K=, Λ=", k, μ, ν, Km, Λm]]; ]; βcrit1[$KernelID] =
Quiet@FindRoot[ϕ1p1[1, β, Λm] ==
BC1[β,
Km, Λm] /. {ϕ[1] -> ϕ1p[
1, β, Λm]}, {β, 0.999, 0.1, 0.9999}
];

Err1[$$KernelID] = (ϕ1p1[1, β, Λm] - BC1[β, Km, Λm]) /. {ϕ[1] -> ϕ1p[ 1, β, Λm]} /. βcrit1[$$KernelID];
If[debug > 0,
Print@ToString@
StringForm[
"FindRoot1:: k=, μ=, ν=, K=, \
Λ=, β1=, Err1=", k, μ, ν,
Km, Λm, Quiet@βcrit1[$$KernelID][[1, 2]], Err1[$$KernelID]]];

Print@StringForm[
"6, kernel /, Err1=", $$KernelID,$$KernelCount,
Err1[$KernelID]]; βcrit2[$$KernelID] = RootSearch[ϕ1p1[1, β, Λm] == BC1[β, Km, Λm] /. {ϕ[1] -> ϕ1p[ 1, β, Λm]}, {β, 0.1, 0.9999} , InitialSamples -> 10 , RootTest -> (Abs[#2] < 10^-4 &) ]; Print@StringForm[ "7, kernel /, βcrit2=",$$KernelID, $$KernelCount, \ βcrit2[$$KernelID]]; If[Length[βcrit2[$$KernelID]] > 0 , Err2[$$KernelID] = (ϕ1p1[1, β, Λm] - BC1[β, Km, Λm]) /. {ϕ[ 1] -> ϕ1p[1, β, Λm]} /. First[βcrit2[$$KernelID]] , Err2[$$KernelID] =.; ]; Print@StringForm[ "8, kernel /, Err2=",$KernelID, $KernelCount, Err2[$KernelID]];
If[debug > 0,
Print@ToString@
StringForm[
"RootSearch2:: k=, μ=, ν=, β2=, Err2=",
k, μ, ν, Quiet@βcrit2[$$KernelID][[1, 1, 2]], Err2[$$KernelID]
]
];

If[Not@NumericQ[Err1[$$KernelID]] || Abs[Err1[$$KernelID]] >
10^-4, βcrit1[$$KernelID] = {β -> β1}]; If[Not@NumericQ[Err2[$$KernelID]] ||
Abs[Err2[$$KernelID]] > 10^-4, βcrit2[$$KernelID] = {{β -> β2}}];
AppendTo[results,
elem = {k, μ, ν,
Km, Λm, βcrit1[$$KernelID][[1, 2]], βcrit2[$$KernelID][[1, 1, 2]], $KernelID}]; Print["elem=", elem]; , {Λm, Λs} ]; (*===================================================================*) results // TableForm results  EDIT #2: I have located a chunck of code in ResearchRoot package that yields the warning message. It is located inside a defintion of RootSearch function. In the code below the source of warning is located between two Print commands added by me: Print["myRootSearch 5: x=", x]; If[invalidQ,$Failed,
(*else*)f=Evaluate[{SamplePoints[#],lhs-rhs/.x->#}]&;
f=ReplacePart[f,Set,{1,0}];
Print["myRootSearch 5.1: f=",f];
u=Function[Evaluate[
withTemp[{setTemp[y,(lhs-rhs)/.x->#],
setTemp[yyp,((lhs-rhs)/D[lhs-rhs,x])/.x->#]},
Hold[
SamplePoints[#]=y,
Which[TrueQ[y==0],0,NumericQ[yyp],yyp
,True,Indeterminate]
]]]];
Print["myRootSearch, 5.2"];
u=ReplacePart[
u/.{withTemp->With,setTemp->Set},{CompoundExpression},{{1,
2,0}},{{1}}];


These Prints produce the following output (from Kernel 2):

myRootSearch 5: x=[Beta]

myRootSearch 5.1: f=(KiamyRootSearchPrivateSamplePoints[#1]=(ParametricFunction[Expression: [Phi] Parameters: {#1,500}

][#1,500]^[Prime])[1]+(7 (192-2 #1) #1 ParametricFunction[Expression: [Phi] Parameters: {#1,500}

][#1,500][1])/(2 (64-#1) (32064-#1))

Thread::tdlen : Objects of unequal length in {0,0}+{0,0}+<<8>>+<<9>> cannot be combined.

myRootSearch, 5.2

As you see, the warning goes between print points 5.1 and 5.2.

EDIT #3:

I have found that error is caused by D[lhs-rhs,x]`.

• If the same exact code works with Do but fails with ParallelDo, then is it possible that the value of some global variable in your code has not been shared to the parallel kernels? As you can see, though, without a minimal working example reproducing your issue, we can only take a wild guess. Would you be able to reproduce the problem with a much smaller test function? Jan 8, 2022 at 17:48
• @MarcoB Indeed, for other more simple equation RootSearch works seamlessly. In my case my equation envolves a function composed of solution of ParametricNDSolve. I will try futher investigate the case. Jan 9, 2022 at 1:56
• @MarcoB: I have added MWE that reproduces the problem. Jan 9, 2022 at 4:26
• @user:460: Ted, please look at my post mathematica.stackexchange.com/questions/261792/… Jan 9, 2022 at 4:54