That's right. fsolve is very likely not thread safe. Most of the complicated system Library (interpreted) commands like solve, limit, fsolve, dsolve, int, etc, are not. (I mean Maple 14, and it's true for Maple 12 as well.)

acer

Extrapolating is somewhat dubious, especially back so far. But are you trying to do something like this?

numer_extrap.mws

acer

How about this, where `ee` is your big expression with the arctan call in it.

stuff := seq([op(t)],t in indets(ee,specfunc(anything,arctan)));

stuff[1];
stuff[2];

acer

You wrote "at least 11K nonzero entries". That's not very large, as a dense 110x110 Matrix has about that much.

I'd want to hear what you'd answer to my earlier query about how you plan to judge acceptable roundoff (ie. according to whether it is unlikely, or...?).

But the following code shows that 5000 individual sums over all entries of dense 110x110 float[8] Arrays can be done on a fast PC in less than a second.

Do you have some special reason to think that roundoff is going to be a problem for your data?

The procedure `P` below is a bit like what's under the hood of the 2D AddAlongDimension routine (it calls the same external BLAS function), except it operates on a single Vector. And its first argument must be a float[8] Vector. The procedure `f2` simply aliases a 2D Array to a 1D Vector and then calls P on that.

P := proc(V::Vector(datatype=float[8]),n::posint)
local extlib, ExtCall, x;
   extlib := ExternalCalling:-ExternalLibraryName("linalg", 'HWFloat');
   ExtCall := ExternalCalling:-DefineExternal('hw_f06ecf', extlib);
   x:=Vector(1,'datatype'='float[8]');
   ExtCall(n,1.0,V,0,1,x,0,0);
   x[1];
end proc:

f2:=proc(A::Array(datatype=float[8]))
local B, m, n;
  #(m,n):=op(map(rhs,[rtable_dims(A)]));
  m:=rhs(rtable_dims(A)[1]);
  n:=rhs(rtable_dims(A)[2]);
  B:=ArrayTools:-Alias(A,[m*n]);
  P(B, m*n);
end proc:

M:=LinearAlgebra:-RandomMatrix(110,110,generator=-1.0..1.0,
                               outputoptions=[datatype=float[8]]):
rtable_options(M,subtype=Array);

CodeTools:-Usage( f2(M) );
memory used=64.82KiB, alloc change=0 bytes, cpu time=0.00s, real time=0.02s

                    -11.847884729671568

CodeTools:-Usage( add(t, t in M) );
memory used=0.60MiB, alloc change=0.62MiB, cpu time=0.02s, real time=0.02s

                    -11.847884729671568

str,st,ba,bu:=time[real](),time(),kernelopts(bytesalloc),kernelopts(bytesused):
for i from 1 to 5000 do
   f2(M);
end do:
time[real]()-str,time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

                0.640, 0.640, 22540256, 26030924


str,st,ba,bu:=time[real](),time(),kernelopts(bytesalloc),kernelopts(bytesused):
for i from 1 to 5000 do
   add(t, t in M);
end do:
time[real]()-str,time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;

              33.462, 33.462, 18477768, 3146484956

acer

ee:=1/(5^(1/2)-1)*(10+2*5^(1/2))^(1/2):
simplify(signum(ee)*arccos(cos(arctan(ee))));

acer

That error from `sort` occurs because BigSLTarray is not actually populated with numbers. It's populated with calls to the random number generating proc that rand(-1000..1000) returns.

Observe how the first example below produces nothing but scalar*proc, while the second example actually generates a random float value. You've accidentally coded `foo` like the first example.

> # first example, like in the posted code
> (1/198781)*rand(-1000 .. 1000);

    1          
  ------ proc()
  198781       
    (proc() option builtin = RandNumberInterface;  end proc;)(6, 2001, 11)
     - 1000
  end proc;

> evalf[50](%);

  0.0000050306618841841021023136014005362685568540252841066 (proc()
    (proc() option builtin = RandNumberInterface;  end proc;)(6, 2001, 11)
     - 1000
  end proc;)

Since BigSLTArray doesn't have actual numeric values, the `sort` routine cannot ascertain whether abs(a)<=abs(b) for `a`, `b` as pairs of entries.

> # second example
> F:=rand(-1000 .. 1000):

> evalf[50](F()/198781);
          -0.00070429266378577429432390419607507759795956353977493

You should be able to get good speed and memory performance by Compiling the routine which adds all the entries. That's probably still true even if you sort the entries (due to a need to avoid roundoff, say) provided that its done inplace. You should be able to Thread the larger process. If I understand you correctly, then at present you are just trying to discover whether the sorting is "necessary". Does your use of random values reflect something in the larger problem, or is it there merely because you are trying to discover how likely the worst roundoff case will be?

Addressing a separate question, in one of the followup responses in this thread: There is no contructor lowercase `float`. It is `Float`. And so floating-point infinity can be created with calls such as Float(infinity) or Float(1,infinity). But since there is no such constructor `float`, then the call float(1,infinity) merely produce an unevaluated function call which Maple does not recognize.

acer

restart:

t_init := 20: # [C] Start temp i mælken
h := 10: # [W/m2*K] varme transmissionskoefficient
A := 0.07: # [m2] Overfladeareal
m := 1: # [kg] massen af mælken
c := 4190: # [j/kgK] specifik varmekapacitet, mælk og karton
t_R := 5: # [C] Temperatur i køleskabet
tau_slut := 3600 * 3: # [s] Simuleringstid
dtau := 10: # [s] Tids step

# Beregning - nogle hjælpestørrelser

i_max := tau_slut/dtau: # Antal elementer (gennemløb i loop'et)

# Tre vektorer oprettes ud udfyldes med nuller

tau := Vector(i_max): # Den numeriske løsning
t := Vector(i_max): # Den numeriske løsning

# Så sættes nogle begyndelsesværdier

t[1] := t_init:
tau[1] := 0: # Starttidspunktet

#Simulering

for i from 2 to i_max do
    tau[i] := tau[i - 1] + dtau; # [s] Den nye tid
    dt := -(h * A)/(m * c) * (t[i - 1] - t_R) * dtau; # [-] Temp. ændringen
    t[i] := t[i - 1] + dt; # [-] Den nye temp.
end do:

plot(<(tau/3600)|t>,gridlines=true,
                    legend="analytisk",legendstyle=[location=right],
                    caption="Mælketemperatur");

acer

> eqn1 := u = V * cos(a) * cos(b):
> eqn2 := v = V*sin(b):
> eqn3 := w = V*sin(a)*cos(b):
> # This below is a roundabout and somewhat arbitrary way to
> # solve for sin(b). By "arbitrary" I mean that it involves
> # visually inspection of the original equations, and deduction
> # of which one(s) to use further.
> #

> res:=solve({eqn1, eqn2, eqn3}, {V,a,b}):
> map(sin, b=eval(b,res)):
> simplify(%);

                                   v             
              sin(b) = --------------------------
                             /  2    2    2    2\
                       RootOf\_Z  - v  - w  - u /

# How is the equations for V in `res` known to be the one to use?
# This approach may work less well in general.
#

> subs(eval(V,res)=V,%);

                                    v
                           sin(b) = -
                                    V

> # Is this next `solve` attempt below not as satisfactory a way
> # to solve for sin(b), because u,v,w appear "too often" in the result?
> # If these are satisfactory, we can stop right here.
> #

> solve({eqn1,eqn2,eqn3},[sin(b),cos(a),V]);

    [[         v sin(a) cos(b)           u sin(a)            w      ]]
    [[sin(b) = ---------------, cos(a) = --------, V = -------------]]
    [[                w                     w          sin(a) cos(b)]]

> sin(b)=eval(sin(b),%[1]);

                             v sin(a) cos(b)
                    sin(b) = ---------------
                                    w       

> # If the above is not satisfactory, because u,v,w appear "too often"
> # in the RHS's, then we could try to force some preference.
> #
> # Notice that the equations are a system of polynomials, if
> # the trig function calls are considered as being "variables"
> # or atomic names.
> # (This is somewhat justified since `a` and `b` appear only
> # as arguments of those trig calls and not otherwise.)
> # We can thus replace all the trig calls like sin(b) with
> # names, temporarily. Then solve. Then revert those names.
> # Hey, it's worth trying...
> #

> G:=proc(eqs::set, X, lesspreferred::list, morepreferred::set)
> local freezeeqs, thaweqs, T;
>    freezeeqs:=map(t->t=freeze(t),
>                   indets(eqs,specfunc(name,{sin,cos,tan})));
>    thaweqs:=map(rhs=lhs,freezeeqs);
>    T:=SolveTools:-PolynomialSystem(subs(freezeeqs,eqs),
>                                    subs(freezeeqs,
>                                         [X,
>                                          op(lesspreferred),
>                                          op(morepreferred minus {X})]),
>                                    'backsub'=false,'engine'='groebner');
>    subs(thaweqs,isolate(T[1][-1],eval(X,freezeeqs)));
> end proc:

> G({eqn1,eqn2,eqn3}, sin(b), [u,v,w], {sin(b),cos(a),V});

                                    v
                           sin(b) = -
                                    V

> G({eqn1,eqn2,eqn3}, cos(a), [u,v,w], {sin(b),cos(a),V});

                                   u    
                       cos(a) = --------
                                V cos(b)

> G({eqn1,eqn2,eqn3}, V, [u,v,w], {sin(b),cos(a),V});

                                 w      
                       V = -------------
                           sin(a) cos(b)

> G({eqn1,eqn2,eqn3}, sin(a), [u,v,w], {sin(b),cos(a),V});

                                   w    
                       sin(a) = --------
                                V cos(b)

> G({eqn1,eqn2,eqn3}, cos(b), [u,v,w], {sin(b),cos(a),V});

                                   w    
                       cos(b) = --------
                                V sin(a)

And one can se that there is still choice in this,

> G({eqn1,eqn2,eqn3}, V, [u,v,w], {sin(b),cos(a),V});

                                 w      
                       V = -------------
                           sin(a) cos(b)

> G({eqn1,eqn2,eqn3}, V, [u,v,w], {sin(a),cos(b),V});

                                 v   
                           V = ------
                               sin(b)

acer

The FromMatlab translator (emulator) has a few problems, but gets close.

I removed the initial transposition of tau and t when created as zero-Arrays. I added the plot call by hand, afterwards. (It could well be that, if the transposition had not had issues, then the plot call could also be translated. I'm not sure.)

Maple did obtain a similar graph as from Matlab.

Note that translated code is not always optimal. Since FromMatlab is better as an emulator, the Maple code equivalents that it produces are sometimes hard to read and not best to learn from. For this source sample, it's the code that creates tau and t which is least legible and (here) unnecessarily convoluted.

Matlab:-FromMatlab("
t_init      =20;      % [C] Start temp i mælken
h           =10;      % [W/m2*K] varme transmissionskoefficient
A           =0.07;    % [m2] Overfladeareal
m           =1;       % [kg] massen af mælken
c           =4190;    % [j/kgK] specifik varmekapacitet, mælk og karton
t_R         =5;       % [C] Temperatur i køleskabet
tau_slut    =3600*3;  % [s] Simuleringstid
dtau        =10;      % [s] Tids step

% Beregning - nogle hjælpestørrelser

i_max       = tau_slut/dtau; % Antal elementer (gennemløb i loop'et)

% Tre vektorer oprettes ud udfyldes med nuller

tau         = zeros(i_max,1); % Den numeriske løsning
t           = zeros(i_max,1); % Den numeriske løsning

% Så sættes nogle begyndelsesværdier

t(1)=t_init;
tau(1) = 0;       % Starttidspunktet

%Simulering
for i=2:i_max
    tau(i)=tau(i-1)+dtau;              % [s] Den nye tid
    dt=-(h*A)/(m*c)*(t(i-1)-t_R)*dtau; % [-] Temp. ændringen
    t(i)=t(i-1)+dt;                    % [-] Den nye temp.
end
                        "):

plot(<(tau/3600)|t>,gridlines=true,
                    legend="analytisk",legendstyle=[location=right],
                    caption="Mælketemperatur");

acer

eval(G(X1,X2,1),G=M1);

The `subs` command does syntactic substitution, without subsequent evaluation. It is often (mis)used in situations where `eval` is the more appropriate command. See both the first and the last paragraphs of the ?subs help-page's Description section.

acer

Regardless of which optimization command is used (eg. The DirectSearch or Optimization packages), the whole operation's timing performance will depend very much on the speed of each individual evaluation of the objective function.

On a mid-range Core2 Duo, the original code posted (as a follow-up comment in this thread) took about 13.5 seconds for a single objective function evaluation.

In the attached worksheet there are a variety of variations on re-coding the objective procedure. The best timing seems to be about 0.5 to 0.75 seconds on that same Core2 Duo.

I added some notes, hopefully illustrating in steps how some realizations and improvements came about, from variant to variant. (It needs the .xls data linked in that follow-up comment.)

F2proc_newer.mw

Note that there is one potential issue with the code that I ought to mention. In the original EIGEXP, the `res` Vector is created with size 1971, and 1971 elements are used from At and TC. But only 1969 elements of `res` are used to form the sum of squares. I presume that the OP knows what's going on with that, and can sort it out if necessary.

I only tested a smattering of individual points. A more thorough check of results between implementations could be done. Small variations past the 10th digit or so are expected, due to Digits, evalhf uses, etc.

acer

If you upload a worksheet containing the code then we might be able to speed up its objective function evaluations.

acer

> ee:=sin(x)^4/cos(x)^2;

                                  4
                            sin(x) 
                            -------
                                  2
                            cos(x) 

> eval(ee,cos=sin/tan);

                              2       2
                        sin(x)  tan(x) 

acer

Tell us just exactly how you'd do it by hand, and we could tell you how to do those steps in Maple.

(Alternatively, someone might just give you the answer in spiffy syntax, which won't help you much for the next time...)

acer

Something simple like this might do for you.

ASTEM97:=module()
option package;
export psat;
local ModuleLoad;
ModuleLoad := proc()
   local dllloc;
   dllloc := "c:/Windows/system/astemdll.dll";
   psat := define_external('psat97mc', 'FORTRAN',
                           'arg1'::('float[8]'),
                           'RETURN'::('float[8]'),
                           'LIB' = dllloc);
   # repeat for other exports   
   NULL;
end proc;
end module:


libname:="c:/TEMP",libname:
LibraryTools:-Create("c:/TEMP/astem97.mla");
savelib(ASTEM97);

The material above only needs doing once, each time you adjust and create the thing. Typically, you would save this build code in a worksheet or text file, so that you could rebuild it as you add more exports, etc. Note that the ModuleLoad routine only runs when the module is read from an .mla archive into a new session. So it won't run in this session, until you restart.

The next part below, after the restart, ie., the setting of libname, can be done in each new session prior to the with(). Or you can adjust the libname in your maple initialization file.

> restart:
> libname:="c:/TEMP",libname:
> with(ASTEM97);

Or you could move the created .mla file to a folder (new, you have to create it in your OS) such as what gets returned by the Maple command,

   cat(kernelopts(homedir),"\\maple\\toolbox\\14\\ASTEM97\\lib");

and then close the GUI and relaunch. Replace the "14" with your Maple version number.

Of course, since I don't have that .dll then all I get is,

> with(ASTEM97);
Error, external linking: error loading external library c:/Windows/system/astemdll.dll: The
specified module could not be found.

You can get fancy, wrapping the define_external calls in try..catch statements if you want it to do something special upon failure. But perhaps try the simple stuff first.

acer