Nice. Now I feel silly for not bothering to think first, about the nature of the equations. acer

Well, solve() should be able to handle (eq1,G=0}. So perhaps there's not need for optimizing in that case. I rather got the impression that the OP had some specific values of the parameters in mind, and that those might not be free to range. I thought that one might use Optimization:-Minimize, being careful with the abs, and either choosing a method not requiring derivatives or supplying them by hand. While GlobalOptimization might require only Lipschitz continuity, Optimization is pickier but sometimes one can coerce it. acer

Well, solve() should be able to handle (eq1,G=0}. So perhaps there's not need for optimizing in that case. I rather got the impression that the OP had some specific values of the parameters in mind, and that those might not be free to range. I thought that one might use Optimization:-Minimize, being careful with the abs, and either choosing a method not requiring derivatives or supplying them by hand. While GlobalOptimization might require only Lipschitz continuity, Optimization is pickier but sometimes one can coerce it. acer

The difference between linalg and LinearAlgebra might be explained by their using differing criteria for choosing between valid candidate pivots. The routine linalg:-ffgausselim selects according to length, ie, length(B[p[j],k]) < length(B[p[i],k]) The routine LinearAlgebra:-LA_Main:-LUDecomposition appears to select the first nonzero entry it can find (in the current leading column). There are reasons to prefer the selection by shortest length. As for the method itself, you can find it implemented in several places in the Maple library. In some places it is implemented as a means to compute the inverse or determinant, eg, LinearAlgebra:-LA_Main:-``MatrixInverse/polynom`. In other places it is implemented to compute a result with a modular technique, or over some particular field, eg. LinearAlgebra:-LA_Main:-`Determinant/modular` or LinearAlgebra:-Generic:-BareissAlgorithm. You can view the source of the Maple library routines named above by calling eval() around their names. You may first need to issue these two commands, kernelopts(opaquemodules=false): interface(verboseproc=3): These sources are mostly straightforward and the core algorithmic implementations within them are not too difficult to understand. They each usually differ in some way, according to their purpose, domain, etc. acer

The difference between linalg and LinearAlgebra might be explained by their using differing criteria for choosing between valid candidate pivots. The routine linalg:-ffgausselim selects according to length, ie, length(B[p[j],k]) < length(B[p[i],k]) The routine LinearAlgebra:-LA_Main:-LUDecomposition appears to select the first nonzero entry it can find (in the current leading column). There are reasons to prefer the selection by shortest length. As for the method itself, you can find it implemented in several places in the Maple library. In some places it is implemented as a means to compute the inverse or determinant, eg, LinearAlgebra:-LA_Main:-``MatrixInverse/polynom`. In other places it is implemented to compute a result with a modular technique, or over some particular field, eg. LinearAlgebra:-LA_Main:-`Determinant/modular` or LinearAlgebra:-Generic:-BareissAlgorithm. You can view the source of the Maple library routines named above by calling eval() around their names. You may first need to issue these two commands, kernelopts(opaquemodules=false): interface(verboseproc=3): These sources are mostly straightforward and the core algorithmic implementations within them are not too difficult to understand. They each usually differ in some way, according to their purpose, domain, etc. acer

Here's a stab at the idea posted above. By altering it so that FFFF_3 accepts A as a second argument (or uses it as a global) then A could be created outside of FFFF_3. This allows one to check that the compiled and uncompiled instances produce the same Y results at a small size like g=10 say.

FFFF_3 := proc(g)
global adder;
local i, Y, A;
  # Generate A with enough entries so that it is unlikely that
  # adder() attempts to access more than that number of values.
  A:=Statistics:-Sample(Statistics:-RandomVariable(':-Poisson'(1)),g*4);
  Y := Vector[row](g,datatype=integer[4]);
  adder(A,Y,g);
  [seq(Y[i], i = 1 .. g)];
end proc:
                                                                                                                                           
addproc := proc(A::Vector(datatype=float[8]),Y::Vector(datatype=integer[4]),g::integer)
local count::integer,j::integer,x::integer;
  count := 1:
  for j from 1 to g do
    x := trunc(A[count]);
    count:=count+1;
    if x=0 then
      Y[j]:=0;
    else
      Y[j]:=trunc(add(A[i],i=count..count+x-1));
    end if;
    count:=count+x;
  end do;
end proc:
                                                                                                                                           
adder := addproc:
T := time(): FFFF_3(2000): time()-T;
                                                                                                                                           
adder:=Compiler:-Compile(addproc):
gc():
                                                                                                                                           
T := time(): FFFF_3(2000): time()-T;

I don't know what the Compiler does with calls to the Maple built-in function add(). Maybe it makes the compiled procedure call back to the Maple kernel for such. But I did not get a big further improvement above by replacing the add() call in addproc() with, say,

      temp := 0.0;
      for i from count to count+x-1 do
        temp := temp+A[i];
      od;
      Y[j]:=trunc(temp);

acer

The (maple kernel) built-in routine trunc() may be faster than floor() for you, where you know that all your values will be greater than zero. But the improvement due to using trunc() instead of floor() may be small. There may also be some other, faster, ways to approach the whole problem. But these too might have their own bottlenecks to figure out. Here's an example idea, which I haven't tested. My guess is that done well it might perform faster, even if a first crack at trying it doesn't seem better. Suppose that, right up front and just the once, you created a very long 1-D Array or Vector from the distribution. Say with g * 10 entries. This is by calling Sample just once, like, A := Statistics:-Sample(Statistics:-RandomVariable((':-Poisson')(1)), g*10); Then you could create some new, inner, proc which can walk A. As it walks A, it keeps a counter. It take a value from A, increments the counter by 1. Take the floor of that, and let that number be assigned to C. Then walk and add up the next C entries, and take the floor and assign to y[i] or do what you will with it. Increment the counter by C. And repeat. The point is that, when calling the routine S returned by Sample(), a new Vector for the results is called each time. These Vectors make Maple garbage, and garbage collection is part of what's taking up time. Such time might be avoided by creating just a single Vector and incurring the Vector creation overhead just once. You can choose the length of Array A so that it is very, very likely to contain enough values from the distribution for your full purpose. You might be able to use add() to add up the C entries from current-counter to current-counter+C. Or you might be able to use for-loops and perhaps even use Maple's Compiler on this new inner procedure. acer

For a Matrix whose entries are all explicitly stored, one practical limit is the amount of available memory. If the amount of available memory is the effective limit, then using a hardware datatype (integer[1], float[8], etc) can allow one to create and use larger Matrices and Vectors. If, on the other hand, not all entries must be stored, then a hard limit in each dimension may be something like 2^(kernelopts(wordsize)-1)-1. This can be relevant when the Matrix is set up to have either so-called empty storage with unstored entries specified only by an indexing function, or sparse storage with a only a specific number of entries' storage allocated up front. For example, on 64bit Linux, M:=Matrix(2^63-1,2^63-1,storage=sparse[1]); and on 32bit Windows, M:=Matrix(2^31-1,2^31-1,storage=sparse[1]); One way to conceptually extend even that limitation could be to use multi-dimensional Arrays instead of Matrices or Vectors. An Array can have from 0 to 63 dimensions (as opposed to Matrices which have just 2 and Vectors which have just 1). Another possible extension might be to have Matrices whose entries are themselves Matrices. Eg, M:=Matrix(2^31-1,2^31-1,storage=sparse[1]); M[1,1]:=Matrix(2^31-1,2^31-1,storage=sparse[1]); acer

Suppose that the variable `sol` were assigned something like this, sol:= fsolve({Eq16, Eq15}, {beta = 4 .. 5, N_bar = 0 .. 1000}); Then, 2-argument eval() should do the trick, to get at the separate values. sol := {N_bar = 18.04394737 + 0.*I, beta = 4.474743382}; eval(N_bar,sol); eval(beta,sol); acer

Suppose that the variable `sol` were assigned something like this, sol:= fsolve({Eq16, Eq15}, {beta = 4 .. 5, N_bar = 0 .. 1000}); Then, 2-argument eval() should do the trick, to get at the separate values. sol := {N_bar = 18.04394737 + 0.*I, beta = 4.474743382}; eval(N_bar,sol); eval(beta,sol); acer

Don't call Sample() and RandomVariable() each time through the loop. Also, when forming x[i], since you're apparently only taking a single (1) value from the distribution then you don't need to call the heavyweight convert(..,`+`) on it. FFFF_2 := proc(g) local i, x, y, S, Samp; Samp:=Statistics:-Sample(Statistics:-RandomVariable(':-Poisson'(1))); for i to g do x[i] := floor(Samp(1)[1]); if 0 < x[i] then y[i] := floor(convert(Samp(x[i]),`+`)); else y[i] := 0 end if end do; S := [seq(y[i], i = 1 .. g)]; end proc; Also, why use local x as a table, when a scalar x would suffice (since you don't use the x[i] for anything else? If that's the case, then replace x[i] with simply x, and save some more. Lastly, don't load the Statistics package outside the proc definition. It's not good programming style. Either explicitly write out the commands, in the proc, as I've done or utilize the `use` functionality. acer

Can LPSolve be used for the entire problem? One might look at the problem as one of finding merely a feasible point. A linear programming problem with only equality constraints. Using a constant objective should allow it to return the first feasible point found. Optimization:-LPSolve(Vector(4),[NoUserValue,NoUserValue,Matrix([[1,2,7,5], [1,1,2,1]]),Vector(2,[3,1])],assume=nonnegative); If so, then could this help with some of the potential numerical difficulties? The help-page ?Optimization,LPSolveMatrixForm shows that this routine has extra options such as feasibilitytolerance and integertolerance. acer

Can LPSolve be used for the entire problem? One might look at the problem as one of finding merely a feasible point. A linear programming problem with only equality constraints. Using a constant objective should allow it to return the first feasible point found. Optimization:-LPSolve(Vector(4),[NoUserValue,NoUserValue,Matrix([[1,2,7,5], [1,1,2,1]]),Vector(2,[3,1])],assume=nonnegative); If so, then could this help with some of the potential numerical difficulties? The help-page ?Optimization,LPSolveMatrixForm shows that this routine has extra options such as feasibilitytolerance and integertolerance. acer

I get amusement from taking these sorts of examples as linear programming problems with a constant objective. The goal is then to find a feasible point. eq1:= q + d = 8: eq2:= 0.25*q + 0.10*d = 1.25: Optimization[LPSolve](1,{eq1,eq2},assume=nonnegint); More seriously, I don't think that one should shy away from suggesting use of isolve() for such examples. The fact that only whole coins may be possible is a genuine part of the problem. Consider cases where more than one answer is possible. For them, the result from isolve({eq2}) is better than that from solve({eq2}). acer

There's no object "sqrt(b^2)". Instead, the call sqrt(b^2) passes b^2 to a somewhat smart routine. The output is the object (b^2)^(1/2). So this is the sort of object at hand, this (b^2)^(1/2) thing. It's a structure, or an expression. It's a power thing. It's not a function, or an unevaluated function call. If you evaluate it at a specific b, then doing so won't do anything clever. It won't call sqrt() again. On the other hand, sqrt() is a function. And it is somewhat clever. So sqrt(2^2) produces 2 right away. Notice that (b^2)^(1/2) and sqrt(b^2) are not the same thing. The first is a structure, or expression. The second is a function call, which happens to return the first (under no assumptions on b). The sqrt() function is not the only mechanism that can construct (b^2)^(1/2). It could be constructed directly, typing it in just as it is. (That's a partial rationale for why evaluating the expression at a particular b doesn't call sqrt() again.) There's nothing clever in (b^2)^(1/2), that can take advantage of b>0 say. The cleverness resided in sqrt(), the function. Evaluating (b^2)^(1/2) at b= will not call sqrt(^2) once more. Once sqrt() returns its result, the opportunity for sqrt() to be clever has gone. In fact, you don't need to create the "special procedure " that you mentioned, because sqrt() is such a thing already. The issue that you are seeing is also present with your own procedure, a:=proc(b); return sqrt(b^2); end proc; To see that, just try, a(b); eval(%,b=2); Just like for sqrt(), once your a() gets called it returns an expression. Evaluating the expression doesn't cause anything clever to re-evaluate or simplify it. That's why a suggestion to hit it with a big hammer like simplify() was made. Also, there is no automatic simplification of ^(1/2) for this positive sort of . Depending on whom you ask, you may get different opinions about how much merit there is in that. A slightly smaller hammer is normal(). That is, 4^(1/2); normal(%); acer